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 first compute the Cholesky factorization \(S = L L^T\) with lapack::pbtrf which is the FLENS port of LAPACK's dpbtrf. We then solve \(Lu=b\) and then \(L^T x = u\) using lapack::pbtrs which is the port of dpbtrs.
Note that we might rename lapack::pbtrf to lapack::potrf and lapack::pbtrs to lapack::potrs.
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::pbtrf(A);
cout << "L^T = " << A.general().upper() << endl;
lapack::pbtrs(A, b);
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;
We factorize \(A\) with lapack::pbtrf
In the upper part of \(A\) now the triangular factor \(L^T\) is stored.
Solve it with lapack::pbtrs
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-pbtrs lapack-pbtrs.cc
Run
$shell> cd flens/examples $shell> ./lapack-pbtrs 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^T = 1.414214 -0.707107 1.224745 -0.816497 1.154701 -0.866025 1.118034 x = 4.000000 7.000000 8.000000 6.000000