[matrix] boost.numeric (matrix library)

From: boost (boost_at_[hidden])
Date: 2001-03-28 18:31:01

Dear  J. Walter and M. Koch & boost'ers,

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 ... ). 

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. This is nice if you have to calculate
sin(x) or atan(x) to calculate rotation matrices. LA_HermiteanMul was 
important
in a few routines special to my problems.

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.

Same applies for the type of the eigen values,
which is also double for HMatrix< complex<double> >.

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.

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.

Best wishes,
Peter

Thanks for reading through this long mail