About the functions LinearCG

General Discussion
8

Dear all,
thank you for all your answers, and for your all your help. Sorry to be late to get back to you. I have been contacted out of the forum by Pierre Jolivet (http://jolivet.perso.enseeiht.fr), who kindly explained to me the main stakes of the LinearCG function, and why my codes were not running. With his authorization, here are the main lines of his explanations.

In the piece of code

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
+ int2d(Th)(f()v)
+ on(1,u=1);
matrix A = a(Vh,Vh,solver=CG);
Vh b,u;
b[] = a(0,Vh,tgv=tgv);
u[] = A^-1b[];

it is implicitly understood that the CG solver is used for the calculation of u[], together with a Jacobi (diagonal) preconditioner. The latter is necessary to “kill” the 1e30 values on the diagonal of A and the rows of b corresponding to Dirichlet boundary conditions in doing the matrix-vector products featured in CG.

Then I tried exactly the same thing, with a Jacobi preconditioner:

real tgv = 1.e30;

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
+ int2d(Th)(f()*v)
+ on(1,u=1);

matrix A = a(Vh,Vh,tgv=tgv);
Vh b,u;
b[] = a(0,Vh,tgv=tgv);

/* Using conjugate gradient /
func real[int] opA(real[int] &p) {
pp[] = Ap;
return(pp[]);
}

/* Preconditioner */
func real[int] opQ(real[int] &p) {
for (int i=0; i<p.n; i++) {
p[i] = p[i] / A(i,i);
}
return§;
}
u[] = 0;
LinearCG(opA,u[],b[],precon=opQ,eps=1.e-6,nbiter=1000);

and it still does not work properly, because the CG algorithm stops after one iteration! The reason is that the initial residual, which is used in the CG algorithm to set the tolerance threshold for convergence is calculated, based upon the provided values of A, u and b, without preconditioning. Again, because of the tgv values, the threshold has an absurd value. To fix this, it is desirable to initialize u no longer at 0 but at 0, together with the correct values at the Dirichlet boundary (and 0 elsewhere), so that no tgv appears in the initial residual Au-b.

Hence, the token

u[] = a(0,Vh,tgv=-1);

before calling LinearCG makes things work fine.
The tgv=-1 means for the vector assembly that only boundary conditions are set without the tgv. When it comes to the matrix assembly, tgv=-1 or tgv=-2 imply that nonsymmetric or symmetric elimination are used, as described in the following tutorial by P. Jolivet:

http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/#pf14

Eventually, Pierre advocates rather the use of LinearGMRES together with a nonsymmetric elimination of the Dirichlet degrees of freedom instead of the use of LinearCG together with symmetric elimination or use of the “tgv trick”.

Thank you very much (and to Pierre) for your help, in any event. This really made me understand a lot of things about the concrete implementation of routines in FreeFem!

1 Like