subroutine bicgstab2 (okprint,l, n, x, rhs, matvec, nonzero, tol, & typestop,mxmv, work, ldw, info) C C Simple BiCGstab(\ell) iterative method, \ell <= 2 C By M.A.Botchev, Jan.'98 C report bugs to botchev@cwi.nl or botchev@excite.com C C Copyright (c) 1998 by M.A.Botchev C Permission to copy all or part of this work is granted, C provided that the copies are not made or distributed C for resale, and that the copyright notice and this C notice are retained. C C This is the "vanilla" version of BiCGstab(\ell) as described C in PhD thesis of D.R.Fokkema, Chapter 3. It includes two enhancements C to BiCGstab(\ell) proposed by G.Sleijpen and H.van der Vorst in C 1) G.Sleijpen and H.van der Vorst "Maintaining convergence C properties of BiCGstab methods in finite precision arithmetic", C Numerical Algorithms, 10, 1995, pp.203-223 C 2) G.Sleijpen and H.van der Vorst "Reliable updated residuals in C hybrid BiCG methods", Computing, 56, 1996, pp.141-163 C C {{ This code based on: C subroutine bistbl v1.0 1995 C C Copyright (c) 1995 by D.R. Fokkema. C Permission to copy all or part of this work is granted, C provided that the copies are not made or distributed C for resale, and that the copyright notice and this C notice are retained. }} C C okprint == (input) LOGICAL. If okprint=.true. residual norm C will be printed to * C l == (input) INTEGER BiCGstab's dimension <= 2 C Set l=2 for highly nonsymmetric problems C n == (input) INTEGER size of the system to solve C x == (input/output) DOUBLE PRECISION array dimension n C initial guess on input, solution on output C rhs == (input) DOUBLE PRECISION array dimension n C right-hand side (rhs) vector C matvec == (input) EXTERNAL name of matrix vector subroutine C to deliver y:=A*x by CALL matvec(n,x,y) C nonzero == (input) LOGICAL tells C BiCGstab(\ell) if initial guess in x is zero or not. C If nonzero is .FALSE., one MATVEC call is saved. C tol == (input/output) DOUBLE PRECISION tolerance for all possible C stopping criteria (see the 'typestop' parameter) C On output, if info=0 or 1, tol is actually achieved C residual reduction or whatever (see the 'typestop' parameter) C typestop== (input) CHARACTER*3 stopping criterion (||.|| denotes C the 2-norm): C typestop='rel' -- relative stopping crit.: ||res|| < tol*||res0|| C typestop='abs' -- absolute stopping crit.: ||res||= n*(2*l+5) C info == (output) INTEGER. info = 0 in case of normal computations C and C info = -m (<0) - means paramater number m has an illegal value C info = 1 - means no convergence achieved (stopping criterion C is not fulfilled) C info = 2 - means breakdown of the algorithm (try to enlarge C parameter l=\ell to get rid of this) C ---------------------------------------------------------- implicit none external matvec integer l, n, mxmv, ldw, info double precision x(n), rhs(n), tol logical okprint,nonzero character*3 typestop double precision work(n,5+2*l) C ----------------------------------------- integer lmax parameter(lmax=2) double precision rwork(lmax+1,3+2*(lmax+1)) logical GoOn, rcmp, xpdt integer ii, i1, jj, kk, nmv double precision alpha,beta,hatgamma,kappa0, kappal,maxval1, & mxnrmr,mxnrmx,omega,rho0,rho1,rnrm0,rnrm,rnrmMax, & sigma,sum1,varrho integer z, zz, y0, yl, y integer rr, r, u, xp, bp double precision zero, one, delta parameter(zero=0d0,one=1d0,delta=1d-2) info = 0 if (l.gt.lmax .or. l.lt.1) info = -2 if (tol.le.zero) info = -8 if (mxmv.lt.0) info = -10 rr = 1 r = rr+1 u = r+(l+1) xp = u+(l+1) bp = xp+1 if (bp*n.gt.ldw) info = -12 z = 1 zz = z+(l+1) y0 = zz+(l+1) yl = y0+1 y = yl+1 if (info.ne.0) return C Initialize first residual if (nonzero) then call matvec (n, x, work(1,r) ) do ii=1,n work(ii,r) = rhs(ii) - work(ii,r) enddo nmv = 1 else do ii=1,n work(ii,r) = rhs(ii) enddo nmv = 0 endif C Initialize iteration loop sum1 = zero do ii=1,n work(ii,rr) = work(ii,r) work(ii,bp) = work(ii,r) work(ii,xp) = x(ii) x(ii) = zero sum1=sum1+ work(ii,r)**2 enddo rnrm0 = sqrt( sum1 ) rnrm = rnrm0 if (typestop.eq.'max') then maxval1=zero do ii=1,n maxval1=max( maxval1, abs( work(ii,r) ) ) enddo rnrmMax = maxval1 endif mxnrmx = rnrm0 mxnrmr = rnrm0 rcmp = .false. xpdt = .false. alpha = zero omega = one sigma = one rho0 = one C Iterate if(typestop.eq.'rel')then GoOn = rnrm.ge.tol*rnrm0 .and. nmv.lt.mxmv else if(typestop.eq.'abs')then GoOn = rnrm.ge.tol .and. nmv.lt.mxmv else if(typestop.eq.'max')then GoOn = rnrmMax.ge.tol .and. nmv.lt.mxmv else info = -9 return end if do while (GoOn) C ===================== C --- The BiCG part --- C ===================== rho0 = -omega*rho0 do kk=1,l sum1 = zero do ii=1,n sum1=sum1+ work(ii,rr)*work(ii,r+kk-1) enddo rho1 = sum1 if (rho0.eq.zero) then info = 2 return endif beta = alpha*(rho1/rho0) rho0 = rho1 do jj=0,kk-1 do ii=1,n work(ii,u+jj) = work(ii,r+jj) - beta*work(ii,u+jj) enddo enddo call matvec(n, work(1,u+kk-1), work(1,u+kk)) nmv = nmv+1 sum1 = zero do ii=1,n sum1=sum1+ work(ii,rr)*work(ii,u+kk) enddo sigma = sum1 if (sigma.eq.zero) then info = 2 return endif alpha = rho1/sigma do ii=1,n x(ii) = alpha*work(ii,u) + x(ii) enddo do jj=0,kk-1 do ii=1,n work(ii,r+jj) = -alpha*work(ii,u+jj+1) + work(ii,r+jj) enddo enddo call matvec (n, work(1,r+kk-1), work(1,r+kk)) nmv = nmv+1 sum1 = zero do ii=1,n sum1=sum1+ work(ii,r)**2 enddo rnrm = sqrt( sum1 ) if (typestop.eq.'max') then maxval1=zero do ii=1,n maxval1=max( maxval1, abs( work(ii,r) ) ) enddo rnrmMax = maxval1 endif mxnrmx = max (mxnrmx, rnrm) mxnrmr = max (mxnrmr, rnrm) enddo C ================================== C --- The convex polynomial part --- C ================================== C Z = R'R do i1=1,l+1 do jj=i1-1,l sum1 = zero do ii=1,n sum1=sum1+ work(ii,r+jj)*work(ii,r+i1-1) enddo rwork(jj+1,z+i1-1) = sum1 rwork(z+i1-1,jj+1) = rwork(jj+1,z+i1-1) enddo enddo do i1=zz,zz+l do ii=1,l+1 rwork(ii,i1) = rwork(ii,i1+(z-zz)) enddo enddo C tilde r0 and tilde rl (small vectors) rwork(1,y0) = -one rwork(2,y0) = rwork(2,z) / rwork(2,zz+1) rwork(l+1,y0) = zero rwork(1,yl) = zero rwork(2,yl) = rwork(2,z+l) / rwork(2,zz+1) rwork(l+1,yl) = -one C Convex combination do ii=1,l+1 rwork(ii,y) = zero enddo do jj=1,l+1 do ii=1,l+1 rwork(ii,y) = rwork(ii,y) + rwork(jj,yl)* & rwork(ii,z+jj-1) enddo enddo sum1 = zero do ii=1,l+1 sum1=sum1+ rwork(ii,yl)*rwork(ii,y) enddo kappal = sqrt( sum1 ) do ii=1,l+1 rwork(ii,y) = zero enddo do jj=1,l+1 do ii=1,l+1 rwork(ii,y) = rwork(ii,y) + rwork(jj,y0)* & rwork(ii,z+jj-1) enddo enddo sum1 = zero do ii=1,l+1 sum1=sum1+ rwork(ii,y0)*rwork(ii,y) enddo kappa0 = sqrt( sum1 ) sum1 = zero do ii=1,l+1 sum1=sum1+ rwork(ii,yl)*rwork(ii,y) enddo varrho = sum1 varrho = varrho / (kappa0*kappal) hatgamma = sign(1d0,varrho)*max(abs(varrho),7d-1) * & (kappa0/kappal) do ii=1,l+1 rwork(ii,y0) = -hatgamma*rwork(ii,yl) + rwork(ii,y0) enddo C Update omega = rwork(l+1,y0) do jj=1,l do ii=1,n work(ii,u) = work(ii,u) - rwork(1+jj,y0)*work(ii,u+jj) x(ii) = x(ii) + rwork(1+jj,y0)*work(ii,r+jj-1) work(ii,r) = work(ii,r) - rwork(1+jj,y0)*work(ii,r+jj) enddo enddo do ii=1,l+1 rwork(ii,y) = zero enddo do jj=1,l+1 do ii=1,l+1 rwork(ii,y) = rwork(ii,y) + rwork(jj,y0)* & rwork(ii,z+jj-1) enddo enddo sum1 = zero do ii=1,l+1 sum1=sum1+ rwork(ii,y0)*rwork(ii,y) enddo rnrm = sqrt( sum1 ) C ================================ C --- The reliable update part --- C ================================ mxnrmx = max (mxnrmx, rnrm) mxnrmr = max (mxnrmr, rnrm) xpdt = (rnrm.lt.delta*rnrm0.and.rnrm0.lt.mxnrmx) rcmp = ((rnrm.lt.delta*mxnrmr.and.rnrm0.lt.mxnrmr) .or.xpdt) if (rcmp) then call matvec (n, x, work(1,r) ) nmv = nmv + 1 do ii=1,n work(ii,r) = work(ii,bp) - work(ii,r) enddo mxnrmr = rnrm if (xpdt) then do ii=1,n work(ii,xp) = x(ii) + work(ii,xp) x(ii) = zero work(ii,bp) = work(ii,r) enddo mxnrmx = rnrm endif endif if(typestop.eq.'rel')then GoOn = rnrm.ge.tol*rnrm0 .and. nmv.lt.mxmv if(okprint) print *, nmv,' ',rnrm/rnrm0 else if(typestop.eq.'abs')then GoOn = rnrm.ge.tol .and. nmv.lt.mxmv if(okprint) print *, nmv,' ',rnrm else if(typestop.eq.'max')then maxval1=zero do ii=1,n maxval1=max(maxval1, abs( work(ii,r) ) ) enddo rnrmMax = maxval1 GoOn = rnrmMax.ge.tol .and. nmv.lt.mxmv if(okprint) print *, nmv,' ',rnrmMax end if enddo C ========================= C --- End of iterations --- C ========================= do ii=1,n x(ii) = work(ii,xp) + x(ii) enddo C --------------------- One matvec can be saved by commenting out this: C (this is to compute the true residual) call matvec (n, x, work(1,r) ) do ii=1,n work(ii,r) = rhs(ii) - work(ii,r) enddo sum1 = zero do ii=1,n sum1=sum1+ work(ii,r)**2 enddo rnrm = sqrt( sum1 ) C --------------------- One matvec can be saved by commenting out this^ if(typestop.eq.'rel')then if (rnrm.gt.tol*rnrm0) info = 1 tol = rnrm/rnrm0 else if(typestop.eq.'abs')then if (rnrm.gt.tol) info = 1 tol = rnrm else if(typestop.eq.'max')then if (rnrmMax.gt.tol) info = 1 tol = rnrmMax end if mxmv = nmv return end