pbcgstab_vdv

 PBCGSTAB_VDV Preconditioned BiCGStab method with FEM preconditioners

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

 preconditioned BiCGStab method

    P^{-1}Ax=P^{-1}b with A and P SPD
    The linear system P z=r is solved by LU 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

0001 function [u,iter,err]=pbcgstab_vdv(A, b, u, tol, maxit ,P,preco);
0002 % PBCGSTAB_VDV Preconditioned BiCGStab method with FEM preconditioners
0003 %
0004 %  [x,iter,res]=pbcgstab_vdv(A, b, x0, tol, itmax,P,preco)
0005 %
0006 % preconditioned BiCGStab method
0007 %
0008 %    P^{-1}Ax=P^{-1}b with A and P SPD
0009 %    The linear system P z=r is solved by LU 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 r=b-A*u;
0043 rtilde0=r;
0044 normar0=r'*r;
0045 iter=0; err=1;
0046 alpha=1;rho=1;rhom1=1;w=1;
0047  if preco~=0
0048  [L,U,PP]=lu(P);
0049  end
0050 
0051 while iter< maxit & err> tol
0052 rho=rtilde0'*r;
0053 if iter==0
0054 p=r;
0055 else
0056 beta=(rho/rhom1)*(alpha/w);
0057 p=r+beta*(p-w*v);
0058 end
0059 if preco==0
0060 pc=p;
0061 else
0062 pc=U\(L\(PP*p));
0063 end
0064 
0065 v=A*pc;
0066 alpha=rho/(v'*rtilde0);
0067 s=r-alpha*v;
0068 if ( norm(s) < tol ),                          % early convergence check
0069         u = u + alpha*pc;
0070         err = norm( s ) / normar0;
0071         break;
0072      end
0073 if preco==0
0074 sc=s;
0075 else
0076 sc=U\(L\(PP*s));
0077 end
0078 t=A*sc;
0079 w=(t'*s)/(t'*t);
0080 u=u+alpha*pc+w*sc;
0081 r=s-w*t;
0082 rhom1=rho;
0083 iter=iter+1;
0084 err=sqrt(r'*r)/normar0;
0085 end