- stuart macgregor wrote:

> I am interested in ublas for linear algebra. In addition to the simple

Atlas bindings work for me in the `real' code (but I am biased ;o)

> operations supported by ublas I will need comprehensive access to lapack

> for the real work.

> I see in the wiki that some work has been done in this area. Is this

> numeric/binding stuff working,

I am also using SuperLU bindings for sparse matrices, but they are

not in boost-sandbox yet.

> and can it be readily extended to the full

Some time ago Toon wrote in boost-ml:

> set of lapack routines?

``Note that the bindings are not complete but the 'framework' is there

and it's just a matter of time to complete them.''

Bindings are mostly a free-time activity. Or, one can say that we

try to generalize some more or less ad-hoc parts of our `regular' work.

So `matter of time' probably means `some (considerable) time' :o(

> Some kind of meta program to generate lapack interfaces from a parameter

Interesting idea.

> spec would be nice (ruby/perl?

> or just a description of what needs to be

As I said, bindings are free-time activity, so any help is

> done) Then I could happily bring up the interfaces I need.

welcomed.

Regarding (c)lapack part of atlas bindings, there are interfaces

for all (interesting) functions which atlas provides (gesv, getrf,

getrs, getri, posv, potrf, potrs, potri; only lauum and trtri are

missing), see

http://math-atlas.sourceforge.net/psdoc/lapackqref.ps.gz

Lapack bindings are Toon's realm, but I hope that my understanding

is good enough to (try to) explain what needs to be done.

As an example, I'll use function `getrf' (which is already there).

As you know, there are four versions of each lapack function

-- float, double, complex, complex double. AFAICS Toon decided

to write bindings only for double and complex double versions,

but other two can be easily added following the same pattern.

I'll use only dgetrf here.

Include files are in

boost/numeric/bindings/lapack

and the `traits.cpp' file is in

libs/numeric/bindings/lapack

(1) `lapack.h' contains declarations of lapack functions:

As there are several conventions for calling Fortran from C/C++

