Re: ublas: hermitean matrices

From: walter_at_[hidden]
Date: 2001-12-04 17:55:25

--- In boost_at_y..., Peter Schmitteckert (boost) <boost_at_s...> wrote:
> Salut,
> 
> On Tuesday 04 December 2001 10:12, walter_at_g... 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.

Agreed. BTW, Fortran BLAS even uses something like 
hermitean_adaptor<matrix<> > :-).
 
> > 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.

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

That is conceivable. Do these border cases justify the introduction 
of a new feature? We then would have to invent something I'd like to 
call return type deduction,  i.e. implement and use a couple of 
traits classes, which describe for example, that the product of a 
hermitean with a real (or a complex with zero imaginary part?) is a 
hermitean and the product of a hermitean with a complex (with nonzero 
imaginary part?) is a general matrix (not to forget the corresponding 
promotion rules).
 
> 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

Theoretically yes, but I described, what we're currently able to 
detect via runtime checks.
 
> > 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 ?

It's the crux with our expression templates, that the right hand side 
is not a concrete matrix. The client code determines the matrix type 
to assign to. Possible solutions are runtime checks or 'return type 
deduction', IMO.
 
[snip]
 
Regards
 
Joerg