subroutine adjust (n,ndim,maxnzz,jcoef,key) 
c
c ... adjust makes adjustments to the jcoef array.
c
c ... parameters -- 
c
c          n      dimension of the matrix.
c          ndim   row dimension of jcoef array in defining routine
c          maxnz  number of columns in jcoef array
c          jcoef  integer matrix representation array
c          key    indicator flag
c                  = 1   remove zeros from jcoef array
c                  = 2   restore zeros to jcoef array
c
c ... specifications for parameters
c
      integer   jcoef(ndim,1) 
c
      maxnz = maxnzz
      if (maxnz .lt. 2) return
      if (key .eq. 2) go to 20
c
c ... change zero elements of jcoef array.
c
      do 15 j = 2,maxnz
         do 10 i = 1,n
 10      if (jcoef(i,j) .le. 0) jcoef(i,j) = i
 15   continue
      return
c
c ... put original zeros back in jcoef array.
c
 20   do 30 j = 2,maxnz
         do 25 i = 1,n
 25      if (jcoef(i,j) .eq. i) jcoef(i,j) = 0
 30   continue
      return
      end 
      subroutine adinfn (nn,ndim,maxnzz,jcoef,coef,nstore,ainf,wksp)
c
c ... adinfn computes an upper bound on the spectral radius of
c     inv(d)*a.
c
c ... parameters -- 
c
c         n       order of system (= nn)
c         ndim    row dimension of coef array in defining routine
c         maxnz   number of columns in coef array (= maxnzz)
c         jcoef   integer matrix representation array
c         coef    matrix representation array
c         nstore  matrix storage mode
c                  = 2   symmetric diagonal format
c                  = 3   nonsymmetric diagonal format
c         ainf    upper bound estimate upon output
c         wksp    workspace vector of length n
c
c ... specifications for parameters
c
      integer   jcoef(2)
      dimension coef(ndim,1), wksp(1)
c
      n = nn
      maxnz = maxnzz
      if (ainf .gt. 0.0) return
      do 10 i = 1,n 
 10   wksp(i) = coef(i,1)
      do 25 jd = 1,maxnz
         do 20 j = 1,maxnz
            if (jcoef(j) .ne. jd) go to 20
            do 15 i = 1,n
 15         wksp(i) = wksp(i) - abs (coef(i,j))
            if (nstore .eq. 3) go to 25 
            len = n - jd
            do 18 i = 1,len
 18         wksp(i+jd) = wksp(i+jd) - abs (coef(i,j))
            go to 25
 20      continue
         go to 30
 25   continue
 30   if (nstore .eq. 2) go to 50
      do 45 jd = 1,maxnz
         do 40 j = 1,maxnz
            if (jcoef(j) .ne. -jd) go to 40
            do 35 i = 1,n
 35         wksp(i) = wksp(i) - abs (coef(i,j))
            go to 45
 40      continue
         go to 50
 45   continue
c
c ... factor.
c
 50   t1 = vmin (n,wksp)
      if (t1 .le. 0.0) t1 = 1.0
      call ainfn (n,ndim,maxnz,jcoef,coef,nstore,ainf,wksp) 
      ainf = ainf/t1
      return
      end 
      subroutine ainfn (nn,ndim,maxnzz,jcoef,coef,nstore,ainf,
     a                  wksp) 
c
c ... ainfn calculates the infinity norm of the matrix a.
c
c ... parameters -- 
c
c         n       order of system (= nn)
c         ndim    row dimension of coef array in defining routine
c         maxnz   number of columns in coef array (= maxnzz)
c         jcoef   integer matrix representation array
c         coef    matrix representation array
c         nstore  matrix storage mode
c                  = 1   purdue format
c                  = 2   symmetric diagonal format
c                  = 3   nonsymmetric diagonal format
c                  = 4   symmetric sparse format
c                  = 5   nonsymmetric sparse format
c         ainf    the infinity norm of the matrix, //a//, upon
c                  output
c         wksp    workspace vector of length n
c
c ... specifications for parameters
c
      integer   jcoef(ndim,2) 
      dimension coef(ndim,1), wksp(1)
c
      n = nn
      maxnz = maxnzz
      if (ainf .gt. 0.0) return
      go to (10,30,55,75,75), nstore
c
c ... ellpack data structure. 
c
 10   do 15 i = 1,n 
 15   wksp(i) = abs (coef(i,1))
      if (maxnz .le. 1) go to 995
      do 25 j = 2,maxnz
         do 20 i = 1,n
 20      wksp(i) = wksp(i) + abs (coef(i,j))
 25   continue
      go to 995
c
c ... symmetric diagonal data structure.
c
 30   do 35 i = 1,n 
 35   wksp(i) = abs (coef(i,1))
      if (maxnz .le. 1) go to 995
      do 50 j = 2,maxnz
         ind = jcoef(j,1)
         len = n - ind
         do 40 i = 1,len
 40      wksp(i) = wksp(i) + abs (coef(i,j))
         do 45 i = 1,len
 45      wksp(i+ind) = wksp(i+ind) + abs (coef(i,j))
 50   continue
      go to 995
c
c ... nonsymmetric diagonal data structure.
c
 55   do 60 i = 1,n 
 60   wksp(i) = abs (coef(i,1))
      if (maxnz .le. 1) go to 995
      do 70 j = 2,maxnz
         ind = jcoef(j,1)
         len = n - iabs(ind)
         ist1 = max0(1,1 - ind)
         ist2 = min0(n,n - ind)
         do 65 i = ist1,ist2
 65      wksp(i) = wksp(i) + abs (coef(i,j))
 70   continue
      go to 995
c
c ... sparse structure.
c
 75   do 80 i = 1,n 
 80   wksp(i) = abs (coef(i,1))
      if (maxnz .le. n) go to 995
      np1 = n + 1
      do 85 k = np1,maxnz
 85   wksp(jcoef(k,1)) = wksp(jcoef(k,1)) + abs (coef(k,1)) 
      if (nstore .eq. 5) go to 995
      do 90 k = np1,maxnz
 90   wksp(jcoef(k,2)) = wksp(jcoef(k,2)) + abs (coef(k,1)) 
c
c ... determine ainf = max (wksp(i)).
c
 995  ainf = vmax (n,wksp)
      return
      end 
      subroutine bdfac (lda,nn,nsizee,nt,nb,a,isym)
c
c ... bdfac computes the factorization of a dense banded matrix.
c
c ... parameters -- 
c
c        lda    leading dimension of array a
c        n      active size of array a
c        nsize  size of an individual subsystem (if multiple systems) 
c                nsize = n upon input if not a multiple system
c        nt     number of diagonals needed to store the super-
c                diagonals
c        nb     number of diagonals needed to store the sub-
c                diagonals
c        a      array
c        isym   symmetry switch
c                = 0   matrix is symmetric
c                = 1   matrix is nonsymmetric
c
c ... specifications for parameters
c
      dimension a(lda,5)
      data lenv / 10 /
c
      n = nn
      maxt = nt
      nsize = nsizee
      nsys = n/nsize
c
c ... branch on symmetry.
c
      if (isym .eq. 1) go to 30
c
c ... symmetric case.
c
c ... diagonal case (maxt = 0).
c
      if (maxt .ne. 0) go to 15
      do 10 i = 1,n 
 10   a(i,1) = 1.0/a(i,1)
      return
c
c ... tridiagonal case (maxt = 1).
c
 15   if (maxt .ne. 1) go to 20
      if (nsys .le. lenv) call tfac (n,a,a(1,2))
      if (nsys .gt. lenv) call tfacm (n,nsize,a,a(1,2))
      return
c
c ... pentadiagonal case (maxt = 2).
c
 20   if (maxt .ne. 2) go to 25
      if (nsys .le. lenv) call pfac (n,a,a(1,2),a(1,3))
      if (nsys .gt. lenv) call pfacm (n,nsize,a,a(1,2),a(1,3))
      return
c
c ... banded case (maxt .gt. 2).
c
 25   if (nsys .le. lenv) call bfac (lda,n,maxt,a,a(1,2))
      if (nsys .gt. lenv) call bfacm (n,nsize,nsys,maxt,a,a(1,2))
      return
c
c ... nonsymmetric case.
c
 30   maxb = nb
c
c ... diagonal case (maxt = maxb = 0).
c
      if (maxt .ne. 0 .or. maxb .ne. 0) go to 40
      do 35 i = 1,n 
 35   a(i,1) = 1.0/a(i,1)
      return
c
c ... tridiagonal case (maxt = maxb = 1).
c
 40   if (maxt .ne. 1 .or. maxb .ne. 1) go to 45
      if (nsys .le. lenv) call tfacn (n,a,a(1,2),a(2,3))
      if (nsys .gt. lenv) call tfacnm (n,nsize,a,a(1,2),a(2,3))
      return
c
c ... pentadiagonal case (maxt = maxb = 2).
c
 45   if (maxt .ne. 2 .or. maxb .ne. 2) go to 50
      if (nsys .le. lenv) call pfacn (n,a,a(1,2),a(1,3),a(2,4),
     a                                  a(3,5))
      if (nsys .gt. lenv) call pfacnm (n,nsize,a,a(1,2),a(1,3),
     a                                   a(2,4),a(3,5))
      return
c
c ... all other cases.
c
 50   if (nsys .le. lenv) call bfacn (lda,n,maxt,maxb,a,a(1,2),
     a                                a(1,maxt+2))
      if (nsys .gt. lenv) call bfacnm (n,nsize,nsys,maxt,maxb,a,a(1,2),
     a                                 a(1,maxt+2))
      return
      end 
      subroutine bdinv (lda,nn,nsizee,nt,nb,fac,isym)
c
c ... bdinv computes the inverse of a dense banded matrix.
c
c ... parameters -- 
c
c        lda    leading dimension of factorization matrix fac
c        n      active size of factorization matrix fac
c        nsize  size of an individual subsystem (if multiple systems) 
c                nsize = n upon input if not a multiple system
c        nt     number of diagonals needed to store the super-
c                diagonals
c        nb     number of diagonals needed to store the sub-
c                diagonals
c        fac    array containing factorization upon input
c        isym   symmetry switch
c                = 0   matrix is symmetric
c                = 1   matrix is nonsymmetric
c
c ... specifications for parameters
c
      dimension fac(lda,3)
      data lenv / 10 /
c
      n = nn
      maxt = nt
      nsize = nsizee
      nsys = n/nsize
c
c ... branch on symmetry.
c
      if (isym .eq. 1) go to 30
c
c ... symmetric case.
c
      if (maxt - 1) 10,20,25
c
c ... diagonal case (maxt = 0).
c
 10   return
c
c ... tridiagonal case (maxt = 1).
c
 20   if (nsys .le. lenv) call tinv (n,fac,fac(1,2))
      if (nsys .gt. lenv) call tinvm (n,nsize,fac,fac(1,2)) 
      return
c
c ... banded case (maxt .ge. 2).
c
 25   call binv (lda,n,maxt+1,fac)
      return
c
c ... nonsymmetric case.
c
 30   maxb = nb
c
c ... diagonal case (maxt = maxb = 0).
c
      if (maxt .ne. 0 .or. maxb .ne. 0) go to 40
      return
c
c ... tridiagonal case (maxt = maxb = 1).
c
 40   if (maxt .ne. 1 .or. maxb .ne. 1) go to 45
      if (nsys .le. lenv) call tinvn (n,fac,fac(1,2),fac(2,3))
      if (nsys .gt. lenv) call tinvnm (n,nsize,fac,fac(1,2),fac(2,3)) 
      return
c
c ... all other cases.
c
 45   call binvn (lda,n,maxt,maxb,fac,fac(1,2),fac(1,maxt+2))
      return
      end 
      subroutine bdsol (lda,nn,nsizee,nt,nb,fac,y,x,isym)
c
c ... bdsol computes the solution to a dense banded matrix. 
c     thus, bdsol finds the solution to   a*x = y,  where fac
c     contains the factorization of the a matrix. 
c
c ... parameters -- 
c
c        lda    leading dimension of array fac
c        n      active size of array fac
c        nsize  size of an individual subsystem (if multiple systems) 
c                nsize = n upon input if not a multiple system
c        nt     number of diagonals needed to store the super-
c                diagonals of the factorization
c        nb     number of diagonals needed to store the sub-
c                diagonals of the factorization
c        fac    array containing the factorization of the matrix
c        y      upon input, y conains the right hand side
c        x      upon output, x contains the solution to  a*x = y
c        isym   symmetry switch
c                = 0   matrix is symmetric
c                = 1   matrix is nonsymmetric
c
c ... specifications for parameters
c
      dimension fac(lda,5), x(1), y(1)
      data lenv / 10 /
c
      n = nn
      maxt = nt
      nsize = nsizee
      nsys = n/nsize
c
c ... branch on symmetry.
c
      if (isym .eq. 1) go to 30
c
c ... symmetric case.
c
c ... diagonal case (maxt = 0).
c
      if (maxt .ne. 0) go to 15
      do 10 i = 1,n 
 10   x(i) = fac(i,1)*y(i)
      return
c
c ... tridiagonal case (maxt = 1).
c
 15   if (maxt .ne. 1) go to 20
      if (nsys .le. lenv) call tsoln (n,fac,fac(1,2),fac(1,2),y,x)
      if (nsys .gt. lenv) call tsolnm (n,nsize,fac,fac(1,2),
     a                                 fac(1,2),y,x)
      return
c
c ... pentadiagonal case (maxt = 2).
c
 20   if (maxt .ne. 2) go to 25
      if (nsys .le. lenv) call psoln (n,fac,fac(1,2),fac(1,3),
     a                                fac(1,2),fac(1,3),y,x)
      if (nsys .gt. lenv) call psolnm (n,nsize,fac,fac(1,2),fac(1,3), 
     a                                 fac(1,2),fac(1,3),y,x)
      return
c
c ... banded case (maxt .ge. 3).
c
 25   if (nsys .le. lenv) call bsol (lda,n,maxt,fac,fac(1,2),y,x)
      if (nsys .gt. lenv) call bsolm (n,nsize,maxt,fac,fac(1,2),y,x)
      return
c
c ... nonsymmetric case.
c
 30   maxb = nb
c
c ... diagonal case (maxt = maxb = 0).
c
      if (maxt .ne. 0 .or. maxb .ne. 0) go to 40
      do 35 i = 1,n 
 35   x(i) = fac(i,1)*y(i)
      return
c
c ... tridiagonal case (maxt = maxb = 1).
c
 40   if (maxt .ne. 1 .or. maxb .ne. 1) go to 45
      if (nsys .le. lenv) call tsoln (n,fac,fac(1,2),fac(2,3),y,x)
      if (nsys .gt. lenv) call tsolnm (n,nsize,fac,fac(1,2),fac(2,3), 
     a                                 y,x)
      return
c
c ... pentadiagonal case (maxt = maxb = 2).
c
 45   if (maxt .ne. 2 .or. maxb .ne. 2) go to 50
      if (nsys .le. lenv) call psoln (n,fac,fac(1,2),fac(1,3),
     a                                fac(2,4),fac(3,5),y,x)
      if (nsys .gt. lenv) call psolnm (n,nsize,fac,fac(1,2),fac(1,3), 
     a                                 fac(2,4),fac(3,5),y,x)
      return
c
c ... all other cases.
c
 50   if (nsys .le. lenv) call bsoln (lda,n,maxt,maxb,fac,fac(1,2),
     a                                fac(1,maxt+2),y,x)
      if (nsys .gt. lenv) call bsolnm (n,nsize,maxt,maxb,fac,
     a                            fac(1,2),fac(1,maxt+2),y,x)
      return
      end 
      subroutine bdsolt (lda,nn,nsizee,nt,nb,fac,y,x)
