Re: [boost] Re: ublas: hermitean matrices

From: boost (boost_at_[hidden])
Date: 2001-12-04 06:51:29

Salut,

On Tuesday 04 December 2001 10:12, walter_at_[hidden] wrote:

> > In my old code, for HMatrix< complex<double> > ,
> > I had a complex array for the off-diagonal elements,
> > and a  double* for the diagonal part. You may call it a premature
> > optimization,
>
> Of course ;-). Shouldn't our storage layout be Fortran BLAS
> compatible?
Currently I 'm afraid it's still important to provide FORTRAN BLAS
compatible layouts, since sometimes one has to call FORTRAN
linear algebra routines. But that shouldn't prevent the
implementation of improved layouts.

> Please, don't mix assignment and computed assignment. I currently see
> the following cases:
>
> Computed assign:
>
> hermitean_matrix *= double_scalar, no problem
> hermitean_matrix *= complex_scalar, precondition: assert (imag
> (complex_scalar) == 0)
>
> Assign:
>
> hermitean_matrix2 = double_scalar * hermitean_matrix1, no problem
> hermitean_matrix2 = complex_scalar * hermitean_matrix1,
> postcondition: assert (imag (hermitean_matrix2 (i, i)) == 0)
> general_matrix = complex_scalar * hermitean_matrix, no problem

But then you have run-time assertions, no complile-time checks.
In addition I'm sure that the case

postcondition: assert (imag (hermitean_matrix2 (i, i)) == 0),
but imag (complex_scalar) != 0,

will be have problems with truncation errors. 

In addition, the pre and post-condition are equivalent, if at least
one diagonal element of the hermitean matrix H is non zero.
Suppose:  complex_scalar c == a + i*b; where a,b are real:

A = c * H = a*H + i*b*H 
imag( A_jj ) = b* H_jj

> general_matrix = complex_scalar * hermitean_matrix, no problem
Maybe this should be the default, while
> hermitean_matrix2 = double_scalar * hermitean_matrix1,
is a special case ?

> Sorry, documentation is slightly outdated and incomplete. Therefore I
> mentioned hermitean_adaptor. For all packed matrix types currently
> supported (banded_matrix, triangular_matrix, symmetric_matrix,
> hermitean_matrix) we provide a proxy class (banded_adaptor,
> triangular_adaptor, symmetric_adaptor, hermitean_adaptor), which
> allows to address the corresponding view of a referenced matrix. So
> if you use hermitean_adaptor<banded_matrix<> >, you already should be
> able to work with an hermitean banded matrix stored in the upper or
> lower triangular part of a banded matrix.

I need a few more days before I can start playing with these matrices.

Best wishes,
Peter