# Linear Algebra With UBLAS

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

Difference (from prior minor revision) (major diff, author diff)

Changed: 91c91,92
 LaGenMatDouble? operator()(const LaIndex?& I, const LaIndex?& J) ;
 LaGenMatDouble? operator()(const LaIndex?& I, const LaIndex?& J);

Changed: 97,98c99,100
 const LaGenMatDouble? A = LaGenMatDouble?::rand(3, 3); A(0, 0) = 1;
 const LaGenMatDouble? A = LaGenMatDouble?::rand(3, 3); A(0, 0) = 1;

## 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:

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'.

``` 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++.
```

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++ ?

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

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