(some compilers require trailing underscore in function names,

some don't, see excelent page

http://www.aei.mpg.de/~jthorn/c2f.html

and, in particular, section `Name mangling'), first a macro

LAPACK_DGETRF is defined, depending on naming convention:

=============================================

#if defined (BIND_FORTRAN_LOWERCASE_UNDERSCORE)

#define LAPACK_DGETRF dgetrf_

#elif defined(BIND_FORTRAN_LOWERCASE)

#define LAPACK_DGETRF dgetrf

#endif

=============================================

(macros BIND_FORTRAN_... are defined in

boost/numeric/bindings/traits/fortran.h).

Now lapack function can be declared:

=============================================

void LAPACK_DGETRF(const int* m, const int* n,

double* a, const int* lda, int* ipiv, int*

info) ;

=============================================

All parameters are pointers.

(2) `traits.hpp' contains class template `traits'.

For each lapack there is a typedef of a function pointer, which

corresponds to the function declaration in `lapack.h':

=============================================

typedef void (*getrf_type)(const int* m, const int* n,

bind_type* a, const int* lda, int* ipiv, int* info) ;

=============================================

but instead of `double*' there is a `bind_type*', which depends

on the template parameter `T' of the class `traits'. `bind_type' is

defined using `value_traits' class (which is defined and specialized

in boost/numeric/bindings/traits/value_traits.hpp).

Class `traits' also contains a (static) function pointer:

=============================================

static getrf_type getrf ;

=============================================

(3) in `traits.cpp' the function pointer is initialized with

(the address of) the lapack function:

=============================================

template<> traits<double>::getrf_type

traits<double>::getrf = LAPACK_DGETRF ;

=============================================

`traits.cpp' must be separately compiled and linked with your

program.

But, in fact, here we have a problem -- see comment [4] below.

(4) `lapack.hpp' contains the C++ interfaces (or bindings).

All bindings (blas, lapack, atlas) are written in terms of traits

classes which are interfaces between these functions and

vector/matrix classes. Namely, idea is to make bindings

independent of the concrete vector/matrix libraries as much

as possible.

Lapack functions require the starting address of the matrix,

number of rows and columns and leading dimension.

So, matrix_traits class, that is, its specializations provide this

functions and some typedefs. For example, size1() (i.e. number

of rows) in matrix_traits specialization for ublas::matrix traits is

static int size1 (matrix_type& m) { return m.size1(); }

and in specialization for Array2D from Roldan Pozo's TNT

static int size1 (matrix_type& m) { return m.dim1(); }

If you use ublas, you don't need to worry about matrix_traits,

because they are already written -- you just have to include

boost/numeric/bindings/traits/ublas_matrix.hpp

for general matrices and matrix proxies and

boost/numeric/bindings/traits/ublas_symmetric.hpp

for symmetric matrices and symmetric adaptors.

(Traits classes for Roldan Pozo's TNT are also

there, but I didn't test them thoroughly.)

Note that bindings don't know and need not know anything

about traits specializations --

boost/numeric/bindings/traits/ublas_matrix.hpp etc.

must be included in your program, not in bindings

headers.

To simplify the use of traits classes, Toon recently introduced

some free accessor function, so instead of writting

matrix_traits<matrix_type>::size1 (m)

you can write simply

matrix_size1 (m)

Now, let's see how all this works in the interface for getrf

(explanations are given as comments; this is Toon's code and I'll

also comment some implementation details -- longer comments,

not really important for functionality, are given after the code):

==============================================

template < typename matrix_type, typename IpivIterator >

int getrf (matrix_type& a, IpivIterator begin_ipiv, IpivIterator end_ipiv)

// instead of 6 parameters of lapack's (sdcz)getrf, we have only

// 3 parameters: matrix to factor and iterator range for pivot vector

// (but see also comment [1] below)

{

using namespace boost::numeric::bindings::traits ;

// why? all names below are fully namespace qualified

typedef typename matrix_type::value_type value_type ;

// see comment [2] below

typedef typename

boost::numeric::bindings::traits::value_traits< value_type>::value_type

bind_type ;

// as in `traits.hpp' in (2) above

// int m = boost::numeric::bindings::traits::size1( a ) ;

// int n = boost::numeric::bindings::traits::size2( a ) ;

// number of matrix rows and columns;

// free accessor functions are used, but in new version they are

// matrix_size1() and matrix_size2(), so this will not compile :o(

// correct calls are:

int m = boost::numeric::bindings::traits::matrix_size1( a ) ;

int n = boost::numeric::bindings::traits::matrix_size2( a ) ;

value_type *a_ptr =

boost::numeric::bindings::traits::matrix_storage( a ) ;

const int lda =

boost::numeric::bindings::traits::leading_dimension( a ) ;

// starting address of matrix storage and its leading dimension

assert( std::distance( begin_ipiv, end_ipiv ) >= std::min( m, n ) );

// pivot vector must be large enough

int info = 0; // return value; in short: is matrix singular?

int& ipiv_ptr = *begin_ipiv ;

// reference to first element of the pivot vector,

// &ipiv_ptr is then its starting address;

// but see also comment [3]

traits<value_type>::getrf( &m, &n, (bind_type*)a_ptr, &lda,

&ipiv_ptr, &info );

// finally, call lapack's function

return info;

}

==============================================

Comments:

[1] I think that this interface [i.e. func (matrix&, begin_iter,

end_iter)] is

mixing of two idioms -- why STL's idiom of iterator ranges [begin, end)

here, and not for other vectors in e.g. blas 1&2 bindings?

Interface of getrf in atlas bindings (atlas/clapack.hpp) is simply:

int getrf (MatrA& a, IVec& ipiv);

[2] It is true that both ublas::matrix<> and TNT::Array2D<> define

value_type, but (in ideal world ;o) bindings should work for other

matrix types which can have val_t or something else.

And that's why we have matrix_traits; so, instead of

typedef typename matrix_type::value_type value_type ;

it should be

typedef typename matrix_traits<matrix_type>::value_type value_type;

[3] See the interface of getrf in atlas bindings (atlas/clapack.hpp)

