Re: [Lapackpp-devel] Error in LaGenMatDouble::copy(...) and other questions

[Christian S. <sti...@tu...>]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Dominik Wagenfuehr schrieb:
> 1. I havent' found the exact error yet, but if I debug my program I will
> get a segmentation fault in
> 
> LaGenMatDouble B(A);
> 
> or
> 
> LaGenMatDouble B;
> B.copy(A);
> 
> A is a (4 x 2)-matrix with some values I have calculated before.
> Unfortunately I cannot debug the lapack-libary so I don't know where the
> error could be. I try to extract the relevant code parts from my program
> to reconstruct it.

That's a serious error. I have spent a lot of checking on the
LaGenMatDouble methods, including the ::copy() method. If you can
extrace the relevant code portion than that would be really helpful.

> 2. Is it possible that you forgot to create the copy-function for
> LaSymmBandMatDouble? I only get
> 
> /usr/local/include/lapackpp/sybmd.h:51: warning: inline function
> 'LaSymmBandMatDouble& LaSymmBandMatDouble::copy(const
> LaSymmBandMatDouble&)' used but never defined
> 
> I have created my own copy-function yet but as I see you have used some
> templates for this. I don't understand much about templates but it would
> be great if you can fix it.

In fact the copy() method is declared "inline" but its definition is in
matrix/src/sybmd.cc, so it is not at all inline. You can fix this
yourself by removing the "inline" keyword in front of copy() in sybmd.h.
 Again, this isn't "my fault" :-) because this code has been untouched
from the original lapackpp authors, see the version history of sybmd.h
at
http://lapackpp.cvs.sourceforge.net/lapackpp/lapackpp/include/sybmd.h?view=log
- - I only started to work on that file after you asked for these issues
yesterday.

> 3. Is there such a simple thing like A := A + alpha*B ? I have only
> found Blas_Mat_Mat_Mult and actually I perform
> 
> Blas_Mat_Mat_Mult(B, I, A, alpha);
> 
> where I as the Identity Matrix. Especially for LaSymmBandMatDouble you
> could save a lot lot lot of time if there would be such a simple
> operation for different types of matrices. I don't know if there exists
> a BLAS-function for it (I haven't found it yet), 

The point with these "simple arithmetics" for the matrices is that
surprisingly few of these are actually available as BLAS functions.
Especially the direct addition of two matrices (with or without scalar
factor) is *not* available in BLAS or LAPACK, see
http://www.netlib.org/blas/blasqr.pdf or
http://www.netlib.org/blas/faq.html or
http://www.cs.colorado.edu/~jessup/lapack/aboutBlas.html

I don't know the reasons for this decision. I'm an engineer, not a
numerical mathematician. I can only *guess* that the plain addition of
matrices for whatever reason is never being computed explicitly in "real
numerics".

The goal of lapackpp, on the other hand, is to be a C++ interface to
BLAS and LAPACK. So if you're asking for functionality that isn't
available in there, then I'm hesitating to put it into lapackpp. Rather,
I'd like to ask you to implement this yourself - I just guess the
LAPACK/BLAS people do have their reasons for deciding what to implement
and what not to implement, and since they are more intelligent than me,
I'd rather stick to their decision. Sorry.

> 4. Same with the operation A = B^T * B. There seems to be no
> BLAS-function for it even if you know that A is symmetric. You could
> save half of the time. And further if you know that B has a special
> structure and A will be banded symmetric...

This function is surely available in BLAS (I've looked this up on the
Quick Reference Guide). It is called DSYRK there (see "man dsyrk" if you
have lapack-manpages installed). The correct place for the C++ function
declaration is include/blas3pp.h and the implementation (=definition) in
blaspp/src/blas3pp.cc, and it turns out this function is already
available as Blas_R1_Update(), search for LaSymmMatDouble and set
alpha=1, beta=0. I have no idea why that name was chosen - maybe I
should add a better function name like Blas_RK_Update().

For the reference, if you discover more such functions that *are*
available in BLAS but not yet in LAPACKPP, here's what we would need to do:

1. Find the BLAS function name (say, dsyrk in level-3-blas)

2. Check whether a declaration of that function name exists in
include/blas3.h (no C++ here, only plain C) - if it doesn't exist, add
the declaration by copy&paste a similar declaration and suitable
reordering of the arguments.

3. Invent a C++ function name (say, Blas_RK_Update) and add a suitable
declaration for that function into include/blas3pp.h.

4. Add the implementation for the C++ function to blaspp/src/blas3pp.cc;
follow the existing checks and argument conversions from the other
existing wrapper functions.

5. (slightly optional) Add some testing code into blaspp/test/tblas++.cc
that will make the tblas++.cc program fail if your new function doesn't
work as expected. This ensures that "make check" will always check your
newly added function.

> PS: I have tested the functions in sybmd.h and sybfd.h and they seem to
> work fine. (Ignoring the little copy-incident above!) I will try to
> comment the functions...

Good. Again I'd invite you to submit any code changes that you consider
an improvements. That is great. Thanks a lot.

Christian
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.2.1 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iQCVAwUBRVxnvGXAi+BfhivFAQJstAP/f9Vf6T/XwqyypyOkP5SgELWpB1YmC0Yd
Fdl672DGaWVExRowrmfWzaXIhED1vXPyU+4rd1tfF28dmAiHa0+1w5SDZVwNIQmm
WrI5lTNRZIwpDfncp0C3sVoDxXWAWhBLlVmwcmSytTxE/b1tf43daOKfAqPNFbDE
q+sMQbhPYO4=
=KT/y
-----END PGP SIGNATURE-----