Re: [ublas] diagonal matrix and scaled axpy

From: Michael Stevens (mail_at_[hidden])
Date: 2005-04-26 09:18:43

Hello Gunter,

On Monday 25 April 2005 14:29, Gunter Winkler wrote:
> Hello,
>
> I was looking for an ublas expression, to compute
>
> y += D * A * x;
>
> where D is a diagonal matrix (stored in a vector of the diagonal elements
> of D) and A is a sparse matrix.

I have been wishing a vector to diagonal matrix adaptor too. It is something I 
need quite regularly.

> In banded.hpp I found the (still to document) diagonal_matrix and
> diagonal_adaptor types which are special banded matrices (with lower and
> upper bandwidth equal to zero).

I never noticed they were undocumented. Definitely something to be fixed.

> Using these types with prod(D,A) is not optimal, since it uses the most
> complicated available prod: (sparse = sparse * sparse).

Are you sure prod dispatches to sparse*sparse. It should be dispatching to 
a 'packed' prod?

> Should we introduce a new matrix proxy? Something like
> DA := scaled_matrix(D, A)
> where D is a vector and A a matrix. DA(i,j) would evaluate to D(i)*A(i,j)
> and could reuse all iterators of A. The proxy should be usable (and
> optimizeable) by axpy_prod, too.

Interesting solution. I'm wondering if it is not possible to define prod so
prod (diagonal_adaptor(vector), matrix) result in such a scaled_matrix object 
automatically.

Michael
-- 
___________________________________
Michael Stevens Systems Engineering
34128 Kassel, Germany
Phone/Fax: [phone]
Navigation Systems, Estimation  and
                 Bayesian Filtering
    http://bayesclasses.sf.net
___________________________________