Demo: Using a linear system direct solver

As a first example we present a simple program that solves a linear system of equations with a tridiagonal matrix. We use for this the Gauss elimination method.

  • As usual, we start by including the principal header file and defining the namespace . We start the program by reading the argument as the system's size. Next, we instantiate the matrix as a tridiagonal one (instance of TrMatrix) and vectors b and x that stand for the right-hand side and the solution vector respectively. 

    #include "OFELI.h"
using namespace OFELI;
                          
int main(int argc, char *argv[])
{
  size_t n=atoi(argv[1]);
  TrMatrix<double> A(n);
  Vect<double> b(n), x(n);
}

  • We next initialize the matrix and the right-hand side vector. For this example, the solution vector satisfies x(i)=i: 

       for (size_t i=2; i<n; i++) {
     A(i,i)   =  2.0;
     A(i,i-1) = -1.0;
     A(i,i+1) = -1.0;
  }
  A(1,1)   =  2.0;
  A(1,2)   = -1.0;
  A(n,n-1) = -1.0;
  A(n,n)   =  2.0;
  b = 0.0;
  b(n) = n+1;

  • To solve this linear system we define an instance of class LinearSolver. The constructor uses the linear system data. Note that the default constructor can be used and then the member functions setMatrix, setRHS, setSolution. We can then define the solver by using the enumerated value DIRECT_SOLVER and solve the system. 

      LinearSolver<double> ls(A,b,x);
  ls.setSolver(DIRECT_SOLVER);
  ls.solve();

  • We end by printing the solution. 

       cout << "Solution:\n" << x << endl; 
   return 0;
}