# Sparse Poisson problem in eigen

Back 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;

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.

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 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;

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;
}
```