Content |
Solving Systems of Linear Equations
In this example we solve a system of linear equations \(Ax = b\) were the coefficient matrix is symmetric positiv definite and banded. We solve the system with lapack::pbsv which is the FLENS interface for LAPACK's dpbsv.
Note that we might rename lapack::pbsv to lapack::posv.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
typedef SbMatrix<BandStorage<T> > SymmetricBandMatrix;
typedef DenseVector<Array<T> > Vector;
typedef SymmetricBandMatrix::IndexType IndexType;
const IndexType n = 4;
SymmetricBandMatrix A(n, Upper, 1);
Vector b(n);
A.general().diag( 1) = -1;
A.general().diag( 0) = 2;
b = 1,
2,
3,
4;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::pbsv(A, b);
cout << "(L\\R) = " << A << endl;
cout << "x = " << b << endl;
}
Comments on Example Code
Define some convenient typedefs for the matrix/vector types of our system of linear equations.
typedef DenseVector<Array<T> > Vector;
We also need an extra vector type for the pivots. The type of the pivots is taken for the system matrix.
Set up the baby problem ...
We allocate a symmetric band matrix with one off-diagonal.
We initialize the positiv definite \(n \times n\) tridiagonal matrix.
A.general().diag( 0) = 2;
And solve it with lapack::pbsv
Compile
Note that we need to link against an external LAPACK implementation:
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -DUSE_CXXLAPACK -framework vecLib -o lapack-pbsv lapack-pbsv.cc
Run
$shell> cd flens/examples $shell> ./lapack-pbsv A = 2.000000 -1.000000 -1.000000 2.000000 -1.000000 -1.000000 2.000000 -1.000000 -1.000000 2.000000 b = 1.000000 2.000000 3.000000 4.000000 (L\R) = 1.414214 -0.707107 -0.707107 1.224745 -0.816497 -0.816497 1.154701 -0.866025 -0.866025 1.118034 x = 4.000000 7.000000 8.000000 6.000000