Content |
Using SuperLU with CCS
This example requires that SuperLU (Version 4.3) is installed on your system.
In this example we will:
-
Setup a general sparse matrix in coordinate storage.
-
Show how to convert the sparse matrix from the coordinate storage format to the compressed column storage format. This format is required by SuperLU.
-
Show how to call the solver from SuperLU. We use quick and dirty hacked wrapper function for this purpose. You certainly can do this more elegant.
Example Code
The following example is basically taken from the SuperLU user manual and ported to FLENS.
#include <flens/flens.cxx>
using namespace flens;
using namespace std;
#include <slu_ddefs.h>
template <typename MA, typename PR, typename PC, typename VB>
int
dgssv(GeCCSMatrix<MA> &A,
DenseVector<PC> &pc,
DenseVector<PR> &pr,
DenseVector<VB> &b)
{
ASSERT(pr.length()==A.numRows());
ASSERT(pc.length()==A.numCols());
superlu_options_t options;
SuperLUStat_t stat;
SuperMatrix A_, L_, U_, B_;
dCreate_CompCol_Matrix(&A_,
A.numRows(), A.numCols(), A.engine().numNonZeros(),
A.engine().values().data(),
A.engine().rows().data(),
A.engine().cols().data(),
SLU_NC, SLU_D, SLU_GE);
dCreate_Dense_Matrix(&B_,
b.length(), 1, b.data(), b.length(),
SLU_DN, SLU_D, SLU_GE);
set_default_options(&options);
options.ColPerm = NATURAL;
StatInit(&stat);
int info;
dgssv(&options, &A_, pc.data(), pr.data(), &L_, &U_, &B_, &stat, &info);
Destroy_SuperMatrix_Store(&A_);
Destroy_SuperMatrix_Store(&B_);
Destroy_SuperNode_Matrix(&L_);
Destroy_CompCol_Matrix(&U_);
StatFree(&stat);
return info;
}
int
main()
{
typedef int IndexType;
typedef IndexBaseZero<IndexType> IndexBase;
typedef CoordStorage<double, CoordColRowCmp, IndexBase> Coord;
const IndexType m = 5;
const IndexType n = 5;
GeCoordMatrix<Coord> A_(m, n);
const double s = 19,
u = 21,
p = 16,
e = 5,
r = 18,
l = 12;
A_(0,0) += s;
A_(1,1) += u;
A_(2,2) += p;
A_(3,3) += e;
A_(4,4) += r;
A_(1,0) += l; A_(2,1) += l; A_(4,0) += l; A_(4,1) += l;
A_(0,2) += u; A_(0,3) += u; A_(3,4) += u;
GeCCSMatrix<CCS<double, IndexBase> > A = A_;
std::cout << "A_ = " << A_ << endl;
std::cout << "A = " << A << endl;
DenseVector<Array<double, IndexBase> > b(m);
b = 1;
DenseVector<Array<IndexType, IndexBase> > pr(m), pc(n);
int info = dgssv(A, pr, pc, b);
if (info==0) {
cout << "x = " << b << endl;
} else {
cout << "SuperLU dgssv: info = " << info << endl;
}
return info;
}
Some Notes
Include the SuperLu stuff for double precision.
Wrapper for the SuperLU solver dgssv. Looks complicated but it's actually pretty simple. We create SuperLU matrix wrapper for our matrices and then call the solver. Anyway, in future we will hide the details in FLENS.
int
dgssv(GeCCSMatrix<MA> &A,
DenseVector<PC> &pc,
DenseVector<PR> &pr,
DenseVector<VB> &b)
{
ASSERT(pr.length()==A.numRows());
ASSERT(pc.length()==A.numCols());
superlu_options_t options;
SuperLUStat_t stat;
SuperMatrix A_, L_, U_, B_;
dCreate_CompCol_Matrix(&A_,
A.numRows(), A.numCols(), A.engine().numNonZeros(),
A.engine().values().data(),
A.engine().rows().data(),
A.engine().cols().data(),
SLU_NC, SLU_D, SLU_GE);
dCreate_Dense_Matrix(&B_,
b.length(), 1, b.data(), b.length(),
SLU_DN, SLU_D, SLU_GE);
set_default_options(&options);
options.ColPerm = NATURAL;
StatInit(&stat);
int info;
dgssv(&options, &A_, pc.data(), pr.data(), &L_, &U_, &B_, &stat, &info);
Destroy_SuperMatrix_Store(&A_);
Destroy_SuperMatrix_Store(&B_);
Destroy_SuperNode_Matrix(&L_);
Destroy_CompCol_Matrix(&U_);
StatFree(&stat);
return info;
}
SuperLU requires an index base of zero. Here we set the default index base of the storage scheme to zero via IndexBaseZero.
typedef IndexBaseZero<IndexType> IndexBase;
typedef CoordStorage<double, CoordColRowCmp, IndexBase> Coord;
Alternative we could specify it through the constructor (see class API).
const IndexType n = 5;
GeCoordMatrix<Coord> A_(m, n);
We setup the matrix from the SuperLU user guide. So here the values.
u = 21,
p = 16,
e = 5,
r = 18,
l = 12;
Matrices in coordinate storage are usually used in FEM for assembling the stiffness matrix. Values will be accumulated. Note that an assignment like 'A(2,3) = ' is not allowed.
A_(1,1) += u;
A_(2,2) += p;
A_(3,3) += e;
A_(4,4) += r;
A_(1,0) += l; A_(2,1) += l; A_(4,0) += l; A_(4,1) += l;
A_(0,2) += u; A_(0,3) += u; A_(3,4) += u;
Convert to compressed column storage.
Just for curiosity: Compare the two formats.
std::cout << "A = " << A << endl;
Setup the right-hand side \(b\). We set all values to \(1\).
b = 1;
Call the SuperLU solver for solving \(Ax = b\). Note that on exit \(b\) is overwritten with the solution \(x\).
int info = dgssv(A, pr, pc, b);
Check the info code
cout << "x = " << b << endl;
} else {
cout << "SuperLU dgssv: info = " << info << endl;
}
Compile
$shell> cd flens/examples $shell> clang++ -Wall -std=c++11 -I../.. -I$HOME/work/SuperLU_4.3/SRC -L$HOME/work/SuperLU_4.3/lib -lsuperlu_4.3 -framework vecLib -o geccs-superlu geccs-superlu.cc
In this case we assume that SuperLU is installed in $HOME/SuperLU_4.3/. So the following options were used:
|
-I$HOME/work/SuperLU_4.3/SRC |
Add the SuperLU headers to the include path. |
|
-L$HOME/work/SuperLU_4.3/lib |
Add the SuperLU library path. |
|
-lsuperlu_4.3 |
Link against the SuperLU library. |
|
-framework vecLib |
This is the Mac OS X way of linking BLAS. On other platforms this might be done more manual, e.g -latlas -lblas. |
Run
$shell> cd flens/examples $shell> ./geccs-superlu A_ = isSorted_ = 1 #0: (0, 0) = 19 #1: (1, 0) = 12 #2: (4, 0) = 12 #3: (1, 1) = 21 #4: (2, 1) = 12 #5: (4, 1) = 12 #6: (0, 2) = 21 #7: (2, 2) = 16 #8: (0, 3) = 21 #9: (3, 3) = 5 #10: (3, 4) = 21 #11: (4, 4) = 18 A = compressed cols: cols: 0 3 6 8 10 12 rows: 0 1 4 1 2 4 0 2 0 3 3 4 values: 19 12 12 21 12 12 21 16 21 5 21 18 x = -0.03125 0.0654762 0.0133929 0.0625 0.0327381