subroutine zbistbl (l, n, x, b, mv, solve, tol,
     $     mxmv, work, ldw, rwork, ldrw, iwork, info)
c
c subroutine zbistbl v1.1 1998    
c 
c Copyright (c) 1995-1998 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
      complex*16 x(n), b(n)
      real*8 tol
      complex*16 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
      complex*16 alpha, beta, omega, rho0, rho1, sigma
      complex*16 varrho, hatgamma
      real*8 rnrm0, rnrm
      real*8 mxnrmx, mxnrmr
      complex*16 kappa0, kappal
c
c     .. Work Aliases ..
c
      integer z, zz, y0, yl, y
      integer rr, r, u, xp, bp
c
c     .. Constants ..
c
      real*8 delta
      parameter (delta = 1d-2)
      complex*16 zzero, zone
      parameter (zzero = (0d0,0d0), zone = (1d0,0d0))
c
c     .. BLAS and LAPACK ..
c
c     subroutine zaxpy
c     subroutine zcopy
c     subroutine zgemv
c     subroutine zgetrf
c     subroutine zgetrs
c     subroutine zlacpy
c     subroutine zlaset
c     subroutine zsymv
c     subroutine zhemv
c     function zdotc
c     function dznrm2
c
      real*8 dznrm2
      complex*16 zdotc
c
c     .. Intrinsic ..
c
      intrinsic abs, max, 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.0d0) 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 zcopy (n, work(1,r), 1, work(1,rr), 1)
      call zcopy (n, work(1,r), 1, work(1,bp), 1)
      call zcopy (n, x, 1, work(1,xp), 1)
      call zlaset ('n', n, 1, zzero, zzero, x, 1)
      rnrm0 = dznrm2 (n, work(1,r), 1)
      rnrm = rnrm0

      mxnrmx = rnrm0
      mxnrmr = rnrm0
      rcmp = .false.
      xpdt = .false.

      alpha = zzero
      omega = zone
      sigma = zone
      rho0 = zone
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 = zdotc (n, work(1,rr), 1, work(1,r+k-1), 1)
            if (rho0.eq.zzero) 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 = zdotc (n, work(1,rr), 1, work(1,u+k), 1)
            if (sigma.eq.zzero) then
               info = 2
               return
            endif
            alpha = rho1/sigma
            call zaxpy (n, alpha, work(1,u), 1, x, 1)
            do j=0,k-1
               call zaxpy (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 = dznrm2 (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 zgemv ('c', n, l+1-(i-1), zone, work(1,r+i-1),
     $           n, work(1,r+i-1), 1, zzero, rwork(i,z+i-1), 1)
            if (l-(i-1).ne.0) then
               call zcopy (l-(i-1), rwork(i+1,z+i-1), 1,
     $              rwork(i,z+i), l+1)
               call zlacgv (l-(i-1), rwork(i,z+i), l+1)
            endif
         enddo
         call zlacpy ('a', l+1, l+1, rwork(1,z), l+1,
     $        rwork(1,zz), l+1)
         call zgetrf (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) = -zone
         call zcopy (l-1, rwork(2,z), 1, rwork(2,y0), 1)
         call zgetrs ('n', l-1, 1, rwork(2,zz+1), l+1, iwork,
     $        rwork(2,y0), l+1, info)
         rwork(l+1,y0) = zzero

         rwork(1,yl) = zzero
         call zcopy (l-1, rwork(2,z+l), 1, rwork(2,yl), 1)
         call zgetrs ('n', l-1, 1, rwork(2,zz+1), l+1, iwork,
     $        rwork(2,yl), l+1, info)
         rwork(l+1,yl) = -zone
c
c        --- Convex combination
c
         call zhemv ('u', l+1, zone, rwork(1,z), l+1,
     $        rwork(1,y0), 1, zzero, rwork(1,y), 1)
         kappa0 = sqrt(zdotc (l+1, rwork(1,y0), 1,
     $        rwork(1,y), 1))

         call zhemv ('u', l+1, zone, rwork(1,z), l+1,
     $        rwork(1,yl), 1, zzero, rwork(1,y), 1)
         kappal = sqrt(zdotc (l+1, rwork(1,yl), 1,
     $        rwork(1,y), 1))

         call zhemv ('u', l+1, zone, rwork(1,z), l+1,
     $        rwork(1,y0), 1, zzero, rwork(1,y), 1)
         varrho =
     $        zdotc (l+1, rwork(1,yl), 1, rwork(1,y), 1)
     $        / (kappa0*kappal)

         hatgamma = 
     $        varrho/abs(varrho)*max(abs(varrho),7d-1)
     $        * (kappa0/kappal)

         call zaxpy (l+1, (-hatgamma), rwork(1,yl), 1,
     $        rwork(1,y0), 1)
c
c        --- Update
c
         omega = rwork(l+1,y0)

         call zgemv ('n', n, l, (-zone), work(1,u+1), n,
     $        rwork(2,y0), 1, zone, work(1,u), 1)
         call zgemv ('n', n, l, zone, work(1,r), n,
     $        rwork(2,y0), 1, zone, x, 1)
         call zgemv ('n', n, l, (-zone), work(1,r+1), n,
     $        rwork(2,y0), 1, zone, work(1,r), 1)

         call zhemv ('u', l+1, zone, rwork(1,z), l+1,
     $        rwork(1,y0), 1, zzero, rwork(1,y), 1)
         rnrm = sqrt (zdotc (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 zaxpy (n, zone, x, 1, work(1,xp), 1)
               call zlaset ('n', n, 1, zzero, zzero, x, 1)
               call zcopy (n, work(1,r), 1, work(1,bp), 1)
               mxnrmx = rnrm
            endif
         endif
      enddo
c
c     =========================
c     --- End of iterations ---
c     =========================
c
      call zaxpy (n, zone, 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 = dznrm2 (n, work(1,r), 1)
      if (rnrm.gt.tol*rnrm0) info = 1
c
c     --- Return
c
      tol = rnrm/rnrm0
      mxmv = nmv

      return 
      end