// different ways to solve the same problem in the same way
border Gamma(t = 0,2*pi) { x = cos(t); y = sin(t); }
mesh Th = buildmesh(Gamma(300));
fespace Vh(Th,P1); 
int n = Vh.ndof;
real[int] uv(n),ellv(n);
Vh u,v,f = 1;
// 1) Standard Poisson, with CG solver
solve Poisson(u,v,solver = CG) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
  -int2d(Th)(f*v)+on(Gamma,u = 1);
cout << "Standard method: " << u[].l2 << endl;
// 2) set the problem and later solve
problem Poisson1(u,v,solver = CG) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
  -int2d(Th)(f*v)+on(Gamma,u = 1);
Poisson1;
cout << "Set and solve later: " << u[].l2 << endl;
// 3) now I try to solve explicitely the linear system
// bilinear form
varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
  +on(Gamma,u = 1);
// compute once and for all the associated matrix
matrix A = a(Vh,Vh,solver=CG); // specify the solver
// the right hand side: it is a bilinear form with a dummy argument
varf ell(unused,v) = int2d(Th)(f*v)+on(Gamma,unused = 1);
ellv = ell(0,Vh);
// it is possbile, for small numbers, to display the matrix and the ellv
// cout << A << endl;
// cout << ellv << endl;
uv = A^(-1)*ellv; // no inversion of A, CG in this case
cout << "Explict linear system solution: " << uv.l2 << endl;
// 4) now LinearCG
// now the func that applies the matrix to a vector
func real[int] a2u(real[int] & u)
{
  real[int] r = A*u;
  return r;
}
// diagonal (Jacobi) preconditioner
real[int] prec(n);
// extract the diagonal of A
prec = A.diag;
// the func that applies the preconditioner
func real[int] P(real[int] & u)
{
//  u ./= prec;
//  return u;
    real[int] r = u./prec;
    return r;
}
// initial vector for LinearCG (it has to satisfy bc)
// trick to force bc
varf boundary(u,v) = on(1,u = 1);
real[int] onboundary = boundary(0,Vh);
uv = 0;
// enforce bc
uv = onboundary ? 1. : uv; // try to comment this
LinearCG(a2u,uv,ellv,nbiter=1000,precon=P);
cout << "LinearCG: " << uv.l2 << endl;
// now I try without the costruction of the matrix. Of course, it is slower
func real[int] newa2u(real[int] & u)
{
  Vh uloc;
  uloc[] = u;
  varf newa(unused,v) = int2d(Th)(dx(uloc)*dx(v)+dy(uloc)*dy(v))
  +on(Gamma,unused = 1); // you can use or not bc
  // if you omit bc, then you NEED the preconditioner
  real[int] r = newa(0,Vh);
  return r;
}
uv = 0;
// enforce bc
uv = onboundary ? 1. : uv; // try to comment this
LinearCG(newa2u,uv,ellv,nbiter=1000);
cout << "LinearCG without matrix: " << uv.l2 << endl;