c
c ... bdsolt computes the transpose solution to a nonsymmetric
c     dense banded matrix.
c     thus, bdsolt finds the solution to   (a**t)*x = y,  where fac
c     contains the factorization of the a matrix. 
c
c ... parameters -- 
c
c        lda    leading dimension of array fac
c        n      active size of array fac
c        nsize  size of an individual subsystem (if multiple systems) 
c                nsize = n upon input if not a multiple system
c        nt     number of diagonals needed to store the super-
c                diagonals of the factorization
c        nb     number of diagonals needed to store the sub-
c                diagonals of the factorization
c        fac    array containing the factorization of the matrix
c        y      upon input, y conains the right hand side
c        x      upon output, x contains the solution to  a*x = y
c
c ... specifications for parameters
c
      dimension fac(lda,5), x(1), y(1)
      data lenv / 10 /
c
      n = nn
      maxt = nt
      maxb = nb
      nsize = nsizee
      nsys = n/nsize
c
c ... nonsymmetric case.
c
c ... diagonal case (maxt = maxb = 0).
c
      if (maxt .ne. 0 .or. maxb .ne. 0) go to 15
      do 10 i = 1,n 
 10   x(i) = fac(i,1)*y(i)
      return
c
c ... tridiagonal case (maxt = maxb = 1).
c
 15   if (maxt .ne. 1 .or. maxb .ne. 1) go to 20
      if (nsys .le. lenv) call tsoln (n,fac,fac(2,3),fac(1,2),y,x)
      if (nsys .gt. lenv) call tsolnm (n,nsize,fac,fac(2,3),fac(1,2), 
     a                                 y,x)
      return
c
c ... pentadiagonal case (maxt = maxb = 2).
c
 20   if (maxt .ne. 2 .or. maxb .ne. 2) go to 25
      if (nsys .le. lenv) call psoln (n,fac,fac(2,4),fac(3,5),
     a                                fac(1,2),fac(1,3),y,x)
      if (nsys .gt. lenv) call psolnm (n,nsize,fac,fac(2,4),fac(3,5), 
     a                                 fac(1,2),fac(1,3),y,x)
      return
c
c ... all other cases.
c
 25   if (nsys .le. lenv) call bsolnt (lda,n,maxt,maxb,fac,fac(1,2),
     a                                fac(1,maxt+2),y,x)
      if (nsys .gt. lenv) call bsontm (n,nsize,maxt,maxb,fac,
     a                            fac(1,2),fac(1,maxt+2),y,x)
      return
      end 
      subroutine bbs (ndim,nn,maxt,t,x) 
c
c ... bbs does a banded back substitution  (i + t)*x = y.
c     t is a rectangular matrix of adjacent super-diagonals.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system
c          maxt   number of columns in t array
c          t      array of active size n by maxt giving the super-
c                  diagonals in the order 1,2,3,...
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(ndim,1), x(1)
c
      n = nn
      do 20 i = n-1,1,-1
         sum = x(i) 
         lim = min0 (maxt,n-i)
         do 15 j = 1,lim
            sum = sum - t(i,j)*x(i+j)
 15      continue
         x(i) = sum 
 20   continue
      return
      end 
      subroutine bbsm (nsize,nsys,maxt,t,x)
c
c ... bbsm does a back solve  (i + t)*x = y  where t is an array
c     containing superdiagonals in order 1,2,... .
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of a single subsystem
c          nsys   number of independent subsystems
c          maxt   number of columns in t array
c          t      array of active size n by maxt containing 
c                  the super-diagonal elements of the factorization
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(nsize,nsys,1), x(nsize,1)
c
      do 25 i = nsize-1,1,-1
         lim = min0 (nsize-i, maxt)
         do 20 j = 1,lim
            ij = i + j
            do 15 l = 1,nsys
 15         x(i,l) = x(i,l) - t(i,l,j)*x(ij,l)
 20      continue
 25   continue
      return
      end 
      subroutine bbst (ndim,nn,maxb,b,x)
c
c ... bbst does a backward substitution  (i + (b**t))*x = y 
c     where the array b represents sub-diagonals.  b corresponds
c     to a banded system.
c
c ... parameters -- 
c
c          ndim   row dimension of b in defining routine
c          n      order of system (= nn)
c          maxb   number of diagonals stored in b 
c          b      array of active size n x maxb giving the
c                  sub-diagonals in the order -1,-2,... .
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension b(ndim,1), x(1)
c
      n = nn
      do 25 i = n,2,-1
         term = x(i)
         lim = min0 (i-1,maxb)
         do 20 j = 1,lim
            x(i-j) = x(i-j) - b(i,j)*term
 20      continue
 25   continue
      return
      end 
      subroutine bbstm (nsize,nsys,maxb,b,x)
c
c ... bbstm does the backward solve (i + (b**t))*x = y  where b
c     contains subdiagonals for multiple banded systems.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  the size of an individual subsystem
c          nsys   the number of subsystems
c          maxb   number of columns in b array
c          b      array of active size n by maxb containing 
c                  sub-diagonals in the order -1,-2,-3,...
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension b(nsize,nsys,1), x(nsize,1)
c
      do 25 i = nsize,2,-1
         lim = min0 (i-1,maxb)
         do 20 j = 1,lim
            do 15 l = 1,nsys
 15         x(i-j,l) = x(i-j,l) - b(i,l,j)*x(i,l) 
 20      continue
 25   continue
      return
      end 
      subroutine bfac (ndim,nn,maxt,d,t)
c
c ... bfac computes a factorization to a single banded
c     symmetric matrix represented by d and t and replaces it.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system (= nn)
c          maxt   number of columns in t array
c          d      vector containing the diagonal elements of a
c          t      array of active size n by maxt containing the
c                  super-diagonals in the order 1,2,3,...
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1)
c
      n = nn
      nm1 = n - 1
      do 20 k = 1,nm1
         pivot = d(k)
         lim = min0 (n-k,maxt)
         do 15 j1 = 1,lim
            term = t(k,j1)/pivot
            jcol1 = k + j1
            d(jcol1) = d(jcol1) - term*t(k,j1)
            if (j1 .eq. lim) go to 15
            j1p1 = j1 + 1
            do 10 j2 = j1p1,lim
               jcol2 = j2 - j1
               t(jcol1,jcol2) = t(jcol1,jcol2) - term*t(k,j2)
 10         continue
 15      continue
 20   continue
      do 25 i = 1,n 
 25   d(i) = 1.0/d(i)
      do 35 j = 1,maxt
         len = n - j
         do 30 i = 1,len
 30      t(i,j) = d(i)*t(i,j) 
 35   continue
      return
      end 
      subroutine bfacm (n,nsize,nsys,maxt,d,t)
c
c ... bfacm computes factorizations to multiple banded
c     symmetric matrices represented by d and t and replaces it.
c
c ... parameters -- 
c
c          n      order of global system (= nn)
c          nsize  order of a single system
c          nsys   number of independent subsystems
c          maxt   number of columns in t array
c          d      vector of length n containing the diagonal
c                  elements of a
c          t      array of active size n by maxt containing the
c                  super-diagonals in the order 1,2,3,...
c
c ... specifications for parameters
c
      dimension d(nsize,1), t(nsize,nsys,1)
c
      nsm1 = nsize - 1
      do 30 k = 1,nsm1
         lim = min0 (nsize-k,maxt)
         do 25 j1 = 1,lim
            jcol1 = k + j1
            do 10 l = 1,nsys
 10         d(jcol1,l) = d(jcol1,l) - (t(k,l,j1)**2)/d(k,l) 
            if (j1 .eq. lim) go to 25
            j1p1 = j1 + 1
            do 20 j2 = j1p1,lim
               jcol2 = j2 - j1
               do 15 l = 1,nsys
                  t(jcol1,l,jcol2) = t(jcol1,l,jcol2)
     a               - t(k,l,j1)*t(k,l,j2)/d(k,l) 
 15            continue
 20         continue
 25      continue
 30   continue
      call vinv (n,d)
      do 35 jj = 1,maxt
         len = n - jj
         call vexopy (len,t(1,1,jj),d,t(1,1,jj),3)
 35   continue
      return
      end 
      subroutine bfacn (ndim,nn,maxt,maxb,d,t,b)
c
c ... bfacn computes a factorization to a single banded
c     nonsymmetric matrix represented by d, t, and b and
c     replaces it.
c
c ... parameters -- 
c
c          ndim   row dimension of t and b in defining routine
c          n      order of system (= nn)
c          maxt   number of diagonals stored in t 
c          maxb   number of diagonals stored in b 
c          d      vector of length n containing the diagonal
c                  elements of a
c          t      array of active size n x maxt giving the
c                  super-diagonals in the order 1,2,3,...
c          b      array of active size n x maxb giving the
c                  sub-diagonals in the order -1,-2,-3,...
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), b(ndim,1)
c
      n = nn
      nm1 = n - 1
      do 35 k = 1,nm1
         pivot = d(k)
         liml = min0 (maxb,n-k)
         limu = min0 (maxt,n-k)
         do 30 ip = 1,liml
            i = k + ip
            term = b(i,ip)/pivot
            do 25 jp = 1,limu 
               term1 = term*t(k,jp)
               l = jp - ip
               if (l) 10,15,20
 10            b(i,-l) = b(i,-l) - term1
               go to 25
 15            d(i) = d(i) - term1
               go to 25
 20            t(i,l) = t(i,l) - term1
 25         continue
 30      continue
 35   continue
c
      do 40 i = 1,n 
 40   d(i) = 1.0/d(i)
      do 50 j = 1,maxt
         len = n - j
         do 45 i = 1,len
 45      t(i,j) = d(i)*t(i,j) 
 50   continue
      do 60 j = 1,maxb
         len = n - j
         do 55 i = 1,len
 55      b(i+j,j) = d(i)*b(i+j,j)
 60   continue
      return
      end 
      subroutine bfacnm (nn,nsize,nsys,maxt,maxb,d,t,b)
c
c ... bfacnm computes a factorization to multiple banded
c     nonsymmetric matrices represented by d, t, and b and
c     replaces it.
c
c ... parameters -- 
c
c          nsize  size of a subsystem
c          nsys   number of independent subsystems
c          maxt   number of diagonals stored in t 
c          maxb   number of diagonals stored in b 
c          n      order of system (= nn)
c          d      vector of length n containing the diagonal
c                  elements of a
c          t      array of active size n x maxt giving the
c                  super-diagonals in the order 1,2,3,...
c          b      array of active size n x maxb giving the
c                  sub-diagonals in the order -1,-2,-3,...
c
c ... specifications for parameters
c
      dimension d(nsize,1), t(nsize,nsys,1), b(nsize,nsys,1)
c
      n = nn
      nsm1 = nsize - 1
      do 50 k = 1,nsm1
         liml = min0 (maxb,nsize-k)
         limu = min0 (maxt,nsize-k)
         do 45 ip = 1,liml
            i = k + ip
            do 40 jp = 1,limu 
               l = jp - ip
               if (l) 10,20,30
 10            do 15 m = 1,nsys
 15            b(i,m,-l) = b(i,m,-l) - b(i,m,ip)*t(k,m,jp)/d(k,m)
               go to 40
 20            do 25 m = 1,nsys
 25            d(i,m) = d(i,m) - b(i,m,ip)*t(k,m,jp)/d(k,m) 
               go to 40
 30            do 35 m = 1,nsys
 35            t(i,m,l) = t(i,m,l) - b(i,m,ip)*t(k,m,jp)/d(k,m)
 40         continue
 45      continue
 50   continue
c
      call vinv (n,d)
      do 55 j = 1,maxt
         len = n - j
         call vexopy (len,t(1,1,j),d,t(1,1,j),3)
 55   continue
      do 60 j = 1,maxb
         len = n - j
         call vexopy (len,b(j+1,1,j),d,b(j+1,1,j),3)
 60   continue
      return
      end 
      subroutine bfs (ndim,nn,maxb,b,x) 
c
c ... bfs does a forward substitution  (i + b)*x = y  where the
c     array b represents sub-diagonals.  b corresponds to a 
c     banded system.
c
c ... parameters -- 
c
c          ndim   row dimension of b in defining routine
c          n      order of system (= nn)
c          maxb   number of diagonals stored in b 
c          b      array of active size n x maxb giving the
c                  sub-diagonals in the order -1,-2,-3,... .
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension b(ndim,1), x(1)
c
      n = nn
      do 15 i = 2,n 
         lim = min0 (i-1,maxb)
         sum = x(i) 
         do 10 j = 1,lim
            sum = sum - b(i,j)*x(i-j)
 10      continue
         x(i) = sum 
 15   continue
      return
      end 
      subroutine bfsm (nsize,nsys,maxb,b,x)
c
c ... bfsm does the forward solve (i + b)*x = y  where b contains
c     subdiagonals for multiple banded systems.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  the size of an individual subsystem
c          nsys   the number of subsystems
c          maxb   number of columns in b array
c          b      array of active size n by maxb containing 
c                  sub-diagonals in the order -1,-2,-3,... .
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension b(nsize,nsys,1), x(nsize,1)
c
      do 20 i = 2,nsize
         lim = min0 (i-1,maxb)
         do 15 j = 1,lim
            do 10 l = 1,nsys
 10         x(i,l) = x(i,l) - b(i,l,j)*x(i-j,l)
 15      continue
 20   continue
      return
      end 
      subroutine bfst (ndim,nn,maxt,t,x)
c
c ... bfst does a banded forward substitution  (i + (t**t))*x = y.
c     t is a rectangular matrix of adjacent super-diagonals.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system
c          maxt   number of columns in t array
c          t      array of active size n by maxt giving the super-
c                  diagonals in the order 1,2,3,...
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(ndim,1), x(1)
c
      n = nn
      nm1 = n - 1
      do 20 i = 1,nm1
         term = x(i)
         lim = min0 (maxt,n-i)
         do 15 j = 1,lim
            x(i+j) = x(i+j) - t(i,j)*term
 15      continue
 20   continue
      return
      end 
      subroutine bfstm (nsize,nsys,maxt,t,x)
c
c ... bfstm does a forward solve  (i + (t**t))*x = y  where t is
c     an array containing superdiagonals in order 1,2,... . 
c     (multiple systems)
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of a single subsystem
c          nsys   number of independent subsystems
c          maxt   number of columns in t array
c          t      array of active size n by maxt containing 
c                  the super-diagonal elements of the factorization
c          x      on input, x contains y
c                 vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(nsize,nsys,1), x(nsize,1)
c
      nsm1 = nsize - 1
      do 20 i = 1,nsm1
         lim = min0 (maxt,nsize-i)
         do 15 j = 1,lim
            ij = i + j
            do 10 l = 1,nsys
 10         x(ij,l) = x(ij,l) - t(i,l,j)*x(i,l)
 15      continue
 20   continue
      return
      end 
      subroutine binv (ndim,nn,maxnz,fact)
c
c ... binv computes an approximate inverse to a single banded
c     symmetric matrix.  fact must contain upon input the output
c     from a factorization routine.
c
c ... parameters -- 
c
c          ndim   row dimension of fact in the defining routine
c          n      order of system (= nn)
c          maxnz  bandwidth of the factorization and inverse
c          fact   array containing factorization diagonals
c                  in the order 0,1,2,3,...
c
c ... specifications for parameters
c
      dimension fact(ndim,2)
c
      n = nn
      nm1 = n - 1
c
c ... general banded matrix.
c
      do 25 ik = 1,nm1
         k = n - ik 
         lim = min0 (ik+1,maxnz)
         sum1= 0.0
         do 15 i = 2,lim
            t1 = fact(k,i)
            sum2= 0.0
            do 10 j = 2,lim
               m1 = min0 (i,j)
               m2 = max0 (i,j)
               l1 = k + m1 - 1
               l2 = m2 - m1 + 1
               sum2 = sum2 - fact(k,j)*fact(l1,l2)
 10         continue
            fact(n,i) = sum2
            sum1 = sum1 - t1*sum2
 15      continue
         fact(k,1) = fact(k,1) + sum1
         do 20 i = 2,lim
 20      fact(k,i) = fact(n,i)
 25   continue
      do 30 i = 2,maxnz
 30   fact(n,i)= 0.0
      return
      end 
      subroutine binvn (ndim,nn,maxt,maxb,d,t,b)
