Linear Algebra With UBLAS

Linear Algebra with 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 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

Bindings

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(not a Boost developer) answers (there are updates from the first answer):

LAPACK++ seems fairly complete (and its API seems a lot more reasonable than the API of bindings to LINPACK), 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. EDIT: things are worse; _both_ operators return a _modifiable_ view, and what's worse is that constness is never enforced. The following code snippet compiles:

      const LaGenMatDouble? A = LaGenMatDouble?::rand(3, 3);
      A(0, 0) = 1;

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.

EDIT: additionally, the LAPACK++ developers do not know about thread safety, and they ignore the very existence of mutable. Look at this comment from source code:

      static int      *info_;   // print matrix info only, not values
      //   originally 0, set to 1, and then
      //   reset to 0 after use.
      // use as in
      //
      //    std::cout << B.info() << std::endl;
      //
      // this *info_ member is unique in that it really isn't
      // part of the matrix info, just a flag as to how
      // to print it.   We've included in this beta release
      // as part of our testing, but we do not expect it
      // to be user accessable.
      // It has to be declared as global static
      // so that we may monitor expresssions like
      // X::(const &X) and still utilize without violating
      // the "const" condition.
      // Because this *info_ is used at most one at a time,
      // there is no harm in keeping only one copy of it,
      // also, we do not need to malloc free space every time
      // we call a matrix constructor.

Back to Effective uBLAS