Content |
CXXLAPACK
CXXLAPACK is a low-level interface for the original Fortran LAPACK. High-level interfaces for LAPACK can be built ontop of CXXLAPACK.
If functionality in FLENS-LAPACK is missing you can use CXXLAPACK as an alternative. In this case you need an external LAPACK implementation on your system. You either can receive it diretly from Netlib or through packages like VecLib, MKL, ACML, etc.
Example Code
#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
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;
b = 20,
-33,
-43,
49;
cerr << "A = " << A << endl;
cerr << "b = " << b << endl;
IndexType info = cxxlapack::getrf(A.numRows(),
A.numCols(),
A.data(),
A.leadingDimension(),
piv.data());
if (info==0) {
cxxlapack::getrs(lapack::getF77Char(NoTrans),
A.numRows(),
IndexType(1),
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.length());
/*
cxxlapack::getrs<IndexType>(lapack::getF77Char(NoTrans),
A.numRows(),
1,
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.stride());
*/
cerr << "x = " << b << endl;
} else {
cerr << "A is numerically singular." << endl;
}
}
Comments on Example Code
Define USE_CXXLAPACK before you include the FLENS headers
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>
The CXXLAPACK functions require the usual information describing the matrices, i.e. dimensions, leading dimension and a data pointer. All these can be retrieved from the FLENS matrix type.
A.numCols(),
A.data(),
A.leadingDimension(),
piv.data());
Note that cxxlapack::getrs has a template parameter for the used index type. That's why we call it with IndexType(1). This makes sure that values A.numRows(), IndexType(1), A.leadingDimension and b.stride() all have the same type.
A.numRows(),
IndexType(1),
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.length());
As an alternative we could specify the template parameter explictly:
cxxlapack::getrs<IndexType>(lapack::getF77Char(NoTrans),
A.numRows(),
1,
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.stride());
*/
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-page01-example1 tut05-page01-example1.cc $shell> ./tut05-page01-example1 A = 2 3 -1 0 -6 -5 0 2 2 -5 6 -6 4 6 2 -3 b = 20 -33 -43 49 x = 1 9 9 9
Another Example: Complex Numbers
We modify the previous example such that the matrix/vector element are complex valued. The functions in CXXLAPACK are overloaded for different element types. That means
-
cxxlapack::getrf is an interface for dgetrf and zgetrf,
-
cxxlapack::getrs is an interface for dgetrs and zgetrs.
Example Code
#include <iostream>
#ifndef USE_CXXLAPACK
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef complex<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 = 3;
Matrix A(n,n);
Vector b(n);
IndexVector piv(n);
A = T(1, 0), T( 1,-1), T( 2,20),
T(0, 1), T( 1, 2), T(-10, 8),
T(0,-1), T(-1, 1), T( 40, 6);
b = T( 1, 0),
T(-1, 1),
T(-2,-1);
cerr << "A = " << A << endl;
cerr << "b = " << b << endl;
IndexType info = cxxlapack::getrf(A.numRows(),
A.numCols(),
A.data(),
A.leadingDimension(),
piv.data());
if (info==0) {
cxxlapack::getrs(lapack::getF77Char(NoTrans),
A.numRows(),
IndexType(1),
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.length());
cerr << "x = " << b << endl;
} else {
cerr << "A is numerically singular." << endl;
}
}
Comments on Example Code
Define USE_CXXLAPACK before you include the FLENS headers
#define USE_CXXLAPACK
#endif
#include <flens/flens.cxx>
Elements are now defined to be complex
We use cxxlapack::getrf as we did in the previous example ...
A.numCols(),
A.data(),
A.leadingDimension(),
piv.data());
... and also cxxlapack::getrs.
A.numRows(),
IndexType(1),
A.data(),
A.leadingDimension(),
piv.data(),
b.data(),
b.length());
Compile and Run
$shell> cd flens/examples $shell> g++ -framework vecLib -std=c++11 -Wall -I../.. -o tut05-page01-example2 tut05-page01-example2.cc $shell> ./tut05-page01-example2 A = (1,0) (1,-1) (2,20) (0,1) (1,2) (-10,8) (0,-1) (-1,1) (40,6) b = (1,0) (-1,1) (-2,-1) x = (0,-1) (0,1) (0,0)