c
c ... binvn computes an approximate inverse to a single banded
c     nonsymmetric matrix.  d, t, and b must contain upon input
c     the output from a factorization routine.
c
c ... parameters -- 
c
c          ndim   row dimension of t and b in the defining routine
c          n      order of system (= nn)
c          maxt   number of columns in t
c          maxb   number of columns in b
c          d      vector of length n containing the diagonal
c                  elements of the factorization
c          t      array of active size n by maxt containing 
c                  the superdiagonals of the factorization
c                  in the order 1,2,3,...
c          b      array of active size n by maxb containing 
c                  the subdiagonals of the factorization
c                   in the order -1,-2,-3,....
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), b(ndim,1)
c
      n = nn
      nm1 = n - 1
c
c ... general banded matrix.
c
      do 75 ik = 1,nm1
         k = n - ik 
c
c ... copy kth row and column into wksp.
c
         limr = min0 (maxt,ik)
         limc = min0 (maxb,ik)
         do 10 j = 1,limr
 10      t(n,j) = t(k,j)
         do 15 j = 1,limc
 15      b(1,j) = b(k+j,j)
c
c ... do computations for kth row.
c
         do 40 j = 1,limr
            sum= 0.0
            lim = min0 (limr,limc+j)
            do 35 i = 1,lim
               kpi = k + i
               l = i - j
               if (l) 20,25,30
 20            sum = sum - t(n,i)*t(kpi,-l)
               go to 35
 25            sum = sum - t(n,i)*d(kpi)
               go to 35
 30            sum = sum - t(n,i)*b(kpi,l)
 35         continue
            t(k,j) = sum
 40      continue
c
c ... do computations for kth column.
c
         do 65 j = 1,limc
            sum= 0.0
            lim = min0 (limc,limr+j)
            kpj = k + j
            do 60 i = 1,lim
               kpi = k + i
               l = i - j
               if (l) 45,50,55
 45            sum = sum - b(1,i)*b(kpj,-l)
               go to 60
 50            sum = sum - b(1,i)*d(kpi)
               go to 60
 55            sum = sum - b(1,i)*t(kpj,l)
 60         continue
            b(kpj,j) = sum
 65      continue
c
c ... compute kth diagonal element.
c
         sum = d(k) 
         lim = min0 (limr,limc)
         do 70 j = 1,lim
 70      sum = sum - t(n,j)*b(k+j,j)
         d(k) = sum 
 75   continue
c
c ... zero out workspace rows.
c
      do 80 j = 1,maxt
 80   t(n,j)= 0.0
      do 85 j = 1,maxb
 85   b(1,j)= 0.0
      return
      end 
      subroutine bmul (ndim,n,maxt,d,t,x,y)
c
c ... bmul computes y = a*x, where x and y are vectors and
c ... a is a banded symmetric matrix.
c
c ... parameters -- 
c
c         ndim          row dimension of array t
c         n             order of matrix 
c         maxt          number of columns in t
c         d             vector of length n giving the
c                        diagonal elements of a
c         t             array of size n by maxt giving the
c                        superdiagonals of a in the order
c                        1,2,....
c         x,y           vectors of order n
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), x(1), y(1)
c
      do 10 i = 1,n 
 10   y(i) = d(i)*x(i)
      if (maxt .le. 0) return 
      do 25 la = 1,maxt
         len = n - la
         do 15 i = 1,len
 15      y(i) = y(i) + t(i,la)*x(i+la)
         do 20 i = 1,len
 20      y(i+la) = y(i+la) + t(i,la)*x(i)
 25   continue
      return
      end 
      subroutine bmuln (ndim,n,maxt,maxb,d,t,b,x,y)
c
c ... bmuln computes y = a*x, where x and y are vectors and 
c ... d, t, and b represent a stored in nonsymmetric band
c ... storage format.
c
c ... parameters -- 
c
c         ndim          row dimension of arrays t and b
c         n             order of array a
c         maxt          number of columns in t array
c         maxb          number of columns in b array
c         d             vector of length n giving the diagonal
c                        elements of a
c         t             array of active size n by maxt giving
c                        the super-diagonals of a in the order
c                        1,2,3,...
c         b             array of active size n by maxb giving
c                        the sub-diagonals of a in the order
c                        -1,-2,-3,....
c         x,y           vectors of order n
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), b(ndim,1), x(1), y(1)
c
      do 10 i = 1,n 
 10   y(i) = d(i)*x(i)
      if (maxt .lt. 1) go to 25
      do 20 j = 1,maxt
         len = n - j
         do 15 i = 1,len
 15      y(i) = y(i) + t(i,j)*x(i+j)
 20   continue
 25   if (maxb .lt. 1) return 
      do 35 j = 1,maxb
         len = n - j
         do 30 i = 1,len
 30      y(i+j) = y(i+j) + b(i+j,j)*x(i)
 35   continue
      return
      end 
      subroutine bmulnt (ndim,n,maxt,maxb,d,t,b,x,y)
c
c ... bmulnt computes y = (a**t)*x, where x and y are vectors and
c ... d, t, and b represent a stored in nonsymmetric band
c ... storage format.
c
c ... parameters -- 
c
c         ndim          row dimension of arrays t and b
c         n             order of array a
c         maxt          number of columns in t array
c         maxb          number of columns in b array
c         d             vector of length n giving the diagonal
c                        elements of a
c         t             array of active size n by maxt giving
c                        the super-diagonals of a in the order
c                        1,2,3,...
c         b             array of active size n by maxb giving
c                        the sub-diagonals of a in the order
c                        -1,-2,-3,...
c         x,y           vectors of order n
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), b(ndim,1), x(1), y(1)
c
      do 10 i = 1,n 
 10   y(i) = d(i)*x(i)
      if (maxt .lt. 1) go to 25
      do 20 j = 1,maxt
         len = n - j
         do 15 i = 1,len
 15      y(i+j) = y(i+j) + t(i,j)*x(i)
 20   continue
 25   if (maxb .lt. 1) return 
      do 35 j = 1,maxb
         len = n - j
         do 30 i = 1,len
 30      y(i) = y(i) + b(i+j,j)*x(i+j)
 35   continue
      return
      end 
      subroutine bsol (ndim,nn,maxt,d,t,y,x)
