[Home]Linear Algebra With UBLAS

BOOST WIKI | RecentChanges | Preferences | Page List | Links List

Showing revision 24

Linear Algebra with uBLAS

Effective uBLAS

uBLAS is, more or less, BLAS (see http://www.netlib.org/blas/). It provides the "building blocks" for storing vector and matricies and vector and matrix operations (scaling, addition, multiplication etc). uBLAS provides the C++ infrastructure on which safe and efficient linear algebra algorithms can be built.

Solve and Factorisation in uBLAS

uBLAS has solve() functions [1], corresponding to the BLAS trsm() functions. These functions are `triangular' solvers. Triangular matrices have a trivial inverse, and therefore their solution is defined a BLAS function. The solve() function can by used for `backward' substitutions after LU or Cholesky factorisations.

uBLAS also has a set of functions for LU factorisation. These implement common LU factorisation algorithms directly in uBLAS. They are defined in the header file <boost/numeric/ublas/lu.hpp>. LU Matrix Inversion demonstrates how these functions can be used to invert matrices.

Interfacing with other Libraries

Linear Algebra is huge topic. There are many thousands of algorithms. Since the 70s large scale libraries of peer-reviewed algorithms have been assembled. These were usually written in Fortran. The most well know are LINPACK and LAPACK (see http://www.netlib.org ). More recently ATLAS has been developed. This has C and Fortran interfaces and includes both reimplemented BLAS functions as well as some important LAPACK algorithms.

Bindings

The best way to interface with these libraries is to use a `bindings' library which *can be used with uBLAS vectors and matrices*. Kresimir Fresl is continuously developing and updating such a library.

For example, for matrix inversion using an LU factorisation you can use ATLAS bindings:

 #include <boost/numeric/bindings/atlas/clapack.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 #include <boost/numeric/bindings/traits/std_vector.hpp>

 namespace ublas = boost::numeric::ublas;
 namespace atlas = boost::numeric::bindings::atlas;

 int main() {
      std::size_t n = 5;
      ublas::matrix<double> A (n, n);
      // fill matrix A

      std::vector<int> ipiv (n);   // pivot vector
      atlas::lu_factor (A, ipiv);  // alias for getrf()
      atlas::lu_invert (A, ipiv);  // alias for getri()
 }

Of course, you must install ATLAS first:

         http://math-atlas.sourceforge.net/

"The ATLAS (Automatically Tuned Linear Algebra Software) project is an ongoing research effort focusing on applying empirical techniques in order to provide portable performance. At present, it provides C and Fortran77 interfaces to a portably efficient BLAS implementation, as well as a few routines from LAPACK."

The `bindings' library is a very useful interface between Boost and uBLAS and traditional linear algebra libraries. It is being actively maintianed with the hope of eventual inclusion in Boost.

Development of the binding interface is located at Boost Sandbox CVS - http://www.boost.org/more/mailing_lists.htm#sandbox

Here you can also find newest, sometimes experimental,versions of the binding. You can access the header files directly at:

http://boost-sandbox.cvs.sourceforge.net/boost-sandbox/boost-sandbox/boost/numeric/bindings/

The starting point for documentation and the many examples is:

       boost-sandbox/boost-sandbox/libs/numeric/bindings

Examples for matrix inversion with ATLAS are in `atlas/ublas_getri.cc' and,if your matrix is SPD, in `atlas/ublas_potri.cc'.

You may also download ready-made snapshot releases of the bindings from http://news.tiker.net/software/boost-bindings.

What about using LAPACK++?

From the LAPACK++ home page:

 NOTE: This package is being superseded by the Template Numerical
 Toolkit (TNT), which utilizes new features of the ANSI C++ specification.
 TNT is a newer design, and will integrate the functionlaity of Lapack++,
 IML++, SparseLib?++, and MV++.
 (TNT home page is: http://math.nist.gov/tnt/)

On the other hand, from the homepage of LAPACK++ v2.1.2: http://lapackpp.sourceforge.net/

 However, they abandoned LAPACK in the year 2000 and stated: "Lapack++ is no longer actively supported.
 The successor to this project is that Template Numerical Toolkit (TNT), see http://math.nist.gov/tnt for details."
 Unfortunately, the project TNT never really took off.
 Therefore this fork from the original LAPACK++ has been started.
 There are a whole number of changes now in here.
 Most notably, this local copy has complex matrices enabled again by adding a custom copy of stdc++'s  complex type (see include/lacomplex.h and include/lacomplex).

So my question, again: What about LAPACK++ ?

Blaisorblade answers:

(I'm not a Boost developer) LAPACK++ seems fairly complete, but it seems not written by C++ experts. Beyond that, it contains its own reimplementation of basic vector operation (indeed, its own reimplementation of BLAS for C++), but it is a really inferior one (it does not definitely use expression templates).

Say, in Lapack++ 2.5.2, LaGenMatDouble? class (a double dense matrix) has two identical methods, having these signatures:

      LaGenMatDouble? operator()(const LaIndex?& I, const LaIndex?& J) ;
      LaGenMatDouble? operator()(const LaIndex?& I, const LaIndex?& J) const;
the non-const overload returns a copy of the selected portion of the matrix, not a view. I wonder if they understand that the second overload is not needed.

In the overload of operator() to access a single element (operator()(int row, int column)), a reference to the element is returned; I wonder whether they'd need to return a proxy reference, since they can share the backing array.

However, I downloaded LAPACK++ yesterday, and gave just a look. I'll compare it with the bindings project, to compare their respective merits.

Back to Effective uBLAS


BOOST WIKI | RecentChanges | Preferences | Page List | Links List
Edit revision 24 of this page | View other revisions | View current revision
Edited December 10, 2007 11:57 am (diff)
Search:
Disclaimer: This site not officially maintained by Boost Developers