[Lapackpp-devel] Lapack++ column ordering, derivations, exceptions, lwork size

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

Dear Jack,

this is great. Thanks for your contributions and thoughts. I am looking 
forward to see your code in CVS soon :-)

We should continue this discussion on the mailing list lapackpp-devel -- 
this has at least the additional advantage that we can automatically use 
sourceforge's mailing list archive. You can subscribe on 
http://lists.sourceforge.net/lists/listinfo/lapackpp-devel

Jacob (Jack) Gryn schrieb:
> I implemented it yesterday with the enum-style method; just changed it to
> use a bool now.  Works fine.  
> 
> I only did this for LaGenMatDouble's I'll do the other data types as well
> and post my changes on CVS shortly.

Sure. Either way -- you can commit it only for LaGenMatDouble for now, 
or you can also commit it for all four important classes at one.

If you commit to CVS, please have a look at the (newly added) file 
README.cvs. It basically emphasizes to add a ChangeLog entry for every 
commit. Apart from that, you are free to commit your code -- I will be 
notified of every commit via the lapackpp-cvs mailing list (which you 
can subscribe as well, if you like to), so I will check every commit 
anyway and might notify if there are any problems. In short, I happily 
encourage you to go ahead and we'll sort things out early enough :-))

> Ideally, there should be better interoptability between the matrix/vector
> data types as there's a lot of duplicate code, and no inheritance.
> For example, it would probably be better to treat a vector as an n x 1
> matrix, inheriting all the properties of a matrix, except that if has a
> coordinate fixed at 1.  This is probably a longer-term thing (maybe for
> 3.0), but since the point is to use C++, why not take advantage of its
> features.  Then, code that, for example, re-orders a matrix upon loading,
> doesn't need to be written 6 times, for each data type.

The LaVectorDouble is a derived class of LaGenMatDouble, isn't it? And 
for LaVectorComplex / LaGenMatComplex as well? But the constructors 
cannot be derived, and, yes, there is some code duplication because of 
this. The re-ordering code should indeed not be rewritten more than 
necessary, and it is by definition only necessary in the matrix classes, 
not in the vector classes. Therefore I thought it is necessary exactly 
twice: In LaGenMatDouble for doubles and in LaGenMatComplex for COMPLEXs.

> Another thing I was wondering, why is it necessary to #define
> LA_COMPLEX_SUPPORT in files where no complex numbers are used?

This define-macro comes from the old v1.1 lapack++ where it was 
*disabled* everywhere. Maybe it should be dropped altogether, so that 
lapackpp will be built with complex support enabled always. If not, then 
this macro should be the switch to turn on or off the declarations of 
the complex data types. In that sense the #define should not be needed 
in files that do not use complex numbers.

 >>>I did some testing with the current developer CVS version..
 >>>I still get "Aborted" if the vector S is of an incorrect size.
 >
 >> Right. The code right now will only accept a vector S of the exactly
 >> correct size. Do you suggest to change this behaviour? I think this
 >> behaviour is ok.
 >
 > I think it's a good idea to check the vector size.  My issue is that
 > the word "Aborted" is kind of useless, when the
 > LaException::enablePrint(true)
 > flag is set.  Maybe I'm just picky.

The word "Aborted" comes from the fact that the library has thrown an 
exception but the application did not catch it. This is the default 
behaviour of C++ (more precisely: gcc) with uncaught exceptions, and we 
cannot do anything about it. Either the application is supposed to have 
a try-catch-block for LaException, or the user has to use a debugger 
anyway which will then show the origina of the exception in its backtrace.

 > I tool a look at the lwork value, according to the documentation at
 > http://netlib2.cs.utk.edu/lapack/patch/src/dgesdd.f
 >
 > LWORK >= 4*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N).
 > (Larger for more performance)
 >
 > I set it to that in both functions that call DGESDD, and it seems to
 > have fixed the problem!  I'm not sure how much higher it should be set
 > to if we want more performance; I guess its ok for now.
 >
 > I have only found the problem with SVDs of LaMatGenDoubles.

If this modified lwork definition fixes your problem, then go ahead and 
commit it. According to the documentation, it is possible to get the 
optimal size of lwork calculated by DGESDD itself (by an additional call 
to DGESDD beforehand). In this case this would be the optimal solution, 
but to be honest, I'm too lazy to implement this :-) If you want to 
commit something like that, just go ahead.

Christian