Sorry, an error occurred while loading the content.

• ## Re: [ublas-dev] Re: [boost] Re: Boost Mathematicians

(15)
• NextPrevious
• Hi Deane, I will answer just to ublas-dev as this is getting a bit detailed. I not sure how to express this in a very clear fashion. A block diagonal matrix is
Message 1 of 15 , Jul 6, 2004
View Source
• 0 Attachment
Hi Deane,

I will answer just to ublas-dev as this is getting a bit detailed.

I not sure how to express this in a very clear fashion. A block diagonal
matrix is a partitioned matrix not a matrix where the elements are matrices.

A block or partitioned matrix is a "special matrix" in the same way as a
triangular matrix is special. We have to be careful not rely to much on the
analogy with a matrix whose element types is a matrix.
I define it using a matrix shape property like other special matrices.

Of course programatically we can use a matrix of matrices. As long as we just
want the data stored or we use addition and subtraction things works fine. As
a matrix in Linear Algebra we have to be careful as many of the necessary
identities do not hold.

In general solving with such a beast will not give the result we expect. Block
diagonal is the only case were the inverse is the same because the cross
terms are all zero.

Does that make sense?
Michael
--
___________________________________
Michael Stevens Systems Engineering

Navigation Systems, Estimation and
Bayesian Filtering
http://bayesclasses.sf.net
___________________________________
• Hi Peter, ... Yes... the term matrix has a multitude of meanings, all of which are confusingly similarly but slightly different. However, the current topic
Message 2 of 15 , Jul 6, 2004
View Source
• 0 Attachment
Hi Peter,

On Jul 6, 2004, at 2:51 AM, Peter Schmitteckert (boost) wrote:
> I agree with you, if you restrict matrices to linear algebra. That's a
> fine
> definition. But I'm used to use the notion 'matrix' in a more general
> sense.

Yes... the term "matrix" has a multitude of meanings, all of which
are confusingly similarly but slightly different. However, the current
topic
of conversation is a linear algebra library. Thus it makes sense for
the purpose
of this conversation to restrict the meaning to that of traditional
linear algebra.

Believe me, I'm not trying to exclude physicists or anyone. I'm highly
in
favor of there being lots of numerical libraries in Boost. There can
even
be many different concepts with the name "Matrix", they'll just be in
different namespaces. All I'm saying is that the concept named "Matrix"
in the linear algebra namespace should have the standard linear algebra
meaning and syntax.

Cheers,
Jeremy

_______________________________________________
Jeremy Siek <jsiek@...>
http://www.osl.iu.edu/~jsiek
Ph.D. Student, Indiana University Bloomington
C++ Booster (http://www.boost.org)
_______________________________________________
• Hi Michael, ... 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
Message 3 of 15 , Jul 6, 2004
View Source
• 0 Attachment
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@...>
http://www.osl.iu.edu/~jsiek
Ph.D. Student, Indiana University Bloomington
C++ Booster (http://www.boost.org)
_______________________________________________
• Hi Jeremy, ... Maybe you could convince Dave to do the same for uBLAS! Michael
Message 4 of 15 , Jul 7, 2004
View Source
• 0 Attachment
Hi Jeremy,

> Then it's a good thing we have Dave to try to implement it :)
Maybe you could convince Dave to do the same for uBLAS!

Michael
Your message has been successfully submitted and would be delivered to recipients shortly.