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