Content |
Solving Systems of Linear Equations
In this example we will solve a system of linear equations \(Ax=b\). For this purpose we first compute the \(LU\) factorization of \(A\) with lapack::trf and then use the triangular solver lapack::trs to finish the job.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
GeMatrix<FullStorage<double> > A(4, 4);
DenseVector<Array<double> > b(4);
DenseVector<Array<int> > piv(4);
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::trf(A, piv);
lapack::trs(NoTrans, A, piv, b);
cout << "X = " << b << endl;
}
Comments on Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
GeMatrix<FullStorage<double> > A(4, 4);
DenseVector<Array<double> > b(4);
DenseVector<Array<int> > piv(4);
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
Compute the \(PLU = A\) factorization with lapack::trf
Solve the system of linear equation \(PLUx = b\) using lapack::trs. Vector \(b\) gets overwritten with vector \(x\).
cout << "X = " << b << endl;
}
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-getrs lapack-getrs.cc
Run
$shell> cd flens/examples $shell> ./lapack-getrs 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
Example with Complex Numbers
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(4, 4);
DenseVector<Array<Complex> > b(4);
DenseVector<Array<int> > piv(4);
A = 2, 3, -I, 0,
-6, -5, 0, 2.+I,
2.*I, -5, 6, -6,
4, 6.*I, 2, -3;
b = 20,
-33.*I,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::trf(A, piv);
lapack::trs(NoTrans, A, piv, b);
cout << "X = " << b << endl;
}
Comments on Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(4, 4);
DenseVector<Array<Complex> > b(4);
DenseVector<Array<int> > piv(4);
Setup matrix \(A\). Note that 2+I or 2*I would not work. That's because the STL does not define the operation int+complex<double> or int*complex<double>.
-6, -5, 0, 2.+I,
2.*I, -5, 6, -6,
4, 6.*I, 2, -3;
b = 20,
-33.*I,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
Compute the \(PLU = A\) factorization with lapack::trf
Solve the system of linear equation \(PLUx = b\) using lapack::trs. Vector \(b\) gets overwritten with vector \(x\).
cout << "X = " << b << endl;
}
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-getrs lapack-complex-getrs.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-getrs A = (2,0) (3,0) (-0,-1) (0,0) (-6,0) (-5,0) (0,0) (2,1) (0,2) (-5,0) (6,0) (-6,0) (4,0) (0,6) (2,0) (-3,0) b = (20,0) (-0,-33) (-43,0) (49,0) X = (16.5967,25.0647) (-9.15476,-5.78498) (32.7744,14.2709) (39.2151,24.624)