for IMHO more elegant way to do this (yes, I'm probably biased ;o)

[4] The biggest problem with blas & lapack bindings is that

although they work with g++ (at least with 3.2, I didn't try later

versions), they cannot be compiled with Comeau C++ 4.3

--- they are not standard conforming.

At first sight, class `traits<>', defined in `lapack/traits.hpp'

(and, similarly, class `blas<>' in `blas/blaspp.hpp') is very elegant

solution for the function overloading -- one class for all four

versions of blas/lapack functions (float, double, complex,

complex double). But function pointers defined in `traits<>'

and `blas<>' have "C++" linkage and blas/lapack functions

have (extern) "C" linkage; therefore, this may, but need not,

work:

template<> traits<double>::getrf_type

traits<double>::getrf = LAPACK_DGETRF ;

Here is a quote from the recent post by Greg Comeau in

comp.lang.c++.moderated (`Re: use of extern "C"'):

===============================================

One reason why this has probably been acceptable in the code

you've written is that your compilers are allowing implicit

conversions between extern "C" and extern "C++" entities,

actually pointers to some such entities, since often the

same internal/underlying conventions (calling, representations, etc)

are used for each, but from a strict standpoint, the standard

does not support that.

===============================================

Steve Clamage gave an example in the same thread:

===============================================

The extern "C" linkage declaration says that C

linkage rules are in effect. If C and C++ functions have different

calling sequences, the code generation will be different depending on

the declared linkage. The same goes for pointer-to function, since

the linkage is part of the pointer type. For example:

int f(int); // C++ linkage

extern "C" int g( int(*)(int) );

...

g(f); // error

Function g is declared to take a pointer to a C-linkage function, and

passing a pointer to a C++-linkage function is not valid. Not all

compilers complain about the error, however.

===============================================

Our example is inverted -- traits<double>::getrf_type is a pointer

to a C++-linkage function and LAPACK_DGETRF is (a pointer)

to a C-linkage function -- but the same rule obviously applies.

I am afraid that blas and lapack bindings need complete rewriting

as in atlas bindings which use less elegant and more conservative

(but also standard conforming) approach -- simple function overloading:

===============================================

int gesv (CBLAS_ORDER const Order, int const N, int const NRHS,

float* A, int const lda, int* ipiv, float* B, int const ldb)

{

return clapack_sgesv (Order, N, NRHS, A, lda, ipiv, B, ldb);

}

int gesv (CBLAS_ORDER const Order, int const N, int const NRHS,

double* A, int const lda, int* ipiv, double* B, int const

ldb)

{

return clapack_dgesv (Order, N, NRHS, A, lda, ipiv, B, ldb);

}

===============================================

> Is there any way to generate vectore and matrices which do not own

the data

> buffer which they provide access to - 'view' or 'sub' objects in

other lin

> alg packages? This would be very valuable to me for interfacing

existing

> code, without duplicating storage and doing unnecessary copies. I

found it

> difficult to determine from the current documentation.

Joerg answered this one ...

> I have one worry: Ublas, and most of the boost modules, seems to me

not for

> the faint hearted. They are clearly powerful, but seem to expect a lot

Well, some boost libraries are certainly frightening, but I don't think

> from a potential user.

that ublas is among them -- its interface seems rather intuitive.

> Are these facilities aimed at library writers rather than users whoes

skill

> focus may lie in other domains?

Some very simple, self-explanatory ublas examples can be found in

subdirectory

libs/numeric/ublas/doc/samples

of the boost distribution.

> I expect that a user level doc may eventually make ublas/lapack much

more

> approachable by a non computer science mathematician or engineer.

I agree. ublas documentation is good and detailed reference

(when you know what you are looking for ;o), but a tutorial with

few well commented examples is probably missing.

Regarding bindings, I can only hope that I'll find some time in the near

future to write some documentation (at least for traits classes and atlas

bindings). In the meantime, there are some examples in

libs/numeric/bindings/atlas/

(in the boost-sandbox); I tried to write meaningful comments ;o).

Begin with ublas_gesv2.cc, then ublas_gesv.cc and

ublas_posv.cc.

Also, don't hesitate to ask further questions.

(BTW, few days ago Toon updated bindings in the sandbox,

so it now contains the latest version -- that is, the note on

Wiki page is obsolete.)

> Forgive me if the above sounds critical and ungracious for free

software,

No, not at all ... I had (and have) the same problem with some other

libraries ;o)

> but I want to use boost and I am worried that it may be too demanding

for

> its benefites to be justified for use by the group of people I work

with.

AFAICS, the main problem is the lack of documentation, in particular

tutorials and introductory examples. Also, as Joerg suggested,

``compiler diagnostics clearly are a pain sometimes'', but this is

a problem with all templated code.

Regards,

fres - Hi Matthias (hi all),

you wrote:

> After some discussions in my group when we had Jeremy Siek visit this

I've been thinking about the lapack bindings for some time now and I believe

> weekend we agreed that we would like to see generic C++ lapack

> interfaces that work with ublas and other matrix libraries submitted to

> boost as a library. One could start with the high level driver

> functions only and add more specialized functions later. I thus would

> like to ask what the current status of the ublas lapack bindings is and

> if somebody is working on a generic version for boost, or if my

> research group could help here?

it's firstly necessary to limit the scope of investigation.

The LAPACK users guide distinguishes the following parts:

<cite>

(1) Linear Equations

(2) Linear Least Squares (LLS) Problems

(3) Generalized Linear Least Squares (LSE and GLM) Problems

(4) Standard Eigenvalue and Singular Value Problems

(a) Symmetric Eigenproblems (SEP)

(b) Nonsymmetric Eigenproblems (NEP)

(c) Singular Value Decomposition (SVD)

(5) Generalized Eigenvalue and Singular Value Problems

(a) Generalized Symmetric Definite Eigenproblems (GSEP)

(b) Generalized Nonsymmetric Eigenproblems (GNEP)

(c) Generalized Singular Value Decomposition (GSVD)

</cite>

This structure seems to suggest some splitting into sublibraries or

namespaces. Which parts are most important?

Golub & van Loan reference the following procedures:

Chapter 3:

(1) _TRSV, _TRSM, _TRCON, _TRRFS, _TRTRS, _TRTRI

(1) _GESV, _GECON, _GERFS, _GESVX, _GETRF, _GETRS, _GETRI, _GEEQU

Chapter 4:

(1) _GBSV, _CGBCON, _GBRFS, _GBSVX, _GBTRF, _GBTRS, _GBEQU

(1) _GTSV, _GTCON, _GTRFS, _GTSVX, _GTTRF, _GBTRS

(1) _POSV, _POCON, _PORFS, _POSVX, _POTRF, _POTRS, _POTRI, _POEQU

(1) _PBSV, _PBCON, _PBRFS, _PBSVX, _PBTRF, _PBTRS

(1) _PTSV, _PTCON, _PTRFS, _PTSVX, _PTTRF, _PTTRS

(1) _SYSV, _SYCON, _SYRFS, _SYSVX, _SYTRF, _SYTRS, _SYTRI

(1) _TBCON, _TBRFS, _TBTRS

Chapter 5:

(?) _LARFG, _LARF, _LARFX, _LARFB, _LARFT

(?) _LARTG, _LARGV, _LARTV, _LARSR, CSROT, CROT, CLACGV

(?) _GEQRF, _GEQPF, _ORMQR, _UNMQR, _ORGQR, _UNGQR

(?) _GERQF, _GELQF, _GEQLF, _TZRQF

(4.c) _GESVD, _BDSQR, _GEBRD, _ORGBR, _GBBRD

(2) _GELS, _GELSS, _GELSX, _GEEQU

Chapter 7:

(4 b) _GEBAL, _GEBAK

(4 b) _GEHRD, _ORMHR, _ORGHR, _UNMHR, _UNGHR

(4 b) _HSEQR, _HSEIN

(4 b) _GEES, _GEESX, _GEEV, _GEEVX

(4 b) _TREVC, _TRSNA, _TREXC, _TRSEN, _TRSYL

(5 b) _GGBAL, _GGHRD, _HGEQZ, _TGEVC, _GGBAK

Chapter 8:

(4 a) _SYEV, _SYEVD, _SYEVX

(4 a) _SYTRD, _SBTRD, _SPTRD

(4 a) _STEQR, _STEDC, _STERF, _PTEQR, _STEBZ, _STEIN

(?) _SYGST, _PBSTF, _SBGST

(4 c) _GESVD, _BDSQR, _GEBRD, _ORGBR, _GBBRD

(5) _GGSVP, _TGSJA

My classification currently is an educated guess only ;-) If I count these

function correctly, I still find an astonishing number of around 130 ;-( To

get any clue of how to proceed I believe we have no other chance than to

follow Kresimir's suggestion to invent names for the corresponding families

of overloaded functions or specialized classes. If nobody else is

interested, I'll try to come up with a proposal after some more scanning...

All the best,

Joerg

P.S.: I didn't see any definitive motion for your idea to wrap up the driver

routines first, so this is kind of an alternative ansatz.

P.P.S.: How are we able to test these wrappers?