Sparse Poisson problem in eigen

OwlgenBack to scientific computing. Lately, I have been using the Eigen libraries for linear algebra. They were recommended by the CGAL project, and they indeed share a common philosophy. Thanks to the rather high C++ level they can accomplish this sort of magic:


  int n = 100;
  VectorXd x(n), b(n);
  SpMat A(n,n);

  fill_A(A);

  fill_b(b);

  // solve Ax = b
  Eigen::BiCGSTAB<SpMat> solver;
  //ConjugateGradient<SpMat> solver;

  solver.compute(A);

  x = solver.solveWithGuess(b,x0);

Notice that A is a sparse matrix! I am next describing how to use this in order to solve the 1D Poisson equation.

jau70

Now, A is a sparse matrix. Care must be taken when filling it, and when addressing its contents. Here is an efficient way of filling it up using a std::vector of triplets. These are (int i, int j, T b)  with T= float or double typically (or complex), giving the value b of the element at row i, column j. Notice I fill it with the celebrated finite difference expression for the Laplacian operator.

void fill_A(SpMat& A) {

 typedef Eigen::Triplet<double> T;

 int n=A.rows();

 std::vector<T> aa; // list of non-zeros coefficients

 double dx2 = double(n)*double(n);

 for(int i=0 ; i < n ; i++) {

 if(i>0) aa.push_back( T(i,i-1, dx2 ));
 // int im1= (i-1 + n) % n;
 // aa.push_back( T(i,im1, dx2 ));

 aa.push_back( T(i,i, -2 * dx2 ) );

 // int ip1= (i+1 + n) % n;
 // aa.push_back( T(i,ip1, dx2 ));

 // cout << im1 << " " << i << " " << ip1 << endl;

 if(i<n-1) aa.push_back( T(i,i+1, dx2 ));

 }

 A.setFromTriplets(aa.begin(), aa.end());

 return;
}

This is how one fills the source term (right hand side of the Poisson equation). It is a regular vector (in the algebra sense), not a sparse one (i.e. it is a “dense” vector):


void fill_b(VectorXd& b) {
 int n=b.size();
 for(int i=0 ; i<n ; i++) {

 double x=double(i+1)/double(n+1);

 b(i)= std::sin(2*M_PI*x);

 std::cout << x << " " << b(i) << std::endl;
 }
 return;
}

We can also provide a vector as an initial guess, since in this case we know
the analytic solution:


void solution(VectorXd& b) {
 int n=b.size();
 for(int i=0 ; i<n ; i++) {
 double x=double(i+1)/double(n+1);

 b(i)= - std::sin(2*M_PI*x) / (4*M_PI*M_PI);
 }
 return;
}

FYI, this is the complete code, with ugly comments and everything. Some comments would apply for another canonical solution, f(x)= x (1-x), with a constant source term equal to 2.


#include <Eigen/IterativeLinearSolvers>
#include <iostream>
#include <vector>

using Eigen::VectorXi;
using Eigen::VectorXd;
using Eigen::SparseMatrix;
using Eigen::ConjugateGradient;

using std::cout;
using std::endl;

const Eigen::IOFormat OctaveFmt(Eigen::StreamPrecision, 0, ", ", ";\n", "", "", "[", "];");

typedef SparseMatrix<double> SpMat;

void fill_A(SpMat& A) {

 typedef Eigen::Triplet<double> T;

 // A.reserve(VectorXi::Constant(n,3));

 std::cout << " Filling A " << std::endl;

 int n=A.rows();

 std::vector<T> aa; // list of non-zeros coefficients

 double dx2 = double(n)*double(n);

 for(int i=0 ; i < n ; i++) {

 if(i>0) aa.push_back( T(i,i-1, dx2 ));

 // int im1= (i-1 + n) % n;
 // aa.push_back( T(i,im1, dx2 ));

 aa.push_back( T(i,i, -2 * dx2 ) );

 // int ip1= (i+1 + n) % n;
 // aa.push_back( T(i,ip1, dx2 ));

 // cout << im1 << " " << i << " " << ip1 << endl;

 if(i<n-1) aa.push_back( T(i,i+1, dx2 ));

 }

 A.setFromTriplets(aa.begin(), aa.end());

 // std::cout << i << std::endl;
// int im1= (i-1 + n) % n;
// int ip1= (i+1 + n) % n;

// std::cout << i << " " << im1 << std::endl;
// std::cout << i << " " << ip1 << std::endl;
// A.insert(i,ip1)= dx2;

 std::cout << " Filled A " << std::endl;

 // for(int i=0 ; i < n ; i++) {

 // if(i>0) cout << A.coeffRef(i,i-1) << " ";

 // cout << A.coeffRef(i,i) << " ";

 // if(i<n-1) cout << A.coeffRef(i,i + 1) << " ";

 // cout << endl;

 // }

 return;
}

void fill_b(VectorXd& b) {
 int n=b.size();
 for(int i=0 ; i<n ; i++) {
 // b(i)=2;

 double x=double(i+1)/double(n+1);

 b(i)= std::sin(2*M_PI*x);

 std::cout << x << " " << b(i) << std::endl;
 }
 return;
}

void solution(VectorXd& b) {
 int n=b.size();
 for(int i=0 ; i<n ; i++) {
 double x=double(i+1)/double(n+1);

 b(i)= - std::sin(2*M_PI*x) / (4*M_PI*M_PI);

 // b(i)= x*(1-x);
 }
 return;
}

int main(void) {

 int n = 100;
 VectorXd x(n), b(n);
 SpMat A(n,n);

 // fill A
 fill_A(A);

 // fill b

 fill_b(b);

 // std::cout << "A= " << A.format(OctaveFmt) << std::endl;
 std::cout << "b= " << b.format(OctaveFmt) << std::endl;

 VectorXd x0(n);
 solution(x0);

 std::cout << "x0= " << x0.format(OctaveFmt) << std::endl;

 VectorXd bb = A*x0;

 std::cout << "bb= " << bb.format(OctaveFmt) << std::endl;

 // solve Ax = b
 Eigen::BiCGSTAB<SpMat> solver;
 //ConjugateGradient<SpMat> solver;

 solver.compute(A);
 if(solver.info()!=Eigen::Success) {
 // decomposition failed
 std::cout << "Failure decomposing\n";
 return 1;
 }

 x = solver.solveWithGuess(b,x0);
 if(solver.info()!=Eigen::Success) {
 // solving failed
 std::cout << "Failure solving\n";
 return 1;
 }
 // solve for another right hand side:
 // x1 = solver.solve(b1);

 std::cout << "#iterations: " << solver.iterations() << std::endl;
 std::cout << "estimated error: " << solver.error() << std::endl;

 std::cout << "x = " << x.format(OctaveFmt) << std::endl;

 // update b, and solve again
 //* x = cg.solve(b);

 return 0;
}

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s