you wrote:

> > > without a symmetric type I can just check for a diagonal matrix

Yep.

> inside the

> > > pinv function and branch, see code below.

>

> That in fact is not true. Defining the interface for a moore-penrose-

> inverse function could look more like this

> (answering myself again ;-) ):

>

> diagonal_matrix<T> pinv(const diagonal_matrix<T>&);

> matrix<T> pinv(const matrix<T>&);

>

> but

>

> banded_matrix<T> pinv(const banded_matrix<T>&)

>

> is not ok, because the pinv of a banded_matrix<T> is a rectangular

> matrix<T>. So branching inside the function ist not nice. At all I

> think it is a nice example, that in a mathematical kind of view,

> there could be a huge difference between banded and diagonal.

> > As far as I've understood you're looking for a type safe diagonal

But the difference in speed could be ;-) BTW: prod(prod()) is dangerous.

> variant of

> > banded_matrix. Would it be helpful, if we'd go to add something like

>

> it would be very nice to have something like that. But I think it is

> only useful if there are more applications with a semantic

> mathematical meaning. In my case the difference in the code is not

> huge:

>

> 1.) without diagonal_matrix:

>

> boost::numeric::ublas::matrix<double> pinv(const

> boost::numeric::ublas::matrix<double>& m)

> {

> using namespace boost::numeric::ublas;

> matrix<double> U,V;

> banded_matrix<double> S;

> svd(m,U,S,V); // assuming we have one

> for (size_t i=0;i<S.size1();i++) if (S(i,i)!=0) S(i,i)=1/S(i,i); //

> pseudoinverse of S stored in S

> return prod(prod(V,S),trans(U));

> }

>

> 2.) with diagonal_matrix:

>

> boost::numeric::ublas::diagonal_matrix<double> pinv

> (boost::numeric::ublas::diagonal_matrix<double>& S)

> {

> for (size_t i=0;i<S.size1();i++) if (S(i,i)!=0) S(i,i)=1/S(i,i);

> return S;

> }

>

> boost::numeric::ublas::matrix<double> pinv(const

> boost::numeric::ublas::matrix<double>& m)

> {

> using namespace boost::numeric::ublas;

> matrix<double> U,V;

> banded_matrix<double> S;

> svd(m,U,S,V);

> return prod(prod(V,pinv(S)),trans(U));

> }

Recommended usage is either prod(prod<some_matrix_type>()) or separating the

products.

> while in general I like more the second one, but the first one is

Yes.

> more efficient :-).

> So, in this case we talk about one line moved out

OK.

> and afterwards moved in again to be more effective.

> OK, to be more enthusiastic: I think a diagonal_matrix would be very

> nice, it would be useful for me anyway.

[snip]

> > and the corresponding diagonal_adaptor to the code base?

Oops. Fixed, thanks.

>

> YES !!!! JUHU !!!

> ...but I think

> > BOOST_UBLAS_INLINE

> > diagonal_matrix (std::size_t size1, std::size_t size2):

> > matrix_type (size1, size2) {}

> could be:

> BOOST_UBLAS_INLINE

> diagonal_matrix (std::size_t size):

> matrix_type (size,size) {}

> or both.

Best,

Joerg