Home
Browse Files

Content

Solving Systems of Linear Equations

In this example we solve a system of linear equations \(Ax = b\) were the coefficient matrix is general, i.e. not necessarily symmetric. We solve the system with lapack::sv whcih is the FLENS port of LAPACK's dgesv.

Example Code

flens/examples/lapack-gesv.cc
#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    lapack::sv(A, piv, b);

    cerr << "x = " << b << endl;
}

Comments on Example Code

#include <iostream>

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{

Define some convenient typedefs for the matrix/vector types of our system of linear equations.

    typedef GeMatrix<FullStorage<T> >           Matrix;
    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.

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

Set up the baby problem ...

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3;

    b = 20,
       -33,
       -43,
        49;

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

And solve it with lapack::sv

    lapack::sv(A, piv, b);

    cerr << "x = " << b << endl;
}

Compile

$shell> cd flens/examples                                                      
$shell> clang++ -std=c++0x -Wall -I../.. -o lapack-gesv lapack-gesv.cc

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-gesv                                                          
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
Home