Hello Joerg,
Hi Patrick,
without a symmetric type I can just check for a diagonal matrix 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 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 // Diagonal matrix class
template // Construction and destruction
BOOST_UBLAS_INLINE
diagonal_matrix ():
matrix_type () {}
BOOST_UBLAS_INLINE
diagonal_matrix (std::size_t size1, std::size_t size2):
matrix_type (size1, size2) {}
BOOST_UBLAS_INLINE
~diagonal_matrix () {} // Assignment
BOOST_UBLAS_INLINE
diagonal_matrix &operator = (const diagonal_matrix &m) {
matrix_type::operator = (m);
return *this;
}
}; and the corresponding diagonal_adaptor to the code base? 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. Are there some recommendations for namespaces for functions like "pinv"
?
Inside the loop, shall I use size_type? I currently haven't enough time to look into a SVD implementation. So I'll
cc this reply to http://groups.yahoo.com/group/ublas-dev. The SVD is not important (i do not know what is inside either). Just keep in
mind that a SVD for a matrix<T>A produces 3 matrices:
matrix<T>U,diagonal_matrix<T>S, matrix<T>V. I should have pointed that out
more clearly. Sorry.
Aeh, is this thread OT here? Shall I just send these topics to the mailing
list uBlas-dev? I will make a reply to the forward as well.
Thank you, greetings
Patrick