Re: [ublas-dev] Re: [boost] Re: Boost Mathematicians
Hi Michael, On Jul 5, 2004, at 4:52 PM, Michael Stevens wrote:
Hi Jeremy,
Thanks for chipping in on this. I don't think we have conversed before, apart from maybe a few patches I submitted for MTL. Although it would be great for everyone to work on uBLAS I think a revived MTL would also be a good thing. The design of MTL was very influential.
Thanks!
I won't comment on the other issues as I have already posted a few thoughts on these matters.
Mmmh, this open pitfalls concerning A*B*C in performance issues for non experts.
No, the right way to handle this is to have the expression template machinery identify that a temporary is needed, and create one automatically.
I wonder if you caught the discussion on this about 6 months ago. I would strongly argue that the automatic creation of temporaries in the subexpression context is a bad thing...
The roles of the expressions we declare (as in declarative programming) is twofold. a) to express the mathematical context b) to express how this should be computed. That is, to be able to express how the expression is to be evaluated. This gives us control over its numerical properties and also of the time and space complexity.
In principle, uBLAS has a one to one relationship between an expression and how it is evaluated. Simply stated this is based on an element by element evaluation of assignment. There are shortcuts in the case of sparse types but the principle holds.
An expression A+B+C can be evaluated in this manner trivially. A*B*C can also be evaluated in this manner, although the element by element rule results in an O() that is not what many users expect. To get the O() that users expect we could automagically insert temporaries. I fear this is a tar pit however. What temporary should be chosen? Each of A B or C could be a either dense, sparse, compressed, banded, symmetric or combinations therefore. The type the user expects or is "the best" is not at all clear.
Although automagically inserting temporaries in non-trivial (especially since it must all be done with template meta-programming) I do not think it is a tar pit. Here is the simply heuristic that I used for MTL 3 to decide where to put temporaries: 1. If an operation needs to read from an argument multiple times, evaluate the argument eagerly into a temporary. For example, a sparse matrix in compressed row format times the sum of two vectors: A * (x + y) In this case, the (x + y) would be evaluated into a temporary before evaluating the multiplication. 2. If an operation writes into its output in multiple passes, then perform the operation in an eager fashion. For example, the sum of a vector and the product of a sparse matrix in compressed column format and a vector. x + (A * y) In this case, the (A * y) would be evaluated into a temporary before evaluating the addition.
I think the only correct solution is to force the user to specify this type. Incidently, the no temporary evaluation should also be an option as in some cases (small sizes) it is desirable. uBLAS currently disallows this to avoid confusing the user. I should be reinstated in some form however.
All the best, Michael
-- ___________________________________ Michael Stevens Systems Engineering
34128 Kassel, Germany Phone/Fax: +49 561 5218038
Navigation Systems, Estimation and Bayesian Filtering http://bayesclasses.sf.net ___________________________________
_______________________________________________ Jeremy Siek <jsiek@osl.iu.edu> http://www.osl.iu.edu/~jsiek Ph.D. Candidate, Indiana University Bloomington C++ Booster (http://www.boost.org) _______________________________________________
Hi Jeremy, On Tuesday 06 July 2004 02:51, Jeremy Graham Siek wrote:
Although automagically inserting temporaries in non-trivial (especially since it must all be done with template meta-programming) I do not think it is a tar pit. Here is the simply heuristic that I used for MTL 3 to decide where to put temporaries:
1. If an operation needs to read from an argument multiple times, evaluate the argument eagerly into a temporary. For example, a sparse matrix in compressed row format times the sum of two vectors:
A * (x + y)
In this case, the (x + y) would be evaluated into a temporary before evaluating the multiplication.
2. If an operation writes into its output in multiple passes, then perform the operation in an eager fashion. For example, the sum of a vector and the product of a sparse matrix in compressed column format and a vector.
x + (A * y)
In this case, the (A * y) would be evaluated into a temporary before evaluating the addition.
OK, this 'eager' evaluation logic seems sound to me. I assume it operates in a bottom up fashion on the expressions replacing operators with temporaries in either of the two cases. Therefore providing a good heuristic to choose when a temporary seems possible. My problem is the exact choice of type for any temporary.... The mixing of sparse and dense is problematic. However your second example shows that we can get a hint from the types involved in many cases. But can this be done in the general case? Let me make problem a little more generic. Even with uBLAS's weakly defined storage concept it is possible to parameterise A or y so it storage is off-line (on file). Although we know a temporary should be used we now need to make a choice of storage type. In principle this is solvable (type traits, type promotion etc), particularly if storage type could be made orthogonal to other type parameters (dense/sparse, matrix/vector). This is however the generic "tar pit" I forsee. Maybe I am being too sceptical in think that if the is no simple (for user to understand even if not to implement) generic solution we should not solve the base cases. I guess I am also sceptical as without a deeper analysis of operators and types the user cannot know when temporaries will automagically appear. One big disadvantages of old OO C++ matrix libraries was the complete lack of control over temporaries. Simply stated I favour a strong correspondence between syntax and semantics. All the best, Michael -- ___________________________________ Michael Stevens Systems Engineering Navigation Systems, Estimation and Bayesian Filtering http://bayesclasses.sf.net ___________________________________
Hi Michael, On Jul 6, 2004, at 4:44 AM, Michael Stevens wrote:
Hi Jeremy,
On Tuesday 06 July 2004 02:51, Jeremy Graham Siek wrote:
Although automagically inserting temporaries in non-trivial
(especially
since it must all be done with template meta-programming) I do not think it is a tar pit. Here is the simply heuristic that I used for MTL 3 to decide where to put temporaries:
1. If an operation needs to read from an argument multiple times, evaluate the argument eagerly into a temporary. For example, a sparse matrix in compressed row format times the sum of two vectors:
A * (x + y)
In this case, the (x + y) would be evaluated into a temporary before evaluating the multiplication.
2. If an operation writes into its output in multiple passes, then perform the operation in an eager fashion. For example, the sum of a vector and the product of a sparse matrix in compressed
column > format and a vector. > > x + (A * y) > > In this case, the (A * y) would be evaluated into a temporary before > evaluating the addition.
OK, this 'eager' evaluation logic seems sound to me. I assume it operates in a bottom up fashion on the expressions replacing operators with temporaries in either of the two cases.
Yes, though there is a least one case when information from above is needed. For example, the outer product of two vectors can be implemented in either a row-wise or column-wise fashion (with respect to the output matrix). So, for example, in A = outer_product(x, y) you need to know the orientation of A to decide how to evaluate the outer product.
Therefore providing a good heuristic to choose when a temporary seems possible. My problem is the exact choice of type for any temporary....
The mixing of sparse and dense is problematic. However your second example shows that we can get a hint from the types involved in many cases. But can this be done in the general case?
I believe so. There are several classes of operations that combine sparsity in different ways. 1) The operation "ands" the sparsity. E.g., addition, subtraction. 2) The operation "ors" the sparsity. E.g., element-wise multiplication. So, you have a table like this: dense "and" dense -> dense dense "and" sparse -> dense sparse "and" dense -> dense sparse "and" sparse -> sparse
Let me make problem a little more generic. Even with uBLAS's weakly defined storage concept it is possible to parameterise A or y so it storage is off-line (on file). Although we know a temporary should be used we now need to make a choice of storage type.
That is indeed a problem, but we can provide a default choice and also provide for user-customization using a traits class. (as I see you mention next)
In principle this is solvable (type traits, type promotion etc), particularly if storage type could be made orthogonal to other type parameters (dense/sparse, matrix/vector). This is however the generic "tar pit" I forsee.
Then it's a good thing we have Dave to try to implement it :)
Maybe I am being too sceptical in think that if the is no simple (for user to understand even if not to implement) generic solution we should not solve the base cases.
I guess I am also sceptical as without a deeper analysis of operators and types the user cannot know when temporaries will automagically appear. One big disadvantages of old OO C++ matrix libraries was the complete lack of control over temporaries. Simply stated I favour a strong correspondence between syntax and semantics.
Yes, that is an important concern. However, I believe the above heuristics are simple enough to provide a strong correspondence. However, we won't really know for sure until we try it.
All the best, Michael
_______________________________________________ Jeremy Siek <jsiek@osl.iu.edu> http://www.osl.iu.edu/~jsiek Ph.D. Student, Indiana University Bloomington C++ Booster (http://www.boost.org) _______________________________________________
participants (3)
-
Jeremy Graham Siek -
Jeremy Siek
-
Michael Stevens