c subroutine bistbl (l, n, x, b, mv, solve, tol, $ mxmv, work, ldw, rwork, ldrw, iwork, info) c 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 THIS WORK IS PROVIDED ON AN "AS IS" BASIS. THE AUTHOR c PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED, c REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS c MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE. c implicit none c c .. Parameters .. c integer l, n, mxmv, ldw, ldrw, iwork(l+1), info double precision x(n), b(n), tol double precision work(n,3+2*(l+1)), rwork(l+1,3+2*(l+1)) c c .. Matrix .. c external mv external solve c c .. Local .. c logical rcmp, xpdt integer i, j, k, nmv double precision alpha, beta, omega, rho0, rho1, sigma double precision varrho, hatgamma double precision rnrm0, rnrm double precision mxnrmx, mxnrmr, kappa0, kappal c c .. Work Aliases .. c integer z, zz, y0, yl, y integer rr, r, u, xp, bp c c .. Constants .. c double precision zero, one, delta parameter (zero = 0d0, one = 1d0, delta = 1d-2) c c .. BLAS and LAPACK .. c c subroutine daxpy c subroutine dcopy c subroutine dgemv c subroutine dgetrf c subroutine dgetrs c subroutine dlacpy c subroutine dlaset c subroutine dsymv c function ddot c function dnrm2 c double precision dnrm2, ddot c c .. Intrinsic .. c intrinsic abs, max, sign, sqrt c c =========================== c .. Executable Statements .. c =========================== c info = 0 if (l.lt.1) info = -1 if (n.lt.1) info = -2 if (tol.le.zero) info = -7 if (mxmv.lt.0) info = -8 rr = 1 r = rr+1 u = r+(l+1) xp = u+(l+1) bp = xp+1 if (bp*n.gt.ldw) info = -10 z = 1 zz = z+(l+1) y0 = zz+(l+1) yl = y0+1 y = yl+1 if (y*(l+1).gt.ldrw) info = -12 if (info.ne.0) return c c --- Initialize first residual c call mv (n, x, work(1,r)) do i=1,n work(i,r) = b(i) - work(i,r) enddo call solve (n, work(1,r)) c c --- Initialize iteration loop c nmv = 0 call dcopy (n, work(1,r), 1, work(1,rr), 1) call dcopy (n, work(1,r), 1, work(1,bp), 1) call dcopy (n, x, 1, work(1,xp), 1) call dlaset ('n', n, 1, zero, zero, x, 1) rnrm0 = dnrm2 (n, work(1,r), 1) rnrm = rnrm0 mxnrmx = rnrm0 mxnrmr = rnrm0 rcmp = .false. xpdt = .false. alpha = zero omega = one sigma = one rho0 = one c c --- Iterate c do while (rnrm.gt.tol*rnrm0 .and. nmv.lt.mxmv) c c ===================== c --- The BiCG part --- c ===================== c rho0 = -omega*rho0 do k=1,l rho1 = ddot (n, work(1,rr), 1, work(1,r+k-1), 1) if (rho0.eq.zero) then info = 2 return endif beta = alpha*(rho1/rho0) rho0 = rho1 do j=0,k-1 do i=1,n work(i,u+j) = work(i,r+j) - beta*work(i,u+j) enddo enddo call mv (n, work(1,u+k-1), work(1,u+k)) call solve (n, work(1,u+k)) nmv = nmv+1 sigma = ddot (n, work(1,rr), 1, work(1,u+k), 1) if (sigma.eq.zero) then info = 2 return endif alpha = rho1/sigma call daxpy (n, alpha, work(1,u), 1, x, 1) do j=0,k-1 call daxpy (n, (-alpha), work(1,u+j+1), 1, $ work(1,r+j), 1) enddo call mv (n, work(1,r+k-1), work(1,r+k)) call solve (n, work(1,r+k)) nmv = nmv+1 rnrm = dnrm2 (n, work(1,r), 1) mxnrmx = max (mxnrmx, rnrm) mxnrmr = max (mxnrmr, rnrm) enddo c c ================================== c --- The convex polynomial part --- c ================================== c c --- Z = R'R c do i=1,l+1 call dgemv ('t', n, l+1-(i-1), one, work(1,r+i-1), $ n, work(1,r+i-1), 1, zero, rwork(i,z+i-1), 1) call dcopy (l-(i-1), rwork(i+1,z+i-1), 1, $ rwork(i,z+i), l+1) enddo call dlacpy ('a', l+1, l+1, rwork(1,z), l+1, $ rwork(1,zz), l+1) call dgetrf (l-1, l-1, rwork(2,zz+1), l+1, $ iwork, info) c c --- tilde r0 and tilde rl (small vectors) c rwork(1,y0) = -one call dcopy (l-1, rwork(2,z), 1, rwork(2,y0), 1) call dgetrs ('n', l-1, 1, rwork(2,zz+1), l+1, iwork, $ rwork(2,y0), l+1, info) rwork(l+1,y0) = zero rwork(1,yl) = zero call dcopy (l-1, rwork(2,z+l), 1, rwork(2,yl), 1) call dgetrs ('n', l-1, 1, rwork(2,zz+1), l+1, iwork, $ rwork(2,yl), l+1, info) rwork(l+1,yl) = -one c c --- Convex combination c call dsymv ('u', l+1, one, rwork(1,z), l+1, $ rwork(1,y0), 1, zero, rwork(1,y), 1) kappa0 = sqrt(ddot (l+1, rwork(1,y0), 1, $ rwork(1,y), 1)) call dsymv ('u', l+1, one, rwork(1,z), l+1, $ rwork(1,yl), 1, zero, rwork(1,y), 1) kappal = sqrt(ddot (l+1, rwork(1,yl), 1, $ rwork(1,y), 1)) call dsymv ('u', l+1, one, rwork(1,z), l+1, $ rwork(1,y0), 1, zero, rwork(1,y), 1) varrho = $ ddot (l+1, rwork(1,yl), 1, rwork(1,y), 1) $ / (kappa0*kappal) hatgamma = $ sign(1d0,varrho)*max(abs(varrho),7d-1) $ * (kappa0/kappal) call daxpy (l+1, (-hatgamma), rwork(1,yl), 1, $ rwork(1,y0), 1) c c --- Update c omega = rwork(l+1,y0) call dgemv ('n', n, l, (-one), work(1,u+1), n, $ rwork(2,y0), 1, one, work(1,u), 1) call dgemv ('n', n, l, one, work(1,r), n, $ rwork(2,y0), 1, one, x, 1) call dgemv ('n', n, l, (-one), work(1,r+1), n, $ rwork(2,y0), 1, one, work(1,r), 1) call dsymv ('u', l+1, one, rwork(1,z), l+1, $ rwork(1,y0), 1, zero, rwork(1,y), 1) rnrm = sqrt (ddot (l+1, rwork(1,y0), 1, $ rwork(1,y), 1)) c c ================================ c --- The reliable update part --- c ================================ 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 mv (n, x, work(1,r)) call solve (n, work(1,r)) do i=1,n work(i,r) = work(i,bp) - work(i,r) enddo mxnrmr = rnrm if (xpdt) then call daxpy (n, one, x, 1, work(1,xp), 1) call dlaset ('n', n, 1, zero, zero, x, 1) call dcopy (n, work(1,r), 1, work(1,bp), 1) mxnrmx = rnrm endif endif enddo c c ========================= c --- End of iterations --- c ========================= c call daxpy (n, one, work(1,xp), 1, x, 1) c c --- Check stopping criterion c call mv (n, x, work(1,r)) do i=1,n work(i,r) = b(i) - work(i,r) enddo call solve (n, work(1,r)) rnrm = dnrm2 (n, work(1,r), 1) if (rnrm.gt.tol*rnrm0) info = 1 c c --- Return c tol = rnrm/rnrm0 mxmv = nmv return end