#include "dolfin.h"

using namespace dolfin;

// PPDE DOLFIN example #2: matrices, Krylov solvers

int main()
{
  dolfin_set("output", "plain text");

  int N = 10;

  Matrix A(N, N);
  Vector b(N);
  Vector x(N);

  cout << "A:" << endl;
  A.show();

  cout << "b:" << endl;
  b.show();

  for(int i = 0; i < N; i++)
  {
    A(i, i) = 2.0;
  }

  for(int i = 0; i < N - 1; i++)
  {
    A(i, i + 1) = -1.0;
    A(i + 1, i) = -1.0;
  }

  for(int i = 0; i < N; i++)
  {
    b(i) = 1.0;
    x(i) = 0;
  }

  cout << "A:" << endl;
  A.show();

  cout << "b:" << endl;
  b.show();

  cout << "x:" << endl;
  x.show();

  KrylovSolver solver;
  solver.solve(A, x, b);

  cout << "x:" << endl;
  x.show();

  Vector r(N);

  A.mult(x, r);
  r -= b;

  cout << "r:" << endl;
  r.show();

  cout << "r norm: " << r.norm() << endl;
}