c
c ... bsol solves a*x = y for a banded and symmetric matrix a. d and
c     t must contain upon input the factorization arrays from bfac.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system
c          maxt   number of columns in t array
c          d      vector of length n containing the diagonal
c                  pivots of the factorization
c          t      array of active size n by maxt giving the super-
c                  diagonals of the factorization in the order
c                  1,2,3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(ndim,1), x(1), y(1), d(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      call bfst (ndim,n,maxt,t,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbs (ndim,n,maxt,t,x)
      return
      end 
      subroutine bsolm (nn,nsize,maxt,d,t,y,x)
c
c ... bsolm solves the system  ax = y  for x, where a is multiple
c     symmetric banded matrices whose factorizations are contained in 
c     d and t.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of a single subsystem
c          maxt   number of columns in t array
c          d      vector of length n containing the diagonal
c                  elements of the factorization
c          t      array of active size n by maxt containing 
c                  the super-diagonal elements of the factorization
c                  in the order 1,2,3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension d(1), t(1), y(1), x(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      nsys = n/nsize
      call bfstm (nsize,nsys,maxt,t,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbsm (nsize,nsys,maxt,t,x)
      return
      end 
      subroutine bsoln (ndim,nn,maxt,maxb,d,t,b,y,x)
c
c ... bsoln solves a*x = y for a banded and nonsymmetric matrix a.
c     d, t, and b must contain upon input the factorization arrays
c     from bfacn.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system
c          maxt   number of columns in t array
c          maxb   number of columns in b array
c          d      vector of length n containing the diagonal
c                  pivots of the factorization
c          t      array of active size n by maxt giving the super-
c                  diagonals of the factorization in the order
c                  1,2,3,...
c          b      array of active size n by maxb giving the sub-
c                  diagonals of the factorization in the order
c                  -1,-2,-3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(ndim,1), x(1), y(1), d(1), b(ndim,1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      call bfs (ndim,n,maxb,b,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbs (ndim,n,maxt,t,x)
      return
      end 
      subroutine bsolnm (nn,nsize,maxt,maxb,d,t,b,y,x)
c
c ... bsolnm solves a*x = y for a banded and nonsymmetric matrix a.
c     d, t, and b must contain upon input the factorization arrays
c     from bfacnm.  (multiple systems)
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of an individual subsystem 
c          maxt   number of columns in t array
c          maxb   number of columns in b array
c          d      vector of length n containing the diagonal
c                  pivots of the factorization
c          t      array of active size n by maxt giving the super-
c                  diagonals of the factorization in the order
c                  1,2,3,...
c          b      array of active size n by maxb giving the sub-
c                  diagonals of the factorization in the order
c                  -1,-2,-3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(1), x(1), y(1), d(1), b(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      nsys = n/nsize
      call bfsm (nsize,nsys,maxb,b,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbsm (nsize,nsys,maxt,t,x)
      return
      end 
      subroutine bsolnt (ndim,nn,maxt,maxb,d,t,b,y,x)
c
c ... bsolnt solves (a**t)*x = y for a banded and nonsymmetric
c     matrix a.  d, t, and b must contain upon input the
c     factorization arrays from bfacn.
c
c ... parameters -- 
c
c          ndim   row dimension of t array in defining routine
c          n      order of system
c          maxt   number of columns in t array
c          maxb   number of columns in b array
c          d      vector of length n containing the diagonal
c                  pivots of the factorization
c          t      array of active size n by maxt giving the super-
c                  diagonals of the factorization in the order
c                  1,2,3,...
c          b      array of active size n by maxb giving the sub-
c                  diagonals of the factorization in the order
c                  -1,-2,-3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(ndim,1), x(1), y(1), d(1), b(ndim,1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      call bfst (ndim,n,maxt,t,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbst (ndim,n,maxb,b,x)
      return
      end 
      subroutine bsontm (nn,nsize,maxt,maxb,d,t,b,y,x)
c
c ... bsontm solves (a**t)*x = y for a banded and nonsymmetric
c     matrix a.  d, t, and b must contain upon input the
c     factorization arrays from bfacnm.  (multiple systems) 
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of an individual subsystem 
c          maxt   number of columns in t array
c          maxb   number of columns in b array
c          d      vector of length n containing the diagonal
c                  pivots of the factorization
c          t      array of active size n by maxt giving the super-
c                  diagonals of the factorization in the order
c                  1,2,3,...
c          b      array of active size n by maxb giving the sub-
c                  diagonals of the factorization in the order
c                  -1,-2,-3,...
c          y      right-hand-side vector
c          x      vector containing solution upon output
c
c ... specifications for parameters
c
      dimension t(1), x(1), y(1), d(1), b(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      nsys = n/nsize
      call bfstm (nsize,nsys,maxt,t,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call bbstm (nsize,nsys,maxb,b,x)
      return
      end 
      subroutine bicol (n,nz,ia,ja,count,father,oppos,propa)
c
c ... bicolor determines whether or not the matrix represented
c     in the sparse (ia,ja) format is bi-colorable.
c     the algorithm used is the union-find algorithm.
c
c ... parameters -- 
c
c        n      number of vertices
c        nz     number of edges (length of ia and ja vectors)
c        ia     integer vector of i values
c        ja     integer vector of j values
c        count  integer workspace vectors of length n each
c        father upon output, count gives the color of each node
c        oppos
c        propa  logical variable indicating on output whether
c                matrix has property a
c
c ... specification of parameters
c
      logical propa 
      integer ia(1), ja(1), count(1), father(1), oppos(1)
      integer v, w, w0, a, b, c, d
c
      do 10 i = 1,n 
         count(i) = 1
         father(i) = 0
         oppos(i) = 0
 10   continue
      do 60 k = 1,nz
         if (ia(k) .eq. ja(k)) go to 60 
c
c ... a = find (ia(k)).
c
         v = ia(k)
 15      if (father(v) .eq. 0) go to 20 
         v = father(v)
         go to 15
 20      w = ia(k)
 25      if (father(w) .eq. 0) go to 30 
         w0 = w
         w = father(w)
         father(w0) = v
         go to 25
 30      a = v
c
c ... b = find (ja(k)).
c
         v = ja(k)
 35      if (father(v) .eq. 0) go to 40 
         v = father(v)
         go to 35
 40      w = ja(k)
 45      if (father(w) .eq. 0) go to 50 
         w0 = w
         w = father(w)
         father(w0) = v
         go to 45
 50      b = v
c
c ... test for a = b.
c
         if (a .ne. b) go to 55
         propa = .false.
         return
c
c ... do unioning.
c
 55      if (oppos(a) .eq. b) go to 60
         if (oppos(b) .eq. 0) then
            c = a
         else
c
c ... c = merge (a,oppos(b)). 
c
            i = a
            j = oppos(b)
            if (count(i) .ge. count(j)) then
               father(j) = i
               count(i) = count(i) + count(j)
               c = i
            else
               father(i) = j
               count(j) = count(i) + count(j)
               c = j
            endif
         endif
         if (oppos(a) .eq. 0) then
            d = b
         else
c
c ... d = merge (b,oppos(a)). 
c
            i = b
            j = oppos(a)
            if (count(i) .ge. count(j)) then
               father(j) = i
               count(i) = count(i) + count(j)
               d = i
            else
               father(i) = j
               count(j) = count(i) + count(j)
               d = j
            endif
         endif
         oppos(c) = d
         oppos(d) = c
 60   continue
c
c ... do coloring.
c
      do 65 i = 1,n 
 65   count(i) = 0
      do 90 i = 1,n 
c
c ... a = find(i).
c
         v = i
 70      if (father(v) .eq. 0) go to 75 
         v = father(v)
         go to 70
 75      w = i
 80      if (father(w) .eq. 0) go to 85 
         w0 = w
         w = father(w)
         father(w0) = v
         go to 80
 85      a = v
         if (count(a) .eq. 0) then
            count(a) = 1
            count(i) = 1
            j = oppos(a)
            if (j .ne. 0) count(j) = 2
         else
            count(i) = count(a)
         endif
 90   continue
      propa = .true.
      return
      end 
      subroutine chgcon (tri,ier)
c
c ... chgcon computes the new estimates for the largest and 
c     smallest eigenvalues (emax and emin) for conjugate gradient
c     acceleration. 
c
c ... parameters -- 
c
c          tri    tridiagonal matrix associated with the eigenvalues
c                    of the conjugate gradient polynomial
c          ier    error code
c
c ... specifications for parameters
c
      dimension tri(2,2)
c
c *** begin -- package common 
c
      common / itcom1 / in, itmax, level, nout, ns1, ns2, ns3,
     a      iplr, iqlr, ntest, is, iacel, idgts, nbl1d, nbl2d
      logical           halt, maxadp, minadp, maxadd, minadd
      common / itcom2 / halt, maxadp, minadp, maxadd, minadd
      common / itcom3 / alpha, beta, zeta, emax, emin, pap, 
     b                  alphao, gamma, sigma, rr, rho, dkq, dkm1,
     b                  ff, rqmin, rqmax, stptst, udnm, ubarnm,
     b                  bnorm, bnorm1
      common / itcom4 / keygs, srelpr, keyzer
c
c *** end   -- package common 
c
c     description of variables in common blocks in main routine
c
      save tl1,tl2,bl1,bl2
      ip = is
      if (ip - 1) 10,20,30
c
c ... ip = 0
c
 10   end = 1.0/alpha
      tri(1,1) = end
      tri(2,1)= 0.0 
      if (maxadp) emax = end
      if (minadp) emin = end
      return
c
c ... ip = 1
c
 20   t1 = 1.0/alpha + beta/alphao
      t2 = beta/(alphao**2)
      tri(1,2) = t1 
      tri(2,2) = t2 
      tsqr = sqrt (t2)
      tl1 = tri(1,1) + tsqr
      tl2 = t1 + tsqr
      bl1 = tri(1,1) - tsqr
      bl2 = t1 - tsqr
      t3 = tri(1,1) + t1
      t4 = sqrt ( (t1-tri(1,1))**2 + 4.0*t2 )
      if (maxadp) emax = (t3 + t4)/2.0
      if (minadp) emin = (t3 - t4)/2.0
      return
c
c ... ip .ge. 2
c
 30   t1 = 1.0/alpha + beta/alphao
      t2 = beta/(alphao**2)
      tsqr = sqrt (t2)
      tri(1,ip+1) = t1
      tri(2,ip+1) = t2
      if (.not. maxadp) go to 40
c
c ... compute new estimate of emax.
c
      tl1 = amax1 (tl1,tl2+tsqr)
      tl2 = t1 + tsqr
      emaxo = emax
      end = amax1 (tl1,tl2)
      e1 = eigvss (ip+1,tri,emaxo,end,2,ier)
      if (ier .ne. 3  .and.  ier .ne. 4) go to 35 
c
c ... poor estimate for emax.  therefore need to stop adaptive
c     procedure and keep old value of emax.
c
      maxadp = .false.
      if (level .ge. 2) write (nout,31) ier,in,emaxo
 31   format (/5x,'estimation of maximum eigenvalue emax halted'
     a        /5x,'routine zbrent returned ier = ',i5
     b        /5x,'adaptive procedure turned off at iteration ',i5
     c        /5x,'final estimate of maximum eigenvalue =',e15.7/)
      go to 40
c
c ... valid emax estimate.  check for small relative change in emax.
c
 35   emax = e1
      if (abs (emax - emaxo) .lt. emax*zeta) maxadp = .false.
c
c ... compute new estimate of emin.
c
 40   if (.not. minadp) return
      bl1 = amin1 (bl1,bl2-tsqr)
      bl2 = t1 - tsqr
      start = amax1 ( 0.0, amin1 (bl1,bl2) )
      emino = emin
      e1 = eigvss (ip+1,tri,start,emino,1,ier)
      if (ier .ne. 3  .and.  ier .ne. 4) go to 45 
c
c ... poor estimate for emin.  therefore need to stop adaptive
c     procedure and keep old value of emin.
c
      minadp = .false.
      if (level .ge. 2) write (nout,41) ier,in,emino
 41   format (/5x,'estimation of minimum eigenvalue emin halted'
     a        /5x,'routine zbrent returned ier = ',i5
     b        /5x,'adaptive procedure turned off at iteration ',i5
     c        /5x,'final estimate of minimum eigenvalue =',e15.7/)
      return
c
c ... valid emin estimate.  check for small relative change in emin.
c
 45   emin = e1
      if (abs (emin - emino) .lt. emin*zeta) minadp = .false.
      return
      end 
      subroutine chgsi (suba,coef,jcoef,wfac,jwfac,nn,z,wksp,
     a                  icode,ier)
c
c ... chgsi adapts on the iteration parameters.
c
c ... parameters -- 
c
c         n         order of system (= nn)
c         z         current pseudo-residual vector
c         wksp      workspace vector of length n
c         icode     output indicator of parameter changes
c                    = 0    estimates of emax, emin not changed
c                    = 1    estimates of emax, emin changed 
c         ier       error code
c
c ... specifications for parameters
c
      external suba 
      integer jcoef(2), jwfac(1)
      dimension z(1), wksp(1), coef(1), wfac(1)
c
c *** begin -- package common 
c
      common / itcom1 / in, itmax, level, nout, ns1, ns2, ns3,
     a      iplr, iqlr, ntest, is, iacel, idgts, nbl1d, nbl2d
      logical           halt, maxadp, minadp, maxadd, minadd
      common / itcom2 / halt, maxadp, minadp, maxadd, minadd
      common / itcom3 / alpha, beta, zeta, emax, emin, pap, 
     b                  alphao, gamma, sigma, rr, rho, dkq, dkm1,
     b                  ff, rqmin, rqmax, stptst, udnm, ubarnm,
     b                  bnorm, bnorm1
      common / itcom4 / keygs, srelpr, keyzer
      common / itcom9 / rhave, zhave, zthave, rcalp, zcalp, ztcalp,
     a           udhav, rdhav, rzhav, rzthav, zdhav, zzthav, ztdhav,
     a                  rdot, rzdot, rztdot, zdot, zztdot, ztdot
      logical rhave, zhave, zthave, rcalp, zcalp, ztcalp,
     a        udhav, rdhav, rzhav, rzthav, zdhav, zzthav, ztdhav
c
c *** end   -- package common 
c
      n = nn
c
      istar = 3
      icode = 0
      if (is .eq. 0) return
      rnrm = sqrt (rzdot)
      rnrmq = sqrt (dkq)
      rnrm1 = sqrt (dkm1)
      qa = rnrm/rnrmq
      t1 = rr**is
      qt = 2.0*sqrt (t1)/(1.0 + t1)
      if (qa .le. qt**ff) return
      if (qa .le. 1.0  .and.  is .le. istar) return
      icode = 1
c
c ... compute rayleigh quotient.
c ...        rq = (z,a*z)/(r,z)
c
      call suba (coef,jcoef,wfac,jwfac,n,z,wksp)
      top= 0.0
      do 10 i = 1,n 
 10   top = top + z(i)*wksp(i)
      if (top .ge. 0.0) go to 15
      ier = -6
      call ershow (ier,'chgsi')
      return
 15   rq = top/rzdot
      kode = 0
      if (rq .gt. rqmax) kode = 1
      rqmin = amin1 (rq,rqmin)
      rqmax = amax1 (rq,rqmax)
      yy = (1.0+t1)*(qa+sqrt (qa*qa-qt*qt))/2.0
      xx = yy**(1.0/float (is))
      if (qa .gt. 1.0) go to 25
      if (kode .eq. 1) go to 25
c
c ... emin adjustment.
c
      eminp = (emax+emin)*(1.0-xx)*(xx-rr)/(2.0*xx*(rr+1.0))
      if (minadp) emin = amin1 (emin,eminp,rqmin) 
      if (maxadp) emax = amax1 (emax,rqmax)
      if (level .ge. 2) write (nout,20) in,rq,eminp,emin,emax
 20   format (/1x,15x,'parameters were changed at iteration',i7/
     a        1x,20x,'rayleigh quotient  ',f15.9/ 
     a        1x,20x,'young estimate     ',f15.9/ 
     a        1x,20x,'emin               ',f15.9/ 
     a        1x,20x,'emax               ',f15.9/)
      return
c
c ... emax adjustment.
c
 25   emaxp = (emax+emin)*(1.0+xx)*(xx+rr)/(2.0*xx*(rr+1.0))
      uu = ((1.0+t1)/(1.0+rr**(is-1))) * (rnrm/rnrm1)
      emaxpp = (emax+emin)*(1.0+uu)*(uu+rr)/(2.0*uu*(rr+1.0))
      if (maxadp) emax = amax1 (emax,1.1*emaxp,1.1*emaxpp,1.1*rqmax)
      if (minadp) emin = rqmin
      if (level .ge. 2) write (nout,30) in,rq,emaxp,emaxpp,emin,emax
 30   format (/1x,15x,'parameters were changed at iteration',i7/
     a        1x,20x,'rayleigh quotient  ',f15.9/ 
     a        1x,20x,'young estimate     ',f15.9/ 
     a        1x,20x,'hageman estimate   ',f15.9/ 
     a        1x,20x,'emin               ',f15.9/ 
     a        1x,20x,'emax               ',f15.9/)
      return
      end 
      subroutine color (nxp,nyp,nzp,nx,ny,nz,pp,p)
c
c ... routine color reproduces a color pattern given by array
c     pp of dimensions nxp x nyp x nzp into the grid color
c     array p of dimensions nx x ny x nz.
c
c ... parameters -- 
c
c       nxp,    integer variables giving the x, y, and z dimensions
c        nyp,    of the pattern array, respectively.
c        nzp
c       nx,ny,  integer variables giving the x, y, and z dimensions
c        nz      of the grid, respectively.
c       pp      integer vector of length  nxp*nyp*nzp
c                giving the color pattern to be repeated
c       p       integer vector of length  nxg*nyg*nzg
c                which contains upon output the grid coloring
c
c ... specifications for parameters
c
      integer pp (nxp,nyp,nzp), p (nx,ny,nz)
c
      do 30 k = 1,nz
         kp = mod (k - 1, nzp) + 1
         do 20 j = 1,ny
            jp = mod (j - 1, nyp) + 1
            do 10 i = 1,nx
               ip = mod (i - 1, nxp) + 1
               p (i,j,k) = pp (ip,jp,kp)
 10         continue
 20      continue
 30   continue
      return
      end 
      subroutine defcon (ndim,nn,maxnz,jcoef,coef,kblsz,iblock,lbhb)
c
c ... define defines block constants for block-structured matrices.
c     (diagonal data structure, constant block size)
c
c ... parameters -- 
c
c         ndim     row dimension of coef array in defining routine
c         nn       size of system
c         maxnz    number of diagonals in coef
c         jcoef    integer vector of size maxnz giving the diagonal
c                   numbers
c         coef     matrix representation array
c         kblsz    constant block size
c         iblock   integer array of size 3 by lbhb
c                   giving block constants upon output
c         lbhb     integer giving the number of diagonal blocks
c                   upon output.
c
c ... specifications for parameters
c
      integer   jcoef(2), iblock(3,3)
      dimension coef(ndim,1)
c
      n = nn
      ipt = 2
      iblock(1,1) = 0
      iblock(1,2) = 0
      iblock(2,1) = 1
      iblock(3,1) = 0
      iblock(3,2) = 0
      do 25 j = 1,maxnz
         jd = jcoef(j)
         do 10 i = 1,n
            if (coef(i,j) .ne. 0.0) go to 15
 10      continue
         go to 25
 15      jcol = i + jd
c
c ... find block for jcol.
c
         ib = (i-1)/kblsz + 1 
         jb = (jcol-1)/kblsz + 1
         id = jb - ib
         if (id .eq. iblock(1,ipt)) go to 20
         ipt = ipt + 1
         iblock(1,ipt) = id
         iblock(3,ipt) = 0
 20      iblock(3,ipt) = iblock(3,ipt) + 1
 25   continue
      lbhb = ipt
c
c ... split zero diagonal block into super and sub diagonals.
c
      jlim = iblock(3,2)
      do 30 j = 1,jlim
         jd = jcoef(j)
         if (jd .lt. 0) go to 35
         iblock(3,1) = iblock(3,1) + 1
         iblock(3,2) = iblock(3,2) - 1
 30   continue
      j = jlim + 1
 35   iblock(2,2) = j
c
c ... form starting positions.
c
      if (lbhb .le. 2) return 
      iblock(2,3) = 1
      if (lbhb .le. 3) return 
      do 40 j = 4,lbhb
 40   iblock(2,j) = iblock(2,j-1) + iblock(3,j-1) 
      return
      end 
      subroutine define (ndim,maxnew,jcnew,coef,ncol,nc,
     a                   iblock,lbhb)
c
c ... define defines block constants for block-structured matrices.
c     (diagonal data structure, nonconstant block size)
c
c ... parameters -- 
c
c         ndim     row dimension of coef array in defining routine
c         maxnew   integer vector giving the number of diagonals
c                   for each distinct block size. 
c         jcnew    integer array of size ncolor*max(maxnew(i))
c                   giving the diagonal numbers for each distinct
c                   block size.
c         coef     matrix representation array
c         ncolor   number of distinct block sizes 
c         nc       integer vector of length ncolor, giving the number 
c                   of nodes for each distinct block size.
c         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants upon output
c         lbhb     integer vector of size ncolor giving the number
c                   of diagonal blocks for each distinct block size
c                   upon output.
c
c ... specifications for parameters
c
      integer   maxnew(ncol), jcnew(ncol,1), nc(ncol), lbhb(ncol),
     a          iblock(3,ncol,3)
      dimension coef(ndim,1)
c
      ncolor = ncol 
      ist = 1
      do 60 k = 1,ncolor
         ncc = nc(k)
         maxnz = maxnew(k)
         ied = ist + ncc - 1
         ipt = 2
         iblock(1,k,1) = 0
         iblock(1,k,2) = 0
         iblock(2,k,1) = 1
         iblock(3,k,1) = 0
         iblock(3,k,2) = 0
         do 35 j = 1,maxnz
            jd = jcnew(k,j)
            do 10 i = ist,ied 
               if (coef(i,j) .ne. 0.0) go to 15
 10         continue
            go to 35
 15         jcol = i + jd
c
c ... find block for jcol.
c
            ib = k
            js = 0
            do 20 ij = 1,ncolor
               js = js + nc(ij)
               if (js .ge. jcol) go to 25
 20         continue
 25         jb = ij 
            id = jb - ib
            if (id .eq. iblock(1,k,ipt)) go to 30 
            ipt = ipt + 1
            iblock(1,k,ipt) = id
            iblock(3,k,ipt) = 0
 30         iblock(3,k,ipt) = iblock(3,k,ipt) + 1 
 35      continue
         lbhb(k) = ipt
c
c ... split zero diagonal block into super and sub diagonals.
c
         jlim = iblock(3,k,2) 
         do 40 j = 1,jlim
            jd = jcnew(k,j)
            if (jd .lt. 0) go to 45
            iblock(3,k,1) = iblock(3,k,1) + 1
            iblock(3,k,2) = iblock(3,k,2) - 1
 40      continue
         j = jlim + 1
 45      iblock(2,k,2) = j
c
c ... form starting positions.
c
         jlim = lbhb(k)
         if (jlim .le. 2) go to 55
         iblock(2,k,3) = 1
         if (jlim .le. 3) go to 55
         do 50 j = 4,jlim
 50      iblock(2,k,j) = iblock(2,k,j-1) + iblock(3,k,j-1)
 55      ist = ied + 1
 60   continue
      return
      end 
      function determ (n,tri,xlmda)
c
c     determ computes the determinant of a symmetric
c     tridiagonal matrix given by tri. det(tri - xlmda*i) = 0
c
c ... parameters -- 
c
c          n      order of tridiagonal system
c          tri    symmetric tridiagonal matrix of order n
c          xlmda  argument for characteristic equation
c
c ... specifications for parameters
c
      dimension tri(2,1)
c
      nm1 = n - 1
      d2 = tri(1,n) - xlmda
      d1 = d2 * (tri(1,nm1) - xlmda) - tri(2,n)
         if (n .eq. 2) go to 20
c
c ... beginning of loop
c
      do 10 l = nm1,2,-1
         d3 = d2
         d2 = d1
         d1 = (tri(1,l-1) - xlmda) * d2 - d3 * tri(2,l)
   10 continue
c
c ... determinant computed
c
   20 determ = d1
c
      return
      end 
      subroutine detsym (ndim,maxnzz,coef,jcoef,nn,isymm)
c
c ... detsym determines if the matrix is symmetric.
c     (purdue storage format) 
c
c ... parameters -- 
c
c         ndim     row dimension of coef in defining routine
c         maxnz    number of columns in coef
c         coef     array of matrix nonzeros
c         jcoef    array of matrix column numbers 
c         n        order of matrix (= nn)
c         isymm    symmetry switch.  upon output, 
c                   isymm = 0  if matrix is symmetric
c                         = 1  if matrix is nonsymmetric
c
c ... specifications for parameters
c
      dimension coef(ndim,2)
      integer   jcoef(ndim,2) 
c
      n = nn
      maxnz = maxnzz
      isymm = 0
      if (maxnz .le. 1) return
      do 20 i = 1,n 
         do 15 j = 2,maxnz
            jcol = jcoef(i,j) 
            if (jcol .eq. i) go to 15
            val = coef(i,j)
            do 10 jj = 2,maxnz
               jcol1 = jcoef(jcol,jj)
               if (jcol1 .ne. i) go to 10
               val1 = coef(jcol,jj)
               if (val1 .eq. val) go to 15
               isymm = 1
               return
 10         continue
            isymm = 1
            return
 15      continue
 20   continue
      return
      end 
      subroutine echall (n,iparm,rparm,icall,icallr,ier)
c
c ... echall initializes the package common blocks from the 
c ... information contained in iparm and rparm.  echall also
c ... prints the values of all parameters in iparm and rparm.
c
c ... parameters -- 
c
c          iparm
c           and
c          rparm  arrays of parameters specifying options and
c                    tolerances
c          icall  indicator of which parameters are being printed
c                    icall = 1,  initial parameters
c                          = 2,  final parameters 
c          icallr  indicator of calling routine
c                          = 1,  called from nspcg
c                          = 2,  called from accelerator
c
c ... specifications for parameters
c
c
c *** begin -- package common 
c
      common / itcom1 / in, itmax, level, nout, ns1, ns2, ns3,
     a      iplr, iqlr, ntest, is, iacel, idgts, nbl1d, nbl2d
      logical           halt, maxadp, minadp, maxadd, minadd
      common / itcom2 / halt, maxadp, minadp, maxadd, minadd
      common / itcom3 / alpha, beta, zeta, emax, emin, pap, 
     b                  alphao, gamma, sigma, rr, rho, dkq, dkm1,
     b                  ff, rqmin, rqmax, stptst, udnm, ubarnm,
     b                  bnorm, bnorm1
      common / itcom4 / keygs, srelpr, keyzer
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / itcom8 / ainf
c
c *** end   -- package common 
c
      logical erflag
      integer   iparm(25)
      dimension rparm(16)
      character*6 inames(25), rnames(16)
      data naiprm, narprm / 11,12 /
      data inames / 'ntest', 'itmax', 'level', 'nout', 'idgts',
     a              'maxadp', 'minadp', 'iomgad', 'ns1', 'ns2', 'ns3',
     a              'nstore', 'iscale', 'iperm', 'ifact', 'lvfill',
     a              'ltrunc', 'ipropa', 'kblsz', 'nbl2d', 'ifctv',
     a              'iqlr', 'isymm', 'ielim', 'ndeg' /
      data rnames / 'zeta', 'emax', 'emin', 'ff', 'fff', 'timit',
     a              'digit1', 'digit2', 'omega', 'alphab', 'betab',
     a              'specr', 'timfac', 'timtot', 'tol', 'ainf' /
c
      if (icall .ne. 1) go to 20
c
c handle accelerator parameters ...
c
      ntest  = iparm(1)
      itmax  = iparm(2)
      level  = iparm(3)
      nout   = iparm(4)
      idgts  = iparm(5)
      maxad  = iparm(6)
      maxadd = (maxad .eq. 1) 
      minad  = iparm(7)
      minadd = (minad .eq. 1) 
      maxadp = maxadd
      minadp = minadd
      iomgad = iparm(8)
      omgadp = (iomgad .eq. 1)
      ns1    = iparm(9)
      ns2    = iparm(10)
      ns3    = iparm(11)
      iqlr   = iparm(22)
      iplr   = iqlr 
c
      zeta   = rparm(1)
      emax   = rparm(2)
      emin   = rparm(3)
      ff     = rparm(4)
      fff    = rparm(5)
      timit  = rparm(6)
      digit1 = rparm(7)
      digit2 = rparm(8)
      omega  = rparm(9)
      alphab = rparm(10)
      betab  = rparm(11)
      specr  = rparm(12)
c
      erflag = .false.
      erflag = erflag .or. ntest .lt. 1 .or. ntest .gt. 10
      erflag = erflag .or. itmax .le. 0 
      erflag = erflag .or. maxad .lt. 0 .or. maxad .gt. 1
      erflag = erflag .or. minad .lt. 0 .or. minad .gt. 1
      erflag = erflag .or. ns1 .lt. 0
      erflag = erflag .or. ns2 .lt. 0
      erflag = erflag .or. emax .lt. 0.0
      erflag = erflag .or. emin .lt. 0.0
      erflag = erflag .or. ff .le. 0.0 .or. ff .gt. 1.0
      if (erflag) go to 999
c
c ... test if eps is too small
c
      temp = 500.0*srelpr
      if (zeta .ge. temp) go to 150
      ier = 2
      call ershow (ier,'echall')
      zeta = temp
      rparm(1) = temp
c
c ... verify n
c
 150  if (n .gt. 0 ) go to 200
      ier = -1
      call ershow (ier,'echall')
      return
c
c now handle preconditioner parameters ...
c
 200  if (icallr .eq. 2) go to 50
      nstore = iparm(12)
      iscale = iparm(13)
      iperm  = iparm(14)
      ifact  = iparm(15)
      lvfill = iparm(16)
      ltrunc = iparm(17)
      ipropa = iparm(18)
      nbl1d  = iparm(19)
      nbl2d  = iparm(20)
      ifctv  = iparm(21)
      iqlr   = iparm(22)
      isymm  = iparm(23)
      ndeg   = iparm(25)
      ainf   = rparm(16)
c
      if (nbl1d .eq. -1) nbl1d = n
      if (nbl2d .eq. -1) nbl2d = n
      kblsz = nbl1d 
      erflag = .false.
      erflag = erflag .or. iqlr .lt. 0 .or. iqlr .gt. 3
      erflag = erflag .or. ipropa .lt. 0 .or. ipropa .gt. 3 
      if (erflag) go to 999
c
c
c ... initialize rest of common variables
c
 50   halt   = .false.
      stptst= 0.0
      udnm   = 1.0
      in     = 0
c
c prepare to do output ...
c
      if (level .le. 2) return
      write (nout,15)
 15   format (/5x,'initial iterative parameters') 
      go to 30
c
 20   if (level .le. 2) return
      write (nout,25)
 25   format (/5x,'final iterative parameters')
c
 30   if (icallr .eq. 2) go to 305
      write (nout,301)
 301  format (5x,'preprocessor and preconditioner parameters')
      ibip = naiprm + 1
      ieip = 25
      ibrp = narprm + 1
      ierp = 16
      go to 300
 305  write (nout,302)
 302  format (5x,'general and acceleration parameters')
      ibip = 1
      ieip = naiprm 
      ibrp = 1
      ierp = narprm 
c
 300  write (nout,35) (i,iparm(i),inames(i),i=ibip,ieip)
 35   format (10x,'iparm(',i2,') =',i15,4x,'(',a6,')'  )
      write (nout,40) (i,rparm(i),rnames(i),i=ibrp,ierp)
 40   format (10x,'rparm(',i2,') =',e15.8,4x,'(',a6,')'  )
      return
c
c error returns ... 
c
c inadmissible option ...
 999  ier = -10
      call ershow (ier,'echall')
      return
      end 
      function eigvss (n,tri,start,end,icode,ier) 
c
c ... eigvss computes a selected eigenvalue of a symmetric
c     tridiagonal matrix for conjugate gradient acceleration.
c     modified imsl routine zbrent used.
c
c ... parameters -- 
c
c          n      order of tridiagonal system
c          tri    symmetric tridiagonal matrix of order n
c          start  initial lower bound of interval containing root
c          end    initial upper bound of interval containing root
c          icode  operation key
c                   = 1   minimum eigenvalue sought
c                   = 2   maximum eigenvalue sought
c          ier    error flag
c
c ... specifications for parameters
c
      dimension tri(2,1)
c
c *** begin -- package common 
c
      common / itcom1 / in, itmax, level, nout, ns1, ns2, ns3,
     a      iplr, iqlr, ntest, is, iacel, idgts, nbl1d, nbl2d
      logical           halt, maxadp, minadp, maxadd, minadd
      common / itcom2 / halt, maxadp, minadp, maxadd, minadd
      common / itcom3 / alpha, beta, zeta, emax, emin, pap, 
     b                  alphao, gamma, sigma, rr, rho, dkq, dkm1,
     b                  ff, rqmin, rqmax, stptst, udnm, ubarnm,
     b                  bnorm, bnorm1
      common / itcom4 / keygs, srelpr, keyzer
c
c *** end   -- package common 
c
      eigvss= 0.0
      itmp = int (-alog10 (abs (zeta))) 
      nsig = max0 (itmp,4)
      maxfn = max0 (itmax,50) 
      eps= 0.0
      a = start
      b = end
      call zbrent (n,tri,eps,nsig,a,b,maxfn,ier)
      if (icode .eq. 1) eigvss = amax1 (a,b)
      if (icode .eq. 2) eigvss = amin1 (a,b)
c
      return
      end 
      subroutine elim1 (nn,ndim,maxnzz,jcoef,coef,rhs,wksp,toll)
c
c ... elim1 removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c     purdue format.
c
c ... parameter list --
c
c         n       dimension of matrix ( = nn)
c         ndim    row dimension of arrays jcoef and coef in the
c                    calling program
c         maxnz   maximum number of nonzero entries per row (=maxnzz) 
c         jcoef   integer array of matrix representation
c         coef    array of sparse matrix representation
c         rhs     right hand side of matrix problem
c         wksp    wksp array of length n
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      integer   jcoef(ndim,1) 
      dimension coef(ndim,1), rhs(1), wksp(1)
c
      n = nn
      maxnz = maxnzz
      tol = toll
      if (n .le. 0  .or.  maxnz .lt. 2) return
c
c ... find maximum off-diagonal elements in absolute value. 
c
      do 10 i = 1,n 
 10   wksp(i)= 0.0
      do 20 j = 2,maxnz
         do 15 i = 1,n
 15      wksp(i) = wksp(i) + abs (coef(i,j))
 20   continue
      do 25 i = 1,n 
 25   wksp(i) = wksp(i) / abs(coef(i,1))
c
c ... eliminate desired rows and columns.
c
      do 35 i = 1,n 
         if (wksp(i) .gt. tol) go to 35 
         rhs(i) = rhs(i)/coef(i,1)
         coef(i,1) = 1.0
         do 30 j = 2,maxnz
            coef(i,j)= 0.0
            jcoef(i,j) = i
 30      continue
 35   continue
      do 45 j = 2,maxnz
         do 40 i = 1,n
            jcol = jcoef(i,j) 
            if (wksp(jcol) .gt. tol) go to 40
            rhs(i) = rhs(i) - coef(i,j)*rhs(jcol) 
            coef(i,j)= 0.0
            jcoef(i,j) = i
 40      continue
 45   continue
      return
      end 
      subroutine elim2 (nn,ndim,maxnzz,jcoef,coef,rhs,wksp,toll)
c
c ... elim2 removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c     symmetric diagonal format.
c
c ... parameter list --
c
c         n       dimension of matrix ( = nn)
c         ndim    row dimension of array coef in the
c                    calling program
c         maxnz   number of diagonals stored
c         jcoef   integer vector of diagonal numbers
c         coef    array of sparse matrix representation
c         rhs     right hand side of matrix problem
c         wksp    wksp array of length n
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      integer   jcoef(1)
      dimension coef(ndim,1), rhs(1), wksp(1)
c
      n = nn
      maxnz = maxnzz
      tol = toll
      if (n .le. 0  .or.  maxnz .lt. 2) return
c
c ... find maximum off-diagonal elements in absolute value. 
c
      do 10 i = 1,n 
 10   wksp(i)= 0.0
      do 25 j = 2,maxnz
         ind = jcoef(j)
         len = n - ind
         do 15 i = 1,len
 15      wksp(i) = wksp(i) + abs (coef(i,j))
         do 20 i = 1,len
 20      wksp(i+ind) = wksp(i+ind) + abs (coef(i,j))
 25   continue
      do 30 i = 1,n 
 30   wksp(i) = wksp(i) / abs(coef(i,1))
c
c ... eliminate desired rows and columns.
c
      do 50 i = 1,n 
         if (wksp(i) .gt. tol) go to 50 
         rhs(i) = rhs(i)/coef(i,1)
         coef(i,1) = 1.0
         do 40 j = 2,maxnz
            jcol = jcoef(j)
            iback = i - jcol
            iforw = i + jcol
            if (iforw .gt. n) go to 35
            if (wksp(iforw) .le. tol) go to 35
            rhs(iforw) = rhs(iforw) - coef(i,j)*rhs(i)
 35         if (iback .lt. 1) go to 40
            rhs(iback) = rhs(iback) - coef(iback,j)*rhs(i)
            coef(iback,j)= 0.0
 40      continue
         do 45 j = 2,maxnz
 45      coef(i,j)= 0.0
 50   continue
      return
      end 
      subroutine elim3 (nn,ndim,maxnzz,jcoef,coef,rhs,wksp,toll)
c
c ... elim3 removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c     nonsymmetric diagonal format.
c
c ... parameter list --
c
c         n       dimension of matrix ( = nn)
c         ndim    row dimension of array coef in the
c                    calling program
c         maxnz   number of diagonals stored
c         jcoef   integer vector of diagonal numbers
c         coef    array of sparse matrix representation
c         rhs     right hand side of matrix problem
c         wksp    wksp array of length n
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      integer   jcoef(1)
      dimension coef(ndim,1), rhs(1), wksp(1)
c
      n = nn
      maxnz = maxnzz
      tol = toll
      if (n .le. 0  .or.  maxnz .lt. 2) return
c
c ... find maximum off-diagonal elements in absolute value. 
c
      do 10 i = 1,n 
 10   wksp(i)= 0.0
      do 20 j = 2,maxnz
         ind = jcoef(j)
         ist1 = max0(1,1 - ind)
         ist2 = min0(n,n - ind)
         do 15 i = ist1,ist2
 15      wksp(i) = wksp(i) + abs (coef(i,j))
 20   continue
      do 25 i = 1,n 
 25   wksp(i) = wksp(i) / abs(coef(i,1))
c
c ... eliminate desired rows and columns.
c
      do 35 i = 1,n 
         if (wksp(i) .gt. tol) go to 35 
         rhs(i) = rhs(i)/coef(i,1)
         coef(i,1) = 1.0
         do 30 j = 2,maxnz
 30      coef(i,j)= 0.0
 35   continue
      do 45 i = 1,n 
         if (wksp(i) .gt. tol) go to 45 
         do 40 j = 2,maxnz
            inew = i - jcoef(j)
            if (inew .lt. 1 .or. inew .gt. n) go to 40
            rhs(inew) = rhs(inew) - coef(inew,j)*rhs(i)
            coef(inew,j)= 0.0 
 40      continue
 45   continue
      return
      end 
      subroutine elim4 (mm,np,ia,ja,a,rhs,wksp,toll)
c
c ... elim4 removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c     symmetric sparse format.
c
c ... parameter list --
c
c         m       number of partitions
c         np      pointer vector to partitions
c         ia      vector of i values
c         ja      vector of j values
c         a       vector of coefficients
c         rhs     right hand side of matrix problem
c         wksp    wksp vector of length n (2n if keygs = 1) 
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      integer ia(1), ja(1), np(2)
      dimension a(1), rhs(1), wksp(1)
c
c *** begin -- package common 
c
      common / itcom4 / keygs, srelpr, keyzer
c
c *** end   -- package common 
c
      m = mm
      n = np(2) - 1 
      nz = np(m+1) - 1
      tol = toll
      np1 = n + 1
c
c ... find sum of absolute values of off-diagonal coefficients.
c
      do 10 i = 1,n 
 10   wksp(i)= 0.0
      if (keygs .eq. 1) go to 30
      do 25 k = 2,m 
         ist = np(k)
         ied = np(k+1) - 1
cdir$ ivdep
         do 15 i = ist,ied
 15      wksp(ia(i)) = wksp(ia(i)) + abs(a(i))
cdir$ ivdep
         do 20 i = ist,ied
 20      wksp(ja(i)) = wksp(ja(i)) + abs(a(i))
 25   continue
      go to 50
 30   do 45 k = 2,m 
         ist = np(k)
         ied = np(k+1) - 1
         len = ied - ist + 1
         call vgathr (len,wksp,ia(ist),wksp(n+1)) 
         do 35 i = ist,ied
 35      wksp(i-ist+1+n) = wksp(i-ist+1+n) + abs(a(i))
         call vscatr (len,wksp(n+1),ia(ist),wksp) 
         call vgathr (len,wksp,ja(ist),wksp(n+1)) 
         do 40 i = ist,ied
 40      wksp(i-ist+1+n) = wksp(i-ist+1+n) + abs(a(i))
         call vscatr (len,wksp(n+1),ja(ist),wksp) 
 45   continue
 50   do 55 i = 1,n 
 55   wksp(i) = wksp(i) / abs(a(i))
c
c ... eliminate desired rows and columns.
c
      do 70 l = 1,n 
         if (wksp(l) .gt. tol) go to 70 
         rhs(l) = rhs(l)/a(l) 
         a(l) = 1.0 
         do 60 k = np1,nz
            i = ia(k)
            j = ja(k)
            if (i .eq. l .and. wksp(j) .gt. tol)
     a              rhs(j) = rhs(j) - a(k)*rhs(i) 
            if (j .ne. l) go to 60
            rhs(i) = rhs(i) - a(k)*rhs(j)
            a(k) = 0.0
 60      continue
         do 65 k = np1,nz
            if (ia(k) .eq. l) a(k) = 0.0
 65      continue
 70   continue
      return
      end 
      subroutine elim5 (mm,np,ia,ja,a,rhs,wksp,toll)
c
c ... elim5 removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c     nonsymmetric sparse format.
c
c ... parameter list --
c
c         m       number of partitions
c         np      pointer vector to partitions
c         ia      vector of i values
c         ja      vector of j values
c         a       vector of coefficients
c         rhs     right hand side of matrix problem
c         wksp    wksp vector of length n (2n if keygs = 1) 
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      integer ia(1), ja(1), np(2)
      dimension a(1), rhs(1), wksp(1)
c
c *** begin -- package common 
c
      common / itcom4 / keygs, srelpr, keyzer
c
c *** end   -- package common 
c
      m = mm
      n = np(2) - 1 
      nz = np(m+1) - 1
      tol = toll
c
c ... find sum of absolute values of off-diagonal coefficients.
c
      do 10 i = 1,n 
 10   wksp(i)= 0.0
      if (keygs .eq. 1) go to 25
      do 20 k = 2,m 
         ist = np(k)
         ied = np(k+1) - 1
cdir$ ivdep
         do 15 i = ist,ied
 15      wksp(ia(i)) = wksp(ia(i)) + abs(a(i))
 20   continue
      go to 40
 25   do 35 k = 2,m 
         ist = np(k)
         ied = np(k+1) - 1
         len = ied - ist + 1
         call vgathr (len,wksp,ia(ist),wksp(n+1)) 
         do 30 i = ist,ied
 30      wksp(i-ist+1+n) = wksp(i-ist+1+n) + abs(a(i))
         call vscatr (len,wksp(n+1),ia(ist),wksp) 
 35   continue
 40   do 45 i = 1,n 
 45   wksp(i) = wksp(i) / abs(a(i))
c
c ... eliminate desired rows and columns.
c
      do 50 i = 1,n 
         if (wksp(i) .gt. tol) go to 50 
         rhs(i) = rhs(i)/a(i) 
         a(i) = 1.0 
 50   continue
      np1 = n + 1
      do 55 k = np1,nz
         if (wksp(ia(k)) .le. tol) a(k) = 0.0
 55   continue
      do 60 k = np1,nz
         j = ja(k)
         if (wksp(j) .gt. tol) go to 60 
         i = ia(k)
         rhs(i) = rhs(i) - a(k)*rhs(j)
         a(k) = 0.0 
 60   continue
      return
      end 
      subroutine ershow (ierr,iname)
c
c ... ershow prints an appropriate error message for the error
c     numbered ier. 
c
c ... parameters -- 
c
c        ier     error number (input)
c                 .gt. 0     for warning errors
c                 .lt. 0     for fatal errors
c        iname   routine name in which error occurred
c
c ... specifications for parameters
c
      character*10 iname
c
c *** begin -- package common 
c
      common / itcom1 / in, itmax, level, nout, ns1, ns2, ns3,
     a      iplr, iqlr, ntest, is, iacel, idgts, nbl1d, nbl2d
      common / itcom3 / alpha, beta, zeta, emax, emin, pap, 
     b                  alphao, gamma, sigma, rr, rho, dkq, dkm1,
     b                  ff, rqmin, rqmax, stptst, udnm, ubarnm,
     b                  bnorm, bnorm1
      common / itcom4 / keygs, srelpr, keyzer
c
c *** end   -- package common 
c
      character*80  fmess(20), wmess(6) 
      data fmess(1)  / 'nonpositive matrix size n' /
      data fmess(2)  / 'insufficient real workspace' /
      data fmess(3)  / 'insufficient integer workspace' /
      data fmess(4)  / 'nonpositive diagonal element' /
      data fmess(5)  / 'nonexistent diagonal element' /
      data fmess(6)  / 'a is not positive definite' /
      data fmess(7)  / 'q is not positive definite' /
      data fmess(8)  / 'unable to permute matrix as requested' /
      data fmess(9)  / 'mdim not large enough to expand matrix' /
      data fmess(10) / 'inadmissible parameter encountered' /
      data fmess(11) / 'incorrect storage mode for block method' /
      data fmess(12) / 'zero pivot encountered during factorization' /
      data fmess(13) / 'breakdown in direction vector calculation' /
      data fmess(14) / 'breakdown in attempt to perform rotation' /
      data fmess(15) / 'breakdown in iterate calculation' / 
      data fmess(16) / 'unimplemented combination of parameters' /
      data fmess(17) / 'error in computing preconditioning polynomial' /
      data fmess(18) / 'unable to perform eigenvalue estimation' /
      data fmess(19) / 'iterative method has gone to sleep' /
      data fmess(20) / 'unknown error' /
      data wmess(1)  / 'failure to converge in itmax iterations' /
      data wmess(2)  / 'zeta too small' /
      data wmess(3)  / 'no convergence in maxfn iterations in zbrent' /
      data wmess(4)  / 'f(a) and f(b) have the same sign in zbrent' / 
      data wmess(5)  / 'negative pivot encountered in factorization' /
      data wmess(6)  / 'unknown warning' /
c
      ier = ierr
      if (ier .eq. 0) return
      if (ier .lt. 0  .and.  level .lt. 0) return 
      if (ier .gt. 0  .and.  level .lt. 1) return 
      if (ier .lt. -19) ier = -20
      if (ier .gt.   5) ier =   6
      if (ier .lt. 0) write (nout,10)
 10   format (//1x,60('*') /
     a          1x,18('*'),' f a t a l    e r r o r ',18('*') /
     a          1x,60('*') /) 
      if (ier .gt. 0) write (nout,20)
 20   format (//1x,60('*') /
     a          1x,16('*'),' w a r n i n g    e r r o r ',16('*') /
     a          1x,60('*') /) 
      write (nout,23) iname
 23   format (1x,'routine ',a10)
      inum = iabs(ier)
      if (ier .gt. 0) go to 30
c
c ... print out fatal errors. 
c
      write (nout,25) fmess(inum)
 25   format (1x,a80)
      go to 999
c
c ... print out warning errors.
c
 30   write (nout,25) wmess(inum)
      if (inum .ne. 2) go to 999
      temp = 500.0*srelpr
      write (nout,35) zeta, srelpr, temp
 35   format (1x,'rparm(1) =',e10.3,' (zeta)'
     a     / 1x, 'a value this small may hinder convergence'
     a     / 1x, 'since machine precision srelpr = ',e10.3
     a     / 1x, 'zeta reset to ',e10.3)
c
c ... print ending line.
c
 999  write (nout,1000)
 1000 format (/1x,60('*')/)
      return
      end 
      subroutine filln (maxnz,jcoef)
c
c ... filln determines the fill-in diagonals for nonsymmetric
c     diagonal storage.
c
c ... parameters -- 
c
c        maxnz   upon input, the number of diagonals
c                upon output, the number of diagonals with fill-in
c        jcoef   upon input, the diagonal numbers 
c                upon output, the diagonal numbers with fill-in
c
c ... specifications for parameters
c
      integer jcoef(2)
c
      maxn = maxnz
      do 20 j1 = 1,maxnz
         do 15 j2 = 1,maxnz
            jd = jcoef(j1) + jcoef(j2)
            if (jcoef(j1)*jcoef(j2) .ge. 0) go to 15
            do 10 j3 = 1,maxn 
               if (jcoef(j3) .eq. jd) go to 15
 10         continue
            maxn = maxn + 1
            jcoef(maxn) = jd
 15      continue
 20   continue
      maxnz = maxn
      return
      end 
      subroutine fills (maxt,jt)
c
c ... fills determines the fill-in diagonals for symmetric
c     diagonal storage.
c
c ... parameters -- 
c
c        maxt    upon input, the number of diagonals in the 
c                 upper triangle
c                upon output, the number of diagonals in the
c                 upper triangle with fill-in
c        jt      upon input, the diagonal numbers in the upper
c                 triangle
c                upon output, the diagonal numbers in the upper
c                 triangle with fill-in 
c
c ... specifications for parameters
c
      integer jt(1) 
c
      maxn = maxt
      do 20 j1 = 1,maxt
         do 15 j2 = 1,maxt
            jd = jt(j1) - jt(j2)
            if (jd .le. 0) go to 15
            do 10 j3 = 1,maxn 
               if (jt(j3) .eq. jd) go to 15
 10         continue
            maxn = maxn + 1
            jt(maxn) = jd
 15      continue
 20   continue
      maxt = maxn
      return
      end 
      subroutine fillnp (ndim,nn,maxcc,jc,c,mwidth,ier)
c
c ... fillnp determines the fill-in structure.
c     (purdue storage, nonsymmetric matrix)
c
c ... parameters -- 
c
c          ndim   row dimension of jc and c arrays
c          n      order of system (= nn)
c          maxc   upon input, maxc is the number of columns in
c                  the c array
c                 upon output, maxc is the number of columns in
c                  the c array with fill-in
c          jc     integer array of active size n by maxc giving the
c                  column numbers of the corresponding elements in c
c          c      array of active size n by maxc giving the 
c                  coefficients of the off-diagonal elements
c          mwidth maximum column width to be allowed for fill-in
c          ier    error code
c                  =    0  no errors detected
c                  =   -2  mwidth too small to accomodate fill-in
c
c ... specifications for parameters
c
      integer   jc(ndim,1)
      dimension c(ndim,1)
c
c
      n = nn
      maxc = maxcc
      maxu = maxc
c
      if (maxc .lt. 1) return 
      nm1 = n - 1
      do 45 k = 1,nm1
         kp1 = k + 1
         do 40 j1 = 1,maxc
         do 35 i = kp1,n
            if (jc(i,j1) .ne. k) go to 35
            do 30 j2 = 1,maxc 
               j = jc(k,j2)
               if (j .le. k .or. j .eq. i) go to 30
               do 10 j3 = 1,maxu
                  if (j .eq. iabs(jc(i,j3))) go to 30
 10            continue
               do 15 j3 = 1,maxu
                  if (jc(i,j3) .ne. i) go to 15
                  jc(i,j3) = -j
                  go to 30
 15            continue
               maxu = maxu + 1
               if (maxu .le. mwidth) go to 20
               ier = -2
               return
 20            do 25 ii = 1,n 
                  jc(ii,maxu) = ii
                  c(ii,maxu)= 0.0
 25            continue
               jc(i,maxu) = -j
 30         continue
 35      continue
 40      continue
 45   continue
c
c ... decode new elements of jt, jb.
c
      do 55 j = 1,maxu
         do 50 i = 1,n
 50      jc(i,j) = iabs(jc(i,j))
 55   continue
      maxcc = maxu
      return
      end 
      subroutine fillsp (ndim,nn,maxtt,jt,t,mwidth,ier)
c
c ... fillsp determines the fill-in structure.
c     (purdue storage, symmetric matrix)
c
c ... parameters -- 
c
c          ndim   row dimension of t and jt arrays
c          n      order of system (= nn)
c          maxt   upon input, maxt is the number of columns in
c                  the t array
c                 upon output, maxt is the number of columns in
c                  the t array with fill-in
c          jt     integer array of active size n by maxt giving the
c                  column numbers of the corresponding elements in t
c          t      array of active size n by maxt giving the 
c                  coefficients of the upper triangle of the matrix
c          mwidth maximum column width of jt and t to be allowed
c          ier    error code
c                  =   0     no error detected
c                  =  -2     mwidth too small to store factor
c
c ... specifications for parameters
c
      dimension t(ndim,1)
      integer   jt(ndim,1)
c
c
      n = nn
      maxt = maxtt
      maxu = maxt
      ier = 0
c
      if (maxt .lt. 1) return 
      nm1 = n - 1
      do 40 k = 1,nm1
         do 35 j1 = 1,maxt
            jcol1 = jt(k,j1)
            if (jcol1 .le. 0 .or. jcol1 .eq. k) go to 35
            do 30 j2 = 1,maxt 
               jcol2 = jt(k,j2)
               if (jcol2 .le. 0 .or. jcol2 .eq. k) go to 30 
               if (jcol2 .le. jcol1) go to 30
               do 10 j3 = 1,maxu
                  if (jcol2 .eq. iabs(jt(jcol1,j3))) go to 30
 10            continue
               do 15 j3 = 1,maxu
                  if (jt(jcol1,j3) .ne. jcol1) go to 15
                  jt(jcol1,j3) = -jcol2 
                  go to 30
 15            continue
               maxu = maxu + 1
               if (maxu .le. mwidth) go to 20
               ier = -2
               return
 20            do 25 i = 1,n
                  jt(i,maxu) = i
                  t(i,maxu) = 0.0
 25            continue
               jt(jcol1,maxu) = -jcol2
 30         continue
 35      continue
 40   continue
c
c ... decode new elements of jt.
c
      do 50 j = 1,maxu
         do 45 i = 1,n
 45      jt(i,j) = iabs(jt(i,j))
 50   continue
      maxtt = maxu
      return
      end 
      subroutine ibfcn1 (lddd,ldtt,n,jd,jt,d,t,ncol,nci,
     a                   iblock,lbhb,iunif,ipropa,ipt,
     a                   omega,wksp,ier)
c
c ... ibfcn1 does an incomplete block factorization of the matrix
c     contained in d and t (version 1, unmodified).
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ic (version 1) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer array of size ncolor by whatever 
c                   giving the diagonal block diagonal numbers for
c                   each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         jt       integer array of size ncolor by whatever 
c                   giving the off-diagonal block diagonal numbers
c                   for each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         ncolor   number of distinct block sizes 
c                   ncolor = 1 if iunif = 1.
c         nci      integer vector of length ncolor, giving the number 
c                   of nodes for each distinct block size.
c                   if iunif = 1, nci(1) is the constant block size.
c         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         lbhb     integer vector of size ncolor giving the number
c                   of diagonal blocks for each distinct block size.
c                   if iunif = 1, lbhb is of length 1.
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0 
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jd(ncol,1), jt(ncol,1), nci(1), lbhb(1),
     a          iblock(3,ncol,2)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   unif, propa
c
      ldd = lddd
      ldt = ldtt
      ncolor = ncol 
      unif = iunif .eq. 1
      propa = ipropa .eq. 1
c
c ... define various constants.
c
      if (unif) go to 15
      klim = ncolor 
      go to 20
 15   kblsz = nci(1)
      na = kblsz
      nb = kblsz
      nc = kblsz
      ii = 1
      kk = 1
      jlim = lbhb(1)
      llim = jlim
      klim = n/kblsz
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      ma = ndt + ndb + 1
c
c ... start factorization.
c
 20   do 95 k = 1,klim
         if (unif) go to 25
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         ma = ndt + ndb + 1
         go to 30
 25      ist = (k - 1)*kblsz + 1
 30      call bdfac (ldd,na,na,ndt,ndb,d(ist,1),1)
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,ndb,wksp,1)
         if (k .eq. klim .or. jlim .le. 2) go to 95
         do 90 i = k+1,klim
            if (unif) go to 35
            ii = i
            llim = lbhb(i)
 35         if (llim .le. 2) go to 90
            do 40 l = 3,llim
               jcol = i + iblock(1,ii,l)
               if (jcol .eq. k) go to 45
 40         continue
            go to 90
 45         mc = iblock(3,ii,l)
            if (unif) go to 50
            nc = ipt(i+1) - ipt(i)
            incc = ipt(k) - ipt(i)
            go to 55
 50         incc = (k - i)*kblsz
 55         istc = ist - incc 
            jstc = iblock(2,ii,l)
            do 85 j = 3,jlim
               jcol = k + iblock(1,kk,j)
               if (jcol .le. k) go to 85
               jdiff = jcol - i
               if (jdiff .ne. 0 .and. propa) go to 85
               do 60 m = 1,llim
                  if (iblock(1,ii,m) .eq. jdiff) go to 65
 60            continue
               go to 85
 65            mb = iblock(3,kk,j)
               istb = ist
               jstb = iblock(2,kk,j)
               if (unif) go to 70
               nb = ipt(jcol+1) - ipt(jcol)
               incb = ipt(jcol) - ipt(k)
               go to 75
 70            incb = (jcol - k)*kblsz
 75            incd = incc + incb
               istd = istc
               jstd = iblock(2,ii,m)
               md = iblock(3,ii,m)
               if (m .eq. 1) go to 80
               call t1prod (na,ldt,ldt,ldt,ncolor,na,nc,nb, 
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jt(ii,jstd),wksp,t(istb,jstb),
     a                      t(istc,jstc),t(istd,jstd))
               go to 85
 80            md = md + iblock(3,ii,2) 
               call t1prod (na,ldt,ldt,ldd,ncolor,na,nc,nb, 
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jd(ii,jstd),wksp,t(istb,jstb),
     a                      t(istc,jstc),d(istd,jstd))
 85         continue
 90      continue
 95   continue
      return
      end 
      subroutine ibfcn2 (lddd,ldtt,n,jd,jt,d,t,ncol,nci,
     a                   iblock,lbhb,iunif,ipropa,ipt,
     a                   omega,wksp,ier)
c
c ... ibfcn2 does an incomplete block factorization of the matrix
c     contained in d and t (version 2, unmodified).
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ic (version 2) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer array of size ncolor by whatever 
c                   giving the diagonal block diagonal numbers for
c                   each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         jt       integer array of size ncolor by whatever 
c                   giving the off-diagonal block diagonal numbers
c                   for each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         ncolor   number of distinct block sizes 
c                   ncolor = 1 if iunif = 1.
c         nci      integer vector of length ncolor, giving the number 
c                   of nodes for each distinct block size.
c                   if iunif = 1, nci(1) is the constant block size.
c         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         lbhb     integer vector of size ncolor giving the number
c                   of diagonal blocks for each distinct block size.
c                   if iunif = 1, lbhb is of length 1.
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0 
c
c ... specifications for parameters
c
      integer   ipt(1), jd(ncol,1), jt(ncol,1), nci(1), lbhb(1),
     a          iblock(3,ncol,2)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   unif, propa
c
      ldd = lddd
      ldt = ldtt
      ncolor = ncol 
      unif = iunif .eq. 1
      propa = ipropa .eq. 1
c
c ... define various constants.
c
      if (unif) go to 15
      klim = ncolor 
      go to 20
 15   kblsz = nci(1)
      na = kblsz
      nb = kblsz
      nc = kblsz
      ii = 1
      kk = 1
      jlim = lbhb(1)
      llim = jlim
      klim = n/kblsz
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      ma = ndt + ndb + 1
c
c ... start factorization.
c
 20   do 95 k = 1,klim
         if (unif) go to 25
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         ma = ndt + ndb + 1
         go to 30
 25      ist = (k - 1)*kblsz + 1
 30      call bdfac (ldd,na,na,ndt,ndb,d(ist,1),1)
         call bdinv (ldd,na,na,ndt,ndb,d(ist,1),1)
         if (k .eq. klim .or. jlim .le. 2) go to 95
         do 90 i = k+1,klim
            if (unif) go to 35
            ii = i
            llim = lbhb(i)
 35         if (llim .le. 2) go to 90
            do 40 l = 3,llim
               jcol = i + iblock(1,ii,l)
               if (jcol .eq. k) go to 45
 40         continue
            go to 90
 45         mc = iblock(3,ii,l)
            if (unif) go to 50
            nc = ipt(i+1) - ipt(i)
            incc = ipt(k) - ipt(i)
            go to 55
 50         incc = (k - i)*kblsz
 55         istc = ist - incc 
            jstc = iblock(2,ii,l)
            do 85 j = 3,jlim
               jcol = k + iblock(1,kk,j)
               if (jcol .le. k) go to 85
               jdiff = jcol - i
               if (jdiff .ne. 0 .and. propa) go to 85
               do 60 m = 1,llim
                  if (iblock(1,ii,m) .eq. jdiff) go to 65
 60            continue
               go to 85
 65            mb = iblock(3,kk,j)
               istb = ist
               jstb = iblock(2,kk,j)
               if (unif) go to 70
               nb = ipt(jcol+1) - ipt(jcol)
               incb = ipt(jcol) - ipt(k)
               go to 75
 70            incb = (jcol - k)*kblsz
 75            incd = incc + incb
               istd = istc
               jstd = iblock(2,ii,m)
               md = iblock(3,ii,m)
               if (m .eq. 1) go to 80
               call t1prod (ldd,ldt,ldt,ldt,ncolor,na,nc,nb,
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jt(ii,jstd),d(ist,1),t(istb,jstb),
     a                      t(istc,jstc),t(istd,jstd))
               go to 85
 80            md = md + iblock(3,ii,2) 
               call t1prod (ldd,ldt,ldt,ldd,ncolor,na,nc,nb,
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jd(ii,jstd),d(ist,1),t(istb,jstb),
     a                      t(istc,jstc),d(istd,jstd))
 85         continue
 90      continue
 95   continue
      return
      end 
      subroutine ibfcn3 (lddd,ldtt,n,jd,jt,d,t,ncol,nci,
     a                   iblock,lbhb,iunif,ipropa,ipt,omega,wksp,
     a                   ier) 
c
c ... ibfcn3 does an incomplete block factorization of the matrix
c     contained in d and t (version 1, modified). 
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ic (version 1) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer array of size ncolor by whatever 
c                   giving the diagonal block diagonal numbers for
c                   each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         jt       integer array of size ncolor by whatever 
c                   giving the off-diagonal block diagonal numbers
c                   for each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         ncolor   number of distinct block sizes 
c                   ncolor = 1 if iunif = 1.
c         nci      integer vector of length ncolor, giving the number 
c                   of nodes for each distinct block size.
c                   if iunif = 1, nci(1) is the constant block size.
c         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         lbhb     integer vector of size ncolor giving the number
c                   of diagonal blocks for each distinct block size.
c                   if iunif = 1, lbhb is of length 1.
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0 
c         omega    relaxation factor between 0 and 1.
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jd(ncol,1), jt(ncol,1), nci(1), lbhb(1),
     a          iblock(3,ncol,2)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   unif, propa
c
      ldd = lddd
      ldt = ldtt
      ncolor = ncol 
      unif = iunif .eq. 1
      propa = ipropa .eq. 1
c
c ... define various constants.
c
      if (unif) go to 15
      klim = ncolor 
      go to 20
 15   kblsz = nci(1)
      na = kblsz
      nb = kblsz
      nc = kblsz
      ii = 1
      kk = 1
      jlim = lbhb(1)
      llim = jlim
      klim = n/kblsz
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      ma = ndt + ndb + 1
c
c ... start factorization.
c
 20   do 100 k = 1,klim
         if (unif) go to 25
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         ma = ndt + ndb + 1
         go to 30
 25      ist = (k - 1)*kblsz + 1
 30      call bdfac (ldd,na,na,ndt,ndb,d(ist,1),1)
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,ndb,wksp,1)
         ip1 = na*ma + 1
         ip2 = ip1 + na - 1
         if (k .eq. klim .or. jlim .le. 2) go to 100
         do 95 i = k+1,klim
            if (unif) go to 35
            ii = i
            llim = lbhb(i)
 35         if (llim .le. 2) go to 95
            do 40 l = 3,llim
               jcol = i + iblock(1,ii,l)
               if (jcol .eq. k) go to 45
 40         continue
            go to 95
 45         mc = iblock(3,ii,l)
            if (unif) go to 50
            nc = ipt(i+1) - ipt(i)
            incc = ipt(k) - ipt(i)
            go to 55
 50         incc = (k - i)*kblsz
 55         istc = ist - incc 
            jstc = iblock(2,ii,l)
            do 90 j = 3,jlim
               jcol = k + iblock(1,kk,j)
               if (jcol .le. k) go to 90
               mb = iblock(3,kk,j)
               istb = ist
               jstb = iblock(2,kk,j)
               if (unif) go to 60
               nb = ipt(jcol+1) - ipt(jcol)
               incb = ipt(jcol) - ipt(k)
               go to 65
 60            incb = (jcol - k)*kblsz
 65            incd = incc + incb
               istd = istc
               jdiff = jcol - i
               if (jdiff .ne. 0 .and. propa) go to 85
               do 70 m = 1,llim
                  if (iblock(1,ii,m) .eq. jdiff) go to 75
 70            continue
               go to 85
 75            jstd = iblock(2,ii,m)
               md = iblock(3,ii,m)
               if (m .eq. 1) go to 80
               call t1prod (na,ldt,ldt,ldt,ncolor,na,nc,nb, 
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jt(ii,jstd),wksp,t(istb,jstb),
     a                      t(istc,jstc),t(istd,jstd))
               call tsumn
     a               (na,nc,nb,na,ldt,ldt,ncolor,ma,mb,mc,md,incb,
     a                incc,incd,jd(kk,1),jt(kk,jstb),jt(ii,jstc),
     a                jt(ii,jstd),wksp,t(istb,jstb),t(istc,jstc),
     a                d(istd,1),omega)
               go to 85
 80            md = md + iblock(3,ii,2) 
               call t1prod (na,ldt,ldt,ldd,ncolor,na,nc,nb, 
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jd(ii,jstd),wksp,t(istb,jstb),
     a                      t(istc,jstc),d(istd,jstd))
               call tsumn
     a               (na,nc,nb,na,ldt,ldt,ncolor,ma,mb,mc,md,incb,
     a                incc,incd,jd(kk,1),jt(kk,jstb),jt(ii,jstc),
     a                jd(ii,jstd),wksp,t(istb,jstb),t(istc,jstc),
     a                d(istd,1),omega)
 85            call rowsum (ldt,na,mb,t(istb,jstb),wksp(ip1),1)
               do 87 iii = ip1,ip2
 87            wksp(iii) = omega*wksp(iii)
               call bdsol (ldd,na,na,ndt,ndb,d(ist,1),wksp(ip1),
     a                     wksp(ip1),1) 
               call vsubd (ldt,ncolor,nc,na,mc,t(istc,jstc),
     a                      jt(ii,jstc),d(istd,1),wksp(ip1),incc)
 90         continue
 95      continue
 100  continue
      return
      end 
      subroutine ibfcn4 (lddd,ldtt,n,jd,jt,d,t,ncol,nci,
     a                   iblock,lbhb,iunif,ipropa,ipt,omega,wksp,
     a                   ier) 
c
c ... ibfcn4 does an incomplete block factorization of the matrix
c     contained in d and t (version 2, modified). 
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ic (version 2) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer array of size ncolor by whatever 
c                   giving the diagonal block diagonal numbers for
c                   each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         jt       integer array of size ncolor by whatever 
c                   giving the off-diagonal block diagonal numbers
c                   for each distinct block size.  jd is 1 by whatever
c                   if iunif = 1.
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         ncolor   number of distinct block sizes 
c                   ncolor = 1 if iunif = 1.
c         nci      integer vector of length ncolor, giving the number 
c                   of nodes for each distinct block size.
c                   if iunif = 1, nci(1) is the constant block size.
c         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         lbhb     integer vector of size ncolor giving the number
c                   of diagonal blocks for each distinct block size.
c                   if iunif = 1, lbhb is of length 1.
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0 
c         omega    relaxation factor between 0 and 1.
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jd(ncol,1), jt(ncol,1), nci(1), lbhb(1),
     a          iblock(3,ncol,2)
      dimension d(lddd,2), t(ldtt,1), wksp(1)
      logical   unif, propa
c
      ldd = lddd
      ldt = ldtt
      ncolor = ncol 
      unif = iunif .eq. 1
      propa = ipropa .eq. 1
c
c ... define various constants.
c
      ip1 = n + 1
      if (unif) go to 15
      klim = ncolor 
      do 13 k = 1,ncolor
         ist = ipt(k) + 1
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         ma = ndt + ndb + 1
         call rowsum (ldd,na,ma,d(ist,1),wksp(ist),1)
 13   continue
      go to 20
 15   kblsz = nci(1)
      na = kblsz
      nb = kblsz
      nc = kblsz
      ii = 1
      kk = 1
      jlim = lbhb(1)
      llim = jlim
      klim = n/kblsz
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      ma = ndt + ndb + 1
      call rowsum (ldd,n,ma,d,wksp,1)
c
c ... start factorization.
c
 20   do 100 k = 1,klim
         if (unif) go to 25
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         ma = ndt + ndb + 1
         go to 30
 25      ist = (k - 1)*kblsz + 1
 30      isu = ist + na - 1
         call bdfac (ldd,na,na,ndt,ndb,d(ist,1),1)
         call bdinv (ldd,na,na,ndt,ndb,d(ist,1),1)
         call bmuln (ldd,na,ndt,ndb,d(ist,1),d(ist,2),d(ist,ndt+2),
     a               wksp(ist),wksp(ip1))
         do 31 iii = ist,isu
            if (wksp(iii) .ne. 0.0) go to 31
            ier = -12
            call ershow (ier,'ibfcn4')
            return
 31      continue
         do 33 iii = ist,isu
 33      d(iii,1) = d(iii,1) + omega*(1.0 - wksp(iii-ist+ip1))/
     a                              wksp(iii)
         ip2 = ip1 + na
         if (k .eq. klim .or. jlim .le. 2) go to 100
         do 95 i = k+1,klim
            if (unif) go to 35
            ii = i
            llim = lbhb(i)
 35         if (llim .le. 2) go to 95
            do 40 l = 3,llim
               jcol = i + iblock(1,ii,l)
               if (jcol .eq. k) go to 45
 40         continue
            go to 95
 45         mc = iblock(3,ii,l)
            if (unif) go to 50
            nc = ipt(i+1) - ipt(i)
            incc = ipt(k) - ipt(i)
            go to 55
 50         incc = (k - i)*kblsz
 55         istc = ist - incc 
            jstc = iblock(2,ii,l)
            do 90 j = 3,jlim
               jcol = k + iblock(1,kk,j)
               if (jcol .le. k) go to 90
               mb = iblock(3,kk,j)
               istb = ist
               jstb = iblock(2,kk,j)
               if (unif) go to 60
               nb = ipt(jcol+1) - ipt(jcol)
               incb = ipt(jcol) - ipt(k)
               go to 65
 60            incb = (jcol - k)*kblsz
 65            incd = incc + incb
               istd = istc
               jdiff = jcol - i
               if (jdiff .ne. 0 .and. propa) go to 85
               do 70 m = 1,llim
                  if (iblock(1,ii,m) .eq. jdiff) go to 75
 70            continue
               go to 85
 75            jstd = iblock(2,ii,m)
               md = iblock(3,ii,m)
               if (m .eq. 1) go to 80
               call t1prod (ldd,ldt,ldt,ldt,ncolor,na,nc,nb,
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jt(ii,jstd),d(ist,1),t(istb,jstb),
     a                      t(istc,jstc),t(istd,jstd))
               call tsumn
     a               (na,nc,nb,ldd,ldt,ldt,ncolor,ma,mb,mc,md,incb,
     a                incc,incd,jd(kk,1),jt(kk,jstb),jt(ii,jstc),
     a                jt(ii,jstd),d(ist,1),t(istb,jstb),t(istc,jstc), 
     a                wksp(istd),1.0)
               go to 85
 80            md = md + iblock(3,ii,2) 
               call t1prod (ldd,ldt,ldt,ldd,ncolor,na,nc,nb,
     a                      ma,mb,mc,md,incb,incc,incd,jd(kk,1),
     a                      jt(kk,jstb),jt(ii,jstc),
     a                      jd(ii,jstd),d(ist,1),t(istb,jstb),
     a                      t(istc,jstc),d(istd,jstd))
 85            call rowsum (ldt,na,mb,t(istb,jstb),wksp(ip1),1)
               call bmuln (ldd,na,ndt,ndb,d(ist,1),d(ist,2),
     a                     d(ist,ndt+2),wksp(ip1),wksp(ip2))
               call vsubd (ldt,ncolor,nc,na,mc,t(istc,jstc),
     a                      jt(ii,jstc),wksp(istd),wksp(ip2),incc)
 90         continue
 95      continue
 100  continue
      return
      end 
      subroutine ibfcs1 (lddd,ldtt,nn,jd,jt,d,t,kblszz,
     a                   iblock,lbhb,ipropa,omega,wksp,ier) 
c
c ... ibfcs1 does an incomplete block factorization of the matrix
c     contained in d and t (version 1, unmodified).
c     symmetric diagonal data structure, natural ordering.
c     block ic (version 1) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer vector giving the diagonal numbers
c                   for the diagonal block
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         kblsz    block size 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         lbhb     number of blocks per block row 
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   jd(1), jt(1), iblock(3,3)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   propa
c
      n = nn
      ldd = lddd
      ldt = ldtt
      na = kblszz
      propa = ipropa .eq. 1
      klim = n/na
      ma = iblock(3,1)
      ndt = ma - 1
c
c ... block tridiagonal case. 
c
      if (lbhb .gt. 3) go to 25
      jblkb = iblock(1,3)
      mb = iblock(3,3)
      incb = jblkb*na
      do 20 k = 1,klim
         ist = (k - 1)*na + 1 
         istd = ist + incb
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         if (istd .gt. n) go to 20
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,0,wksp,0)
         call t2prod (na,na,ldt,ldt,ldd,ma,mb,mb,ma,
     a                incb,incb,0,jd,jt,jt,jd,wksp,t(ist,1),t(ist,1), 
     a                d(istd,1))
 20   continue
      return
c
c ... general block structure.
c
 25   do 50 k = 1,klim
         ist = (k - 1)*na + 1 
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         if (k .eq. klim) go to 50
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,0,wksp,0)
         jjlim = min0(lbhb,klim-k+2)
         do 45 jjc = 3,jjlim
            jblkc = iblock(1,jjc)
            jstc = iblock(2,jjc)
            mc = iblock(3,jjc)
            incc = jblkc*na
            istd = ist + incc 
            if (istd .gt. n) go to 45
            do 40 jjb = 3,jjlim
               jblkb = iblock(1,jjb)
               jstb = iblock(2,jjb)
               mb = iblock(3,jjb)
               incb = jblkb*na
               jdiff = jblkb - jblkc
               if (jdiff .lt. 0) go to 40
               if (jdiff .ne. 0 .and. propa) go to 40
               do 30 jjd = 1,jjlim
                  if (jdiff .eq. iblock(1,jjd)) go to 35
 30            continue
               go to 40
 35            jblkd = iblock(1,jjd)
               jstd = iblock(2,jjd)
               md = iblock(3,jjd)
               incd = jblkd*na
               if (jjd .ne. 1) call t2prod
     a               (na,na,ldt,ldt,ldt,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jt(jstd),wksp,t(ist,jstb),t(ist,jstc),
     a                t(istd,jstd))
               if (jjd .eq. 1) call t2prod
     a               (na,na,ldt,ldt,ldd,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jd,wksp,t(ist,jstb),t(ist,jstc),
     a                d(istd,1))
 40         continue
 45      continue
 50   continue
      return
      end 
      subroutine ibfcs2 (lddd,ldtt,nn,jd,jt,d,t,kblszz,
     a                   iblock,lbhb,ipropa,omega,wksp,ier) 
c
c ... ibfcs2 does an incomplete block factorization of the matrix
c     contained in d and t (version 2, unmodified).
c     symmetric diagonal data structure, natural ordering.
c     block ic (version 2) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer vector giving the diagonal numbers
c                   for the diagonal block
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         kblsz    block size 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         lbhb     number of blocks per block row 
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c
c ... specifications for parameters
c
      integer   jd(1), jt(1), iblock(3,3)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   propa
c
      n = nn
      ldd = lddd
      ldt = ldtt
      na = kblszz
      propa = ipropa .eq. 1
      klim = n/na
      ma = iblock(3,1)
      ndt = ma - 1
c
c ... block tridiagonal case. 
c
      if (lbhb .gt. 3) go to 25
      jblkb = iblock(1,3)
      mb = iblock(3,3)
      incb = jblkb*na
      do 20 k = 1,klim
         ist = (k - 1)*na + 1 
         istd = ist + incb
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         call bdinv (ldd,na,na,ndt,0,d(ist,1),0)
         if (istd .gt. n) go to 20
         call t2prod (na,ldd,ldt,ldt,ldd,ma,mb,mb,ma,
     a                incb,incb,0,jd,jt,jt,jd,d(ist,1),t(ist,1),
     a                t(ist,1),d(istd,1))
 20   continue
      return
c
c ... general block structure.
c
 25   do 50 k = 1,klim
         ist = (k - 1)*na + 1 
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         call bdinv (ldd,na,na,ndt,0,d(ist,1),0)
         if (k .eq. klim) go to 50
         jjlim = min0(lbhb,klim-k+2)
         do 45 jjc = 3,jjlim
            jblkc = iblock(1,jjc)
            jstc = iblock(2,jjc)
            mc = iblock(3,jjc)
            incc = jblkc*na
            istd = ist + incc 
            if (istd .gt. n) go to 45
            do 40 jjb = 3,jjlim
               jblkb = iblock(1,jjb)
               jstb = iblock(2,jjb)
               mb = iblock(3,jjb)
               incb = jblkb*na
               jdiff = jblkb - jblkc
               if (jdiff .lt. 0) go to 40
               if (jdiff .ne. 0 .and. propa) go to 40
               do 30 jjd = 1,jjlim
                  if (jdiff .eq. iblock(1,jjd)) go to 35
 30            continue
               go to 40
 35            jblkd = iblock(1,jjd)
               jstd = iblock(2,jjd)
               md = iblock(3,jjd)
               incd = jblkd*na
               if (jjd .ne. 1) call t2prod
     a               (na,ldd,ldt,ldt,ldt,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jt(jstd),d(ist,1),t(ist,jstb),t(ist,jstc),
     a                t(istd,jstd))
               if (jjd .eq. 1) call t2prod
     a               (na,ldd,ldt,ldt,ldd,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jd,d(ist,1),t(ist,jstb),t(ist,jstc),
     a                d(istd,1))
 40         continue
 45      continue
 50   continue
      return
      end 
      subroutine ibfcs3 (lddd,ldtt,nn,jd,jt,d,t,kblszz,
     a                   iblock,lbhb,ipropa,omegaa,wksp,ier)
c
c ... ibfcs3 does an incomplete block factorization of the matrix
c     contained in d and t (version 1, modified). 
c     symmetric diagonal data structure, natural ordering.
c     block ic (version 1) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer vector giving the diagonal numbers
c                   for the diagonal block
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         kblsz    block size 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         lbhb     number of blocks per block row 
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         omega    relaxation factor between 0. and 1.
c                   = 0   no modification
c                   = 1   full modification
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   jd(1), jt(1), iblock(3,3)
      dimension d(lddd,1), t(ldtt,1), wksp(1)
      logical   propa
c
      n = nn
      ldd = lddd
      ldt = ldtt
      na = kblszz
      omega = omegaa
      propa = ipropa .eq. 1
      klim = n/na
      ma = iblock(3,1)
      ndt = ma - 1
c
c ... block tridiagonal case. 
c
      if (lbhb .gt. 3) go to 25
      ip1 = na*ma + 1
      ip2 = ip1 + na - 1
      jblkb = iblock(1,3)
      mb = iblock(3,3)
      incb = jblkb*na
      do 20 k = 1,klim
         ist = (k - 1)*na + 1 
         istd = ist + incb
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         if (istd .gt. n) go to 20
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,0,wksp,0)
         call t2prod (na,na,ldt,ldt,ldd,ma,mb,mb,ma,
     a                incb,incb,0,jd,jt,jt,jd,wksp,t(ist,1),
     a                t(ist,1),d(istd,1))
         call tsum (na,na,ldt,ldt,ma,mb,mb,ma,incb,incb,
     a              0,jd,jt,jt,jd,wksp,t(ist,1),t(ist,1),
     a              d(istd,1),d(istd,1),wksp(ip1),1,omega)
         call rowsum (ldt,na,mb,t(ist,1),wksp(ip1),1)
         do 15 iii = ip1,ip2
 15      wksp(iii) = omega*wksp(iii)
         call bdsol (ldd,na,na,ndt,0,d(ist,1),wksp(ip1),
     a               wksp(ip1),0)
         call vsubdt (ldt,1,na,na,mb,t(ist,1),jt,d(istd,1), 
     a                wksp(ip1),incb)
 20   continue
      return
c
c ... general block structure.
c
 25   ip1 = na*ma + 1
      ip2 = ip1 + na - 1
      do 60 k = 1,klim
         ist = (k - 1)*na + 1 
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         if (k .eq. klim) go to 60
         call mcopy (ldd,na,na,ma,d(ist,1),wksp)
         call bdinv (na,na,na,ndt,0,wksp,0)
         jjlim = min0(lbhb,klim-k+2)
         do 55 jjc = 3,jjlim
            jblkc = iblock(1,jjc)
            jstc = iblock(2,jjc)
            mc = iblock(3,jjc)
            incc = jblkc*na
            istd = ist + incc 
            if (istd .gt. n) go to 55
            do 50 jjb = 3,jjlim
               jblkb = iblock(1,jjb)
               jstb = iblock(2,jjb)
               mb = iblock(3,jjb)
               incb = jblkb*na
               istdd = ist + incb
               if (istdd .gt. n) go to 50
               jdiff = jblkb - jblkc
               if (jdiff .lt. 0) go to 50
               if (jdiff .ne. 0 .and. propa) go to 40
               do 30 jjd = 1,jjlim
                  if (jdiff .eq. iblock(1,jjd)) go to 35
 30            continue
               go to 40
 35            jblkd = iblock(1,jjd)
               jstd = iblock(2,jjd)
               md = iblock(3,jjd)
               incd = jblkd*na
               if (jjd .ne. 1) call t2prod
     a               (na,na,ldt,ldt,ldt,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jt(jstd),wksp,t(ist,jstb),t(ist,jstc),
     a                t(istd,jstd))
               if (jjd .eq. 1) call t2prod
     a               (na,na,ldt,ldt,ldd,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jd,wksp,t(ist,jstb),t(ist,jstc),
     a                d(istd,1))
               if (jjd .ne. 1) call tsum
     a               (na,na,ldt,ldt,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jt(jstd),wksp,t(ist,jstb),t(ist,jstc),
     a                d(istd,1),d(istdd,1),wksp(ip1),0,omega)
               if (jjd .eq. 1) call tsum
     a               (na,na,ldt,ldt,ma,mb,mc,md,incb,
     a                incc,incd,jd,jt(jstb),jt(jstc),
     a                jd,wksp,t(ist,jstb),t(ist,jstc),
     a                d(istd,1),d(istdd,1),wksp(ip1),1,omega)
c
 40            call rowsum (ldt,na,mb,t(ist,jstb),wksp(ip1),1)
               do 42 iii = ip1,ip2
 42            wksp(iii) = omega*wksp(iii)
               call bdsol (ldd,na,na,ndt,0,d(ist,1),wksp(ip1),
     a                     wksp(ip1),0) 
               call vsubdt (ldt,1,na,na,mc,t(ist,jstc),jt(jstc),
     a                      d(istd,1),wksp(ip1),incc)
               if (jdiff .eq. 0) go to 50
               call rowsum (ldt,na,mc,t(ist,jstc),wksp(ip1),1)
               do 45 iii = ip1,ip2
 45            wksp(iii) = omega*wksp(iii)
               call bdsol (ldd,na,na,ndt,0,d(ist,1),wksp(ip1),
     a                     wksp(ip1),0) 
               call vsubdt (ldt,1,na,na,mb,t(ist,jstb),jt(jstb),
     a                      d(istdd,1),wksp(ip1),incb)
 50         continue
 55      continue
 60   continue
      return
      end 
      subroutine ibfcs4 (lddd,ldtt,nn,jd,jt,d,t,kblszz,
     a                   iblock,lbhb,ipropa,omegaa,wksp,ier)
c
c ... ibfcs4 does an incomplete block factorization of the matrix
c     contained in d and t (version 2, modified). 
c     symmetric diagonal data structure, natural ordering.
c     block ic (version 2) preconditioning.
c
c ... parameters -- 
c
c         ldd      row dimension of d array
c         ldt      row dimension of t array
c         n        size of system
c         jd       integer vector giving the diagonal numbers
c                   for the diagonal block
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         kblsz    block size 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         lbhb     number of blocks per block row 
c         ipropa   property a switch
c                   = 0   matrix does not have block property a
c                   = 1   matrix has block property a
c         omega    relaxation factor between 0. and 1.
c                   = 0   no modification
c                   = 1   full modification
c         wksp     real workspace vector
c
c ... specifications for parameters
c
      integer   jd(1), jt(1), iblock(3,3)
      dimension d(lddd,2), t(ldtt,1), wksp(1)
      logical   propa
c
      n = nn
      ldd = lddd
      ldt = ldtt
      na = kblszz
      omega = omegaa
      propa = ipropa .eq. 1
      klim = n/na
      ma = iblock(3,1)
      ndt = ma - 1
c
c ... block tridiagonal case. 
c
      if (lbhb .gt. 3) go to 25
      ip1 = n + 1
      ip2 = ip1 + na
      jblkb = iblock(1,3)
      mb = iblock(3,3)
      incb = jblkb*na
      call rowsum (ldd,n,ma,d,wksp,0)
      do 20 k = 1,klim
         ist = (k - 1)*na + 1 
         isu = k*na 
         istd = ist + incb
         call bdfac (ldd,na,na,ndt,0,d(ist,1),0)
         call bdinv (ldd,na,na,ndt,0,d(ist,1),0)
         call bmul (ldd,na,ndt,d(ist,1),d(ist,2),wksp(ist),wksp(ip1)) 
         do 10 ii = ist,isu
            if (wksp(ii) .ne. 0.0) go to 10
            ier = -12
            call ershow (ier,'ibfcs4')
            return
 10      continue
         do 15 ii = ist,isu
 15      d(ii,1) = d(ii,1) + omega*(1.0 - wksp(ii-ist+ip1))/
     a                              wksp(ii)
         if (istd .gt. n) go to 20
         call t2prod (na,ldd,ldt,ldt,ldd,m