Content |
Solving a Symmetric Positive Definite System of Linear Equations
In this example we solve a system of linear equations \(Sx = b\) were the coefficient matrix is symmetric and positiv definite. We first compute the Cholesky factorization \(S = L L^T\) with lapack::potrf which is the FLENS port of LAPACK's dpotrf. We then solve \(Lu=b\) and then \(L^T x = u\) using lapack::potrs which is the port of dpotrs.
Example Code
#include <iostream>
using namespace flens;
using namespace std;
int
main()
{
GeMatrix<FullStorage<double> > A(3, 3);
DenseVector<Array<double> > b(3);
A = 2, 1, 0,
0, 2, 1,
0, 0, 2;
b = 1,
2,
3;
auto S = A.upper().symmetric();
cout << "S = " << S << endl;
cout << "b = " << b << endl;
int info = lapack::potrf(S);
if (info!=0) {
cerr << "S is singular" << endl;
return info;
}
lapack::potrs(S, b);
cout << "inv(S)*b = " << b << endl;
return 0;
}
Comments on Example Code
Setup the raw data.
0, 2, 1,
0, 0, 2;
\(S\) is a symmetric matrix view constructed from the upper triangular part of \(A\). Note that \(S\) only references data from \(A\).
Computes the Cholesky factorization \(S = L L^T\) where \(L\) is upper triangular. Matrix \(S\) (i.e. the lower triangular part of \(A\)) gets overwritten with \(L^T\).
If info is not zero the factorization could not be computed as the matrix was singular.
cerr << "S is singular" << endl;
return info;
}
Compute \(x = S^{-1}b\). Vector \(b\) gets overwritten with the solution vector \(x\).
cout << "inv(S)*b = " << b << endl;
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-potrs lapack-potrs.cc
Run
$shell> cd flens/examples $shell> ./lapack-potrs S = 2 1 0 1 2 1 0 1 2 b = 1 2 3 inv(S)*b = 0.5 0 1.5
Example with Complex Numbers
Example Code
#include <iostream>
using namespace flens;
using namespace std;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(3, 3);
DenseVector<Array<Complex> > b(3);
A = 2, I, 0,
0, 2, 1.+I,
0, 0, 3;
b = 1,
2,
3;
auto H = A.upper().hermitian();
cout << "H = " << H << endl;
cout << "b = " << b << endl;
int info = lapack::potrf(H);
if (info!=0) {
cerr << "H is singular" << endl;
return info;
}
lapack::potrs(H, b);
cout << "inv(H)*b = " << b << endl;
return 0;
}
Comments on Example Code
#include <iostream>
using namespace flens;
using namespace std;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(3, 3);
DenseVector<Array<Complex> > b(3);
Setup the raw data.
0, 2, 1.+I,
0, 0, 3;
b = 1,
2,
3;
\(H\) is a symmetric matrix view constructed from the upper triangular part of \(A\). Note that \(H\) only references data from \(A\).
cout << "H = " << H << endl;
cout << "b = " << b << endl;
Computes the Cholesky factorization \(H = L L^T\) where \(L\) is upper triangular. Matrix \(H\) (i.e. the lower triangular part of \(A\)) gets overwritten with \(L^T\).
If info is not zero the factorization could not be computed as the matrix was singular.
cerr << "H is singular" << endl;
return info;
}
Compute \(x = H^{-1}b\). Vector \(b\) gets overwritten with the solution vector \(x\).
cout << "inv(H)*b = " << b << endl;
return 0;
}
Compile
Note that an external LAPACK implementation is needed for the complex variant of lapack::potrs
$shell> cd flens/examples $shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-potrs lapack-complex-potrs.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-potrs H = (2,0) (0,1) (0,0) (0,-1) (2,0) (1,1) (0,-0) (1,-1) (3,0) b = (1,0) (2,0) (3,0) inv(H)*b = (0.2,-0.6) (1.2,-0.6) (0.8,0.6)