precg

 PRECG Preconditioned Conjugate Gradient method with FEM preconditioners

  [x,iter,res]=precg(A, b, x0, tol, itmax,P,preco)

 preconditioned conjugate gradient

    P^{-1}Ax=P^{-1}b with A and P SPD
    The linear system P z=r is solved by Cholesky factorization. P is
    factorized once.

 Input:  A: matrix (n x n)
         b: right hand side (n x 1)
         x0: intial datum (n x 1)
         tol: tolerance for the stopping test
         itmax: maximum number of iterations
         preco : choice of the preconditioner:
            1 : P=K_FE            (4.4.45), pag. 221, CHQZ2
                A=K_GNI
            2 : P=M_FE^{-1}K_FE   (4.4.46), pag. 221, CHQZ2
                A=M_GNI^{-1}K_GNI
            3 : P=M_FE,d^{-1}K_FE   (4.4.47), pag. 221, CHQZ2
                A=M_GNI^{-1}K_GNI
            4 : P=M_FE^{-1/2}K_FE M_FE^{-1/2}  (4.4.48), pag. 221, CHQZ2
                A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
            5 : P=KM_FE<d^{-1/2}K_FE M_FE,d^{-1/2}_FE (4.4.49), pag. 221, CHQZ2
                 A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
         P: preconditioner

 Output: x: solution (n x 1)
         iter: iterations needed to satisfy stopping test
         res: normalized 2-norm residual at the last iteration iter


 Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
                    "Spectral Methods. Fundamentals in Single Domains"
                    Springer Verlag, Berlin Heidelberg New York, 2006.

0001 function [x,iter,res]=precg(A, b, x0, tol, itmax,P,preco)
0002 % PRECG Preconditioned Conjugate Gradient method with FEM preconditioners
0003 %
0004 %  [x,iter,res]=precg(A, b, x0, tol, itmax,P,preco)
0005 %
0006 % preconditioned conjugate gradient
0007 %
0008 %    P^{-1}Ax=P^{-1}b with A and P SPD
0009 %    The linear system P z=r is solved by Cholesky factorization. P is
0010 %    factorized once.
0011 %
0012 % Input:  A: matrix (n x n)
0013 %         b: right hand side (n x 1)
0014 %         x0: intial datum (n x 1)
0015 %         tol: tolerance for the stopping test
0016 %         itmax: maximum number of iterations
0017 %         preco : choice of the preconditioner:
0018 %            1 : P=K_FE            (4.4.45), pag. 221, CHQZ2
0019 %                A=K_GNI
0020 %            2 : P=M_FE^{-1}K_FE   (4.4.46), pag. 221, CHQZ2
0021 %                A=M_GNI^{-1}K_GNI
0022 %            3 : P=M_FE,d^{-1}K_FE   (4.4.47), pag. 221, CHQZ2
0023 %                A=M_GNI^{-1}K_GNI
0024 %            4 : P=M_FE^{-1/2}K_FE M_FE^{-1/2}  (4.4.48), pag. 221, CHQZ2
0025 %                A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
0026 %            5 : P=KM_FE<d^{-1/2}K_FE M_FE,d^{-1/2}_FE (4.4.49), pag. 221, CHQZ2
0027 %                 A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
0028 %         P: preconditioner
0029 %
0030 % Output: x: solution (n x 1)
0031 %         iter: iterations needed to satisfy stopping test
0032 %         res: normalized 2-norm residual at the last iteration iter
0033 %
0034 %
0035 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0036 %                    "Spectral Methods. Fundamentals in Single Domains"
0037 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0038 
0039 %   Written by Paola Gervasio
0040 %   $Date: 2007/04/01$
0041 
0042 
0043 x=x0;
0044 res=[];
0045 r=b-A*x;
0046 if(preco==0)
0047 %
0048 z=r;
0049 elseif(preco==1 | preco==4 | preco ==5)
0050 %
0051 P1=chol(P);
0052 z=P1\((P1)'\r);
0053 
0054 end
0055 rr0=z'*r;
0056 normar0=sqrt(r'*r);
0057 p=z;
0058 err=normar0; iter=0;
0059 
0060 while  err > tol & iter < itmax
0061 
0062 v=A*p;
0063 alpha=rr0/(p'*v);
0064 x=x+alpha*p;
0065 r=r-alpha*v;
0066 if(preco==0)
0067 z=r;
0068 elseif(preco==1 | preco==4 | preco==5)
0069 z=P1\((P1)'\r);
0070 end
0071 rr=z'*r;
0072 beta=rr/rr0;
0073 p=z+beta*p;
0074 err=sqrt(r'*r)/normar0;
0075 %res=[res;err];
0076 iter=iter+1;
0077 rr0=rr;
0078 %fprintf('It. GC %d, err %13.6e \n',iter,err)
0079 end
0080 res=err;
0081 if err>tol & iter>=itmax
0082 fprintf('The method does not converge to prefixed tolerance in %d iterations\n',iter)
0083 end