[Lapackpp-devel] Re: lapack++ conversion

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

Dear Walter,

thank you for your interesting question. Before I answer this, let me explain 
why lapack++ is in the state that it currently is: I'm not a math expert. I 
don't know much about numerics, or at least not more than what we learnt in 
the "math for engineering" lectures. I am a researcher in communication, and 
I am also an application programmer in C and C++ (but I have never learnt or 
programmed in Fortran). But I took over lapackpp because I needed a library 
for complex-valued matrices in C++, and the existing lapack++ code seemd 
"almost" usable. Eventually I managed to get it to do those complex-valued 
matrix operations that I needed, but everything else was kept pretty much the 
same and that is why the library looks how it currently looks... 

Am Mittwoch, 26. April 2006 19:33 schrieb Walter Schreppers:
> I was wondering how the lapack++ conversion from the original fortran
> lapack is being done.
> Are the wrapper functions hand written or are their tools involved
> (like swig) which automatically generate these?

Are you referring to the data type conversions, or the function call 
conversions? 

The data types are accessed as follows: Any matrix is passed as a plain C 
array into the fortran function call, i. e. a "double" matrix is passed as a 
"double*" pointer that points to the start of the array. Same with all other 
matrix data types. This uses the fact that fortran and C "happen" to use the 
same storage format for all basic number types. (Complex numbers are the 
notable exception from this rule, though.)

The function call conversions are actually all using hand-written wrappers. 
The C declarations of the fortran-lapack functions are in include/lapack*.h. 
The call of the correct fortran-lapack function then happens in the various 
C++ functions that accept the lapack++ matrix types as arguments, e.g. in 
blaspp/src/blas1pp.cc. Those again are all hand-written. I don't think either 
one could be created automatically -- there is too much weird argument 
passing involved. 

The C declarations of the fortran-lapack functions originally have been 
created by f2c. But some of them seemed wrong and were subsequently 
hand-corrected by myself according to the required arguments to the fortran 
functions that were documented in the lapack man pages.

The current structure of lapack++ leads to the unfortunate result that for 
each desired fortran-lapack operation, I need to (1) write or retrieve from 
somewhere the C declaration, and then (2) write a C++ function that accept 
the lapack++ matrix as argument. And both steps for each and every number 
data type, if this is needed. Quite a lot of hand work. 

I thought about template functions here to save some typing (e.g. to replace 
step (2)), but actually it doesn't seem possible to replace any of this by 
templates. The reason is it would require "template specializations", and 
those are only allowed for class templates but *not* for function templates 
in C++. It would require template specializations, because eventually, for 
each number type, a different fortran-lapack function needs to be called. In 
order to properly distinguish this, one would need a separate specialized 
template for each number type where the correct fortran-lapack function will 
be called for that number type. But those function template specializations 
are not allowed in C++ (the reason for that being some ambiguity issues, if I 
remember correctly). There would be the option to put the call to lapack 
operations in a class... but the original intention of lapack++ was a direct 
representation of lapack, and the whole lapack concept is about operations 
that are called as functions. (Actually thinking about it, there might be one 
potential template solution: The lapack operations could be put into a class 
as static methods... I will think about this a bit more.)

> The reason I ask is that we would be interested in creating a
> multiprecision lapack library. We have already done this successfully
> with the GSL routines using a precompiler
> (http://www.win.ua.ac.be/~wschrep/precompiler ).

Interesting. You could use your precompiler to add yet another matrix class in 
lapack++, only with your MpIeee number type instead of Double/Int/Whatever. 
The double matrix is defined in the files include/{vd,gmd}.h and 
matrix/src/{vd,gmd}.cc , which would need to be converted by the precompiler 
-- and actually the implementation of various member methods is done by C 
macros and template functions, all in matrix/src/gmtmpl.cc. Simply for saving 
typing, and having only one function implementation for all number types.

> We wanted to tackle the lapack routines by modifying f2c (but this is
> tedious work since the f2c implementation is quite cryptic due to old
> style in coding). We hoped there might be an easier way using lapack++
> and some form of templates to get the job done.

Do you want to call fortran-lapack functions that have been compiled in 
Fortran for MpIeee? For those, you would still need the C header 
declarations. Those might indeed be converted by your precompiler, but you 
would need to check them by hand. If you have them, you could try to convert 
the rest of the functions in lapack++ to use those MpIeee lapack-functions. 
It shouldn't be too difficult, but I'm not sure if that is what you want.

> I tried TNT some time ago but it still contained many hardware double
> variables in the code giving bad results when trying to use our
> multiprecision number class MpIeee .

Might also be an issue in lapackpp. I don't know.

All the best,

Christian