[eigen] initial sparse PCG solver

[Daniel Lowengrub <lowdanie@xxxxxxxxx>]

Hi,
I started working on a preconditioned conjugate gradient solver for
the Sparse module.  Before working on the interface and eigen
integration I decided to try and get the basic functionality setup.  I
attached a file with the implementation of the algorithm and a small
test program.  The test program uses the simple jacobi preconditioner.
 Obviously better preconditioners will have to be implemented in order
for the PCG algorithm to be feasible, but for now I'd like to focus on
the algorithm itself.
Any feedback on the implementation would be appreciated.
The implementation used here is based on this paper:
www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf

Daniel

#include <iostream>
#include <Eigen/Sparse>
#include <Eigen/Core>

#define EPSILON 0.01

typedef Eigen::Matrix<double, 10, 10> Matrix10d;
typedef Eigen::Matrix<double, 10, 1> Vector10d;
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXd;

typedef Eigen::SparseMatrix<double> SparseMatd;
typedef Eigen::SparseVector<double> SparseVecd;

//---------------------------------------------------------
// Iteratively solves Ax = b using the preconditioner M
//--------------------------------------------------------- 
SparseVecd conjugateGradient(SparseMatd& A, SparseMatd& invM, SparseVecd& b, SparseVecd& x0, int maxIter)
{
  SparseVecd x = x0; //initial approximation
  SparseVecd r = b - A * x0; //initial residual
  SparseVecd dir = invM * r; //initial search direction
  SparseVecd s, q;
  double absNew = r.dot(dir); //the square of the absolute value of r scaled by invM
  double absInit = absNew;  //the initial absolute value
  double absOld;
  double alpha, beta; //some values used later in the algorithm
  int i = 0;

  while ((i < maxIter) && (absNew > EPSILON*EPSILON*absInit))
    {
      q = A * dir;  //store this handy value
      alpha = absNew / dir.dot(q); //the amount we travel on dir
      x += alpha * dir; //update solution
      r -= alpha * q; //update residue
      s = invM * r; //another useful value
      absOld = absNew; 
      absNew = r.dot(s); //update the absolute value of r
      beta = absNew / absOld; //calculate the Gram-Schmidit value used to create the new search direction
      dir = s + beta * dir; //update search direction
      i++;
    }
  return x; //return the closest approximation
}

int main()
{  
  Matrix10d A, E, J;
  Vector10d d, x;
  VectorXd b, x0;
  SparseMatd sA, sJ;
  SparseVecd sb, sx0;

  E = Matrix10d::Random();
  E = E.part<Eigen::Upper>();
  A = E.transpose() * E; //construct a random positive definite matrix
  
  d = A.diagonal();
  d = d.array().inverse();
  J = Matrix10d::Identity();
  J.diagonal() = d; //set J to the inverse of the diagonal of A as a simple preconditioner

  b = Vector10d::Random();
  x0 = Vector10d::Zero(); //start from 0 since we hae no information

  //create sparse matricies to pass to the CG function
  sA = A.sparseView();
  sJ = J.sparseView();
  sb = b.sparseView();
  sx0 = x0.sparseView();

  x = (conjugateGradient(sA, sJ, sb, sx0, 10)).toDense();

  Vector10d residue = b - A*x;
  std::cout <<"the residue norm is: "<< std::endl << residue.norm() <<std::endl;
  std::cout <<"the answer is: "<< std::endl << x <<std::endl;

  return 1;
}