Content

High-Level Interface for CXXLAPACK

You can easily create a high-level interface for CXXLAPACK that can be used like FLENS-LAPACK. This means that a user can use LAPACK without dealing with things like pointers, strides and leading dimensions. Only matrices and vectors get passed to the high-level interface. The low-level stuff happens inside of these functions.

We hope that at the end of this page you understand

Example Code

Instead of calling cxxlapack::getrf and cxxlapack::getrs directly we wrap them into functions that receive FLENS matrix/vector types. Inside of this wrapper we then extract information like stride, leading dimension and pointers and delegate functionality to the CXXLAPACK layer.

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

using namespace flens;

namespace mylapack {

template <typename MA, typename VPIV>
typename GeMatrix<MA>::IndexType
trf(GeMatrix<MA> &A, DenseVector<VPIV> &piv)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    return cxxlapack::getrf<IndexType>(A.numRows(),
                                       A.numCols(),
                                       A.data(),
                                       A.leadingDimension(),
                                       piv.data());
}

template <typename MA, typename VPIV, typename VB>
void
trs(Transpose trans, const GeMatrix<MA> &A, const DenseVector<VPIV> &piv,
    DenseVector<VB> &b)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

// namespace mylapack

A system of linear equations can now be solved easily. We use the functions mylapack::trf and mylapack:trs just like their FLENS-LAPACK counter parts:

#include <iostream>

#include "tut05-mylapack-version1.h"

using namespace flens;
using namespace std;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3.1;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    IndexType info = mylapack::trf(A, piv);

    if (info==0) {
        mylapack::trs(NoTrans, A, piv, b);

        cerr << "x = " << b << endl;
    } else {
        cerr << "A is numerically singular." << endl;
    }
}

Comments on Example Code

Let's have a closer look at the wrapper and how it gets used.

High-Level Interface (tut05-mylapack-version1.h)

Define USE_CXXLAPACK before you include the FLENS headers

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

Our high-level interface gets its own namespace

namespace mylapack {

We define our high-level interface for getrf it simply calls the CXXLAPACK interface getrf for the LAPACK functions dgetrf/zgetrf.

template <typename MA, typename VPIV>
typename GeMatrix<MA>::IndexType
trf(GeMatrix<MA> &A, DenseVector<VPIV> &piv)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    return cxxlapack::getrf<IndexType>(A.numRows(),
                                       A.numCols(),
                                       A.data(),
                                       A.leadingDimension(),
                                       piv.data());
}

Analogously we define a high-level interface for getrs. Note that here the right hand side b is a vector. So you might want to implement another variant of this function were the right hand side is a general matrix

template <typename MA, typename VPIV, typename VB>
void
trs(Transpose trans, const GeMatrix<MA> &A, const DenseVector<VPIV> &piv,
    DenseVector<VB> &b)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

Main Routine (tut05-page02-example1.cc)

Include our code snippet containing the high-level interface

#include "tut05-mylapack-version1.h"

We now can use our interface just the way we are using FLENS-LAPACK ...

    IndexType info = mylapack::trf(A, piv);

... for a user only the namspace differs.

        mylapack::trs(NoTrans, A, piv, b);

Compile and Run

In this example we link against vecLib which contains a LAPACK implementation. Note that the way we link against vecLib using the framework option is a Mac OS X specific feature.

$shell> cd flens/examples                                                      
$shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page02-example1 tut05-page02-example1.cc                                       
$shell> ./tut05-page02-example1                                                
A = 
            2             3            -1             0 
           -6            -5             0             2 
            2            -5             6            -6 
            4             6             2          -3.1 
b = 
           20            -33            -43             49 
x = 
     0.908163        9.18367        9.36735        9.18367 

Trouble (Part 1): Dealing with Views

We only apply minor modifications to the example code above:

#include <iostream>

#define USE_CXXLAPACK
#include <flens/flens.cxx>

#include "tut05-mylapack-version1.h"

using namespace flens;
using namespace std;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Underscore<IndexType>  _;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3.1;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    IndexType info = mylapack::trf(A, piv);

    if (info==0) {
        mylapack::trs(NoTrans, A, piv, b(_(1,n)));
//
//      This would work (but it might be inconvenient):
//
//      auto _b = b(_(1,n));
//      mylapack::trs(NoTrans, A, piv, _b);
//

        cerr << "x = " << b << endl;
    } else {
        cerr << "A is numerically singular." << endl;
    }
}

We use our high-level interface that we defined before.

#include "tut05-mylapack-version1.h"

We define an underscore object to ease the creation of matrix/vector views.

    Underscore<IndexType>  _;

Passing b(_(1,n)) as last argument to trs will cause a compile-time error. This is because the vector view b(_(1,n)) is a temporary object but trs expects a lvalue reference.

        mylapack::trs(NoTrans, A, piv, b(_(1,n)));
