Re: [matrix] boost.numeric (matrix library)

From: walter_at_[hidden]
Date: 2001-03-28 17:28:23

Hello, 

Peter Schmitteckert wrote:

> I had a short look at your matrix classes.
> Like many other I have also  written 
> my own expression template library for matrices.
> I had to do some linear algebra (eigenvalues,
> eigen vectors, linear system of equation, on dense and sparse
> matrices, dyadics vectors ... ). 

Many things to come :-).

> One problem was to write a generic routine for diagonalitzation of
> dense symmetric/hermitean matrices, e.g. HMatrix<T>
> In my implementation I introduced a helper class which defines 
> some functions on the type T,
> e.g.
>   double LA_Abs( const double& z)
>   {
>     return fabs(z);
>   }
>   double LA_HermiteanMul( const double& a,  const double& b)
>   {
>     return a*b;
>   }
> 
>   double LA_SquaredAbs( const  complex<double>& a)
>   {
>     return a.real() * a.real() + a.imag() * a.imag();
>   }
> 
>   double LA_HermiteanMul( const  complex<double>& a,  const  
complex<double>& 
> b)
>   // defined as ( a*b + b*a) /2
>   {
>     return  a.real() * b.real() - a.imag() * b.imag();
>   }

> Note that the return type of x = LA_SquaredAbs( const  
complex<double>& a) is
> a double and not a complex number. 

Yep. We call this one currently type_traits<std::complex<double> 
>::norm_type;

> LA_HermiteanMul was 
> important
> in a few routines special to my problems.

BTW, where is the typo in LA_HermiteanMul, comment or code?

> In the same way I defined Hermitean Matrices in a way that
> the Diagonal Elements for HMatrix< complex<double> > are
> of type double.   One has to be careful with 
>   complex<double>& HMatrix::operator(i,j)
> since this reference doesn't exist in my implementation,
> as the diagonal elements are stored in a double* array.
> But it works fine.

General complex matrices and corresponding operations are supported 
with BLAS-like functionality. If you immediately need a Hermitean 
specialization, have a look at our band_matrix (banded_matrix?) 
specialization. However to take full advantage of a Hermitean 
specialization, we first have to refine our iterator concept.

I am not sure, whether I like your idea to store the diagonal 
elements in a separate array (vector?).

> I would be interested in switching to boost and to do
> some contributions to a matrix/linear algebra library
> if Hermitean matrices will be supported.

IMHO Hermitean matrices should be supported.

> Please note, that I even have to use more abtract classes as
> complex numbers, e.g. small matrices which share the property
> that the norm type and the eigen value type differ from the
> matrix element type.

Even more things to come :-).

> Uups, I just found your 
> 
> template<class E, class F>
>     class vector_scalar_unary
> {
> ...
> }

These things need a lot of documentation, I fear.

Best regards

Joerg Walter