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
C Parameter list (by M.Botchev, October, 1997):
C
C l      == (input) INTEGER BiCGstab's dimension
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 b      == (input) DOUBLE PRECISION array dimension n
C           right-hand side (rhs) vector
C mv     == (input) EXTERNAL name of matrix vector subroutine
C           to deliver y:=A*x by CALL mv(n,x,y)
C solve  == (input) EXTERNAL name of subroutine to perform
C           a preconditioner solve by CALL solve(n,x). This call 
C           delivers solution to M*x=x, M is the left preconditioning 
C           matrix, x is rhs on input and solution on output.
C tol    == (input) DOUBLE PRECISION tolerance for the relative
C           residual norm stopping criterion: ||r_k|| <= tol*||r_0||
C mxmv   == (input/output) INTEGER.  On input: maximum number of matrix 
C           vector multiplications allowed to be done.  On input: 
C           actual number of matrix vector multiplications done
C work   == (workspace) DOUBLE PRECISION array dimension (n,3+2*(l+1))
C ldw    == (input) INTEGER leading dimension of work (??)
C rwork  == (workspace) DOUBLE PRECISION array dimension (l+1,3+2*(l+1))
C ldrw   == (input) INTEGER leading dimension of rwork (??)
C iwork  == (workspace) INTEGER array dimension (l+1)
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 to get rid of this)
C Three parameters were added to link the code to VAC:
C oktest == (input) LOGICAL.  Set to .true. if intermediate residuals should 
C           be printed
C nonzer == (input) LOGICAL, indicates whether the initial guess is zero
C           vector.  If (nonzer.eq..false.) one matrix-vector 
C           multiplication will be saved
C stc    == (input) CHARACTER*3, stopping criterion parameter:  
C           stc='abs' -- use absolute stopping criterion ||r_k|| <= tol,
C           otherwise use relative stopping criterion, as indicated above

      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

*         print *,'2: work(1,r) = ',work(1,r),' ',work(2,r),' ',work(3,r)
*         print *,'3:         x = ',x(1),' ',x(2),' ',x(3)
*         pause

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
*         do i=1,l+1
*            write(*,'(20f9.2)') (rwork(i,j),j=1,3+2*(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))
            nmv = nmv + 1
	    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
c         print *,'rnrm =',rnrm
      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