//
//      This would work (but it might be inconvenient):
//
//      auto _b = b(_(1,n));
//      mylapack::trs(NoTrans, A, piv, _b);
//

So let's have a look at the compiler errors:

$shell> cd flens/examples                                                      
$shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page02-example2 tut05-page02-example2.cc                                       
tut05-page02-example2.cc: In function 'int main()':
tut05-page02-example2.cc:57:49: error: invalid initialization of non-const reference of type 'flens::DenseVector, std::allocator > >&' from an rvalue of type 'flens::DenseVector >::View {aka flens::DenseVector, std::allocator > >}'
In file included from tut05-page02-example2.cc:9:0:
tut05-mylapack-version1.h:40:1: error: in passing argument 4 of 'void mylapack::trs(cxxblas::Transpose, const flens::GeMatrix&, const flens::DenseVector&, flens::DenseVector&) [with MA = flens::FullStorage; VPIV = flens::Array; VB = flens::ArrayView, std::allocator >]'
$shell> ./tut05-page02-example2                                                
/tmp/FLENS/shell.sh: line 3: ./tut05-page02-example2: No such file or directory

Actually the compiler even tells us how we can solve the problem!

Let's Fix It (Part 1)

The C++11 standard allows the definition of functions that receive non-const rvalue objects. So we change the definition of mylapack::trf and mylapack:trs:

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

using namespace flens;

namespace mylapack {

template <typename MA, typename VPIV>
typename GeMatrix<MA>::IndexType
trf(GeMatrix<MA> &A, DenseVector<VPIV> &piv)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    return cxxlapack::getrf<IndexType>(A.numRows(),
                                       A.numCols(),
                                       A.data(),
                                       A.leadingDimension(),
                                       piv.data());
}

template <typename MA, typename VPIV, typename VB>
void
trs(Transpose trans, const GeMatrix<MA> &A, const DenseVector<VPIV> &piv,
    DenseVector<VB> &&b)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

// namespace mylapack

Define USE_CXXLAPACK before you include the FLENS headers

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

Note the && in front of b. This means that b can now be non-const rvalue objects. Such objects arise when we create views using the underscore operator, e.g. by b(_(1,n)).

template <typename MA, typename VPIV, typename VB>
void
trs(Transpose trans, const GeMatrix<MA> &A, const DenseVector<VPIV> &piv,
    DenseVector<VB> &&b)
{
    typedef typename GeMatrix<MA>::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

We modify the second example such that it now includes our new versions of mylapack::trf and mylapack::trs:

#include <iostream>

#define USE_CXXLAPACK
#include <flens/flens.cxx>

#include "tut05-mylapack-version2.h"

using namespace flens;
using namespace std;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Underscore<IndexType>  _;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3.1;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    IndexType info = mylapack::trf(A, piv);

    if (info==0) {
        mylapack::trs(NoTrans, A, piv, b(_(1,n)));

        cerr << "x = " << b << endl;
    } else {
        cerr << "A is numerically singular." << endl;
    }
}

Now it compiles like a charm:

$shell> cd flens/examples                                                      
$shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page02-example3 tut05-page02-example3.cc                                       
$shell> ./tut05-page02-example3                                                
A = 
            2             3            -1             0 
           -6            -5             0             2 
            2            -5             6            -6 
            4             6             2          -3.1 
b = 
           20            -33            -43             49 
x = 
     0.908163        9.18367        9.36735        9.18367 

Trouble (Part 2): Let's go nuts!!!

Let's modify our very first example (where views did not occur at all) such that it also includes our new version of mylapack.

#include <iostream>

#define USE_CXXLAPACK
#include <flens/flens.cxx>

#include "tut05-mylapack-version2.h"

using namespace flens;
using namespace std;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3.1;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    IndexType info = mylapack::trf(A, piv);

    if (info==0) {
        mylapack::trs(NoTrans, A, piv, b);

        cerr << "x = " << b << endl;
    } else {
        cerr << "A is numerically singular." << endl;
    }
}

We now include the second version of mylapack.

#include "tut05-mylapack-version2.h"

Like in the first example we again pass a non-view object. So here b is a non-const lvalue again.

        mylapack::trs(NoTrans, A, piv, b);

For some fucking stupid reasons (long story) this will not compile:

$shell> cd flens/examples                                                      
$shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page02-example4 tut05-page02-example4.cc                                       
tut05-page02-example4.cc: In function 'int main()':
tut05-page02-example4.cc:51:41: error: cannot bind 'Vector {aka flens::DenseVector >}' lvalue to 'flens::DenseVector >&&'
In file included from tut05-page02-example4.cc:9:0:
tut05-mylapack-version2.h:32:1: error:   initializing argument 4 of 'void mylapack::trs(cxxblas::Transpose, const flens::GeMatrix&, const flens::DenseVector&, flens::DenseVector&&) [with MA = flens::FullStorage; VPIV = flens::Array; VB = flens::Array]'

