int n1 = 20;
mesh Th = square(n1,n1);
fespace Xh(Th,P1);
Xh u,v,delta,u0,w;
int n = Xh.ndof;
real mu = 1e-2;
real rho = 1;
real tol = 1e-9;
real[int] ret(n),rhs(n),deltav(n);
u0 = x*(x-1)*y*(y-1); // initial guess
// without preconditioner
func real[int] F(real[int] & u)
{
  Xh uloc;
  uloc[] = u;
  varf Fvar(unused,v) = int2d(Th)(-mu*(dx(uloc)*dx(v)+dy(uloc)*dy(v))
				  -rho*(uloc^2*v)+v)
    +on(1,2,3,4,unused=0);
  ret = Fvar(0,Xh);
  return ret;
}
func real[int] JF(real[int] & delta)
{
  Xh deltaloc;
  deltaloc[] = delta;
  varf JFvar(unused,v) = int2d(Th)(mu*(dx(deltaloc)*dx(v)+dy(deltaloc)*dy(v))
				   +2*rho*(deltaloc*u*v))
    +on(1,2,3,4,unused=0);
  ret = JFvar(0,Xh);
  return ret;
}
real cpu = clock();
u = u0;
rhs = F(u[]);
deltav = 0; // initial guess
LinearCG(JF,deltav,rhs); 
cout << "deltav " << deltav.linfty << endl;
while (deltav.linfty > tol)
  {
    u[] = u[]+deltav;
    rhs = F(u[]);
    LinearCG(JF,deltav,rhs);
    cout << "deltav " << deltav.linfty << endl;
  }
cout << "norm(u) = " << u[].l2 << " CPU s. " << clock()-cpu << endl;
// with preconditioner
cpu = clock();
varf AR0(delta,v) = int2d(Th)(mu*(dx(delta)*dx(v)+dy(delta)*dy(v))
			      +rho*(delta*u0*v))
  +on(1,2,3,4,delta=0);
matrix AR0mat = AR0(Xh,Xh,solver=Cholesky,factorize=1);
func real [int] P(real[int] & w)
{
  ret = AR0mat^-1*w;
  return ret;
}
u = u0;
rhs = F(u[]);
deltav = 0; // initial guess
LinearCG(JF,deltav,rhs,precon=P);
cout << "deltav " << deltav.linfty << endl;
while (deltav.linfty > tol)
  {
    u[] = u[]+deltav;
    rhs = F(u[]);
    LinearCG(JF,deltav,rhs,precon=P); 
    cout << "deltav " << deltav.linfty << endl;
  }
cout << "norm(u) = " << u[].l2 << " CPU s. " << clock()-cpu << endl;