Re: [Boost-Users] ublas: questions from a newcomer

From: Kresimir Fresl (fresl_at_[hidden])
Date: 2003-03-25 11:24:30

stuart macgregor wrote:

 > I am interested in ublas for linear algebra.  In addition to the simple
 > 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,

Atlas bindings work for me in the `real' code (but I am biased ;o)
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
 > set of lapack routines?

Some time ago Toon wrote in boost-ml:

``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
 > spec would be nice (ruby/perl?

Interesting idea.

 > or just a description of what needs to be
 > done)  Then I could happily bring up the interfaces I need.

As I said, bindings are free-time activity, so any help is
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
 > from a potential user.

Well, some boost libraries are certainly frightening, but I don't think
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