Let's Fit It (Final Part)

The trick is exploiting a special template argument deduction rule for function templates. We don't want to go through the whole C++11 standard here. Just consider this declaration:

 template <typename T>                                                   
     void                                                                
     dummy(T &&t);                                                       

Note that here T is not part of a user defined type (i.e. we don't have something like GeMatrix<T> here). In this case we can pass either a lvalue or rvalue object to dummy:

We realize a few things:

However, we are not completely done yet. We don't want unrestricted template parameters. For this reason we provide some traits:

Analogously we provide traits for all other matrix and vector types:

This can be combined with the RestrictTo trait as we show in the final version for our mylapack functions

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

using namespace flens;

namespace mylapack {

template <typename MA, typename VPIV>
typename RestrictTo<IsGeMatrix<MA>::value
                 && IsIntegerDenseVector<VPIV>::value,
         typename RemoveRef<MA>::Type::IndexType>::Type
trf(MA &&A, VPIV &&piv)
{
    typedef typename RemoveRef<MA>::Type::IndexType  IndexType;
    return cxxlapack::getrf<IndexType>(A.numRows(),
                                       A.numCols(),
                                       A.data(),
                                       A.leadingDimension(),
                                       piv.data());
}

template <typename MA, typename VPIV, typename VB>
typename RestrictTo<IsGeMatrix<MA>::value
                 && IsIntegerDenseVector<VPIV>::value
                 && IsDenseVector<VB>::value,
         void>::Type
trs(Transpose trans, const MA &A, const VPIV &piv, VB &&b)
{
    typedef typename MA::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

// namespace mylapack

Define USE_CXXLAPACK before you include the FLENS headers

#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>

The trait is of form RestrictTo<Cond,B>::Type. If Cond is true then RestrictTo<Cond,B>::Type equals B. This will be the return type. If Cond is false then the trait can not be instantiated. In this case the function can not be instantiated either and the whole thing gets ignored.

Note that we also need the RemoveRef trait (yeah, C++ brain fuck!). If MA is of type T & then RemoveRef<T>::Type is of type T.

template <typename MA, typename VPIV>
typename RestrictTo<IsGeMatrix<MA>::value
                 && IsIntegerDenseVector<VPIV>::value,
         typename RemoveRef<MA>::Type::IndexType>::Type
trf(MA &&A, VPIV &&piv)
{
    typedef typename RemoveRef<MA>::Type::IndexType  IndexType;
    return cxxlapack::getrf<IndexType>(A.numRows(),
                                       A.numCols(),
                                       A.data(),
                                       A.leadingDimension(),
                                       piv.data());
}

For trs things are like above. But now the return type is simpler and for traits-newbies it might be easier to read.

template <typename MA, typename VPIV, typename VB>
typename RestrictTo<IsGeMatrix<MA>::value
                 && IsIntegerDenseVector<VPIV>::value
                 && IsDenseVector<VB>::value,
         void>::Type
trs(Transpose trans, const MA &A, const VPIV &piv, VB &&b)
{
    typedef typename MA::IndexType IndexType;
    cxxlapack::getrs<IndexType>(lapack::getF77Char(trans),
                                A.numRows(),
                                1,
                                A.data(),
                                A.leadingDimension(),
                                piv.data(),
                                b.data(),
                                b.length());
}

We can now conclude with an example using mylapack

#include <iostream>

#define USE_CXXLAPACK
#include <flens/flens.cxx>

// We now include the final version of `mylapack`.
#include "tut05-mylapack-version3.h"

using namespace flens;
using namespace std;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const Underscore<IndexType> _;

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3.1;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    IndexType info = mylapack::trf(A, piv);

    if (info==0) {
        const auto b_old = b;

//      We call trs with a lvalue ...
        mylapack::trs(NoTrans, A, piv, b);

        b = b_old;

//      ... and a rvalue. Just to see if now both cases work.
        mylapack::trs(NoTrans, A, piv, b(_(1,n)));

        cerr << "x = " << b << endl;
    } else {
        cerr << "A is numerically singular." << endl;
    }
}

... and finally ...

$shell> cd flens/examples                                                      
$shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page02-example5 tut05-page02-example5.cc                                       
$shell> ./tut05-page02-example5                                                
A = 
            2             3            -1             0 
           -6            -5             0             2 
            2            -5             6            -6 
            4             6             2          -3.1 
b = 
           20            -33            -43             49 
x = 
     0.908163        9.18367        9.36735        9.18367