subroutine nmcalc (coef,jcoef,wfac,jwfac,icall,subq,nn,
     a                   rhs,ubar,wksp,ier)
      implicit double precision (a-h, o-z)
c
c ... nmcalc calculates the quantities
c
c       bnorm   = sqrt (rhs,rhs)
c       bnorm1  = any other norm of rhs needed for the stopping test
c       ubarnm  = sqrt (ubar,ubar)
c
c     which are needed in the stopping tests.
c
c     the stopping tests are --
c
c    (1)  (emax/emin) * sqrt ( (r ,zt)/(rhs,inv(q)*rhs) )
c    (2)  ( 1.0/emin) * sqrt ( (zt,zt)/(u,u) )
c    (3)  (emax/emin) * sqrt ( (zt,zt)/(inv(q)*rhs,inv(q)*rhs) )
c    (4)                sqrt ( (zt,zt)/(inv(q)*rhs,inv(q)*rhs) )
c    (5)                sqrt ( (r ,r )/(rhs,rhs) )
c    (6)                sqrt ( (u-ubar,u-ubar)/(ubar,ubar) )
c    (7)  (emax/emin) * sqrt ( (r,z)/(rhs,inv(ql)*rhs) )
c    (8)  ( 1.0/emin) * sqrt ( (z,z)/(u,u) )
c    (9)  (emax/emin) * sqrt ( (z,z)/(inv(ql)*rhs,inv(ql)*rhs) )
c   (10)                sqrt ( (z,z)/(inv(ql)*rhs,inv(ql)*rhs) )
c
c ... parameters -- 
c
c        icall      key for initial or secondary call
c                    = 1   initial call 
c                    = 2   later call (needed if q is changed)
c        subq       preconditioning routine
c        n          order of system
c        rhs        right hand side
c        ubar       known solution
c        wksp       workspace vector of length n
c        ier        error code
c                    =  0  no error detected
c                    = -7  q is not positive definite
c
c ... specifications for parameters
c
      dimension rhs(1), ubar(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(1)
      external subq 
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 / srelpr, keyzer, keygs
c
c *** end   -- package common 
c
c
      n = nn
      nteste = ntest
      if (ntest .gt. 6) nteste = ntest - 6
      go to (10,50,20,20,30,40), nteste 
c
c ... bnorm1: sqrt(b,q(inv)b).
c
 10   call subq (coef,jcoef,wfac,jwfac,n,rhs,wksp)
      sum = vdot (n,rhs,wksp) 
      if (sum .ge. 0.0d0) go to 15
      ier = -7
      call ershow (ier,'nmcalc')
      return
 15   bnorm1 = max ( sqrt(sum),srelpr )
      return
c
c ... bnorm1: sqrt(q(inv)b,q(inv)b).
c
 20   call subq (coef,jcoef,wfac,jwfac,n,rhs,wksp)
      sum = vdot (n,wksp,wksp)
      bnorm1 = max ( sqrt(sum),srelpr )
      return
c
c ... bnorm.
c
 30   if (icall .eq. 2) return
      sum = vdot (n,rhs,rhs)
      bnorm = max ( sqrt(sum),srelpr )
      bnorm1 = bnorm
      return
c
c ... ubarnm.
c
 40   if (icall .eq. 2) return
      sum = vdot (n,ubar,ubar)
      ubarnm = max ( sqrt(sum),srelpr )
      return
c
c ... exit.
c
 50   return
      end 
      subroutine omgchg (ssorcp,coef,jcoef,wfac,jwfac,n,p,r)
      implicit double precision (a-h, o-z)
c
c ... omgchg changes alphab and betab for a new estimate of omega.
c
c ... parameters -- 
c
c         n       order of system (= nn)
c         p       vector from acceleration algorithm
c         r       workspace vector from acceleration algorithm
c
c ... specifications for parameters
c
      dimension p(1), r(1), coef(1), jcoef(2), wfac(1), jwfac(1)
      external ssorcp
c
c *** begin -- package common 
c
      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
      logical           omgadp
      common / itcom5 / omega, alphab, betab, fff, specr, omgadp
c
c *** end   -- package common 
c
c
c ... update alphab and betab.
c
      call ssorcp (coef,jcoef,wfac,jwfac,n,p,r,pdp,pldup)
      alphab = min (alphab, (pap/pdp) - 1.0d0)
      betab  = max (betab , pldup/pdp)
      return
      end 
      subroutine out (nn,v,iswt,noutt)
      implicit double precision (a-h, o-z)
c
c     out effects printing of residual and solution
c     vectors - called from perror1
c
c ... parameters -- 
c
c          v      vector of length n
c          iswt   labelling information 
c          nout   output device number (= noutt)
c
c ... specifications for parameters
c
      dimension v(nn)
c
         n = nn
         nout = noutt
         if (n .le. 0) return 
c
         kupper = min (n, 4) 
         if (iswt .eq. 1) write (nout,10)
 10      format (//5x,'residual vector')
         if (iswt .eq. 2) write (nout,15)
 15      format (//5x,'solution vector')
         write (nout,20) (i,i=1,kupper) 
 20      format (10x,4i15)
         write (nout,25)
 25      format (10x,65('-') /)
c
         do 35 j = 1,n,4
            kupper = min (j+3,n)
            jm1 = j - 1
            write (nout,30) jm1,(v(k),k=j,kupper) 
 30         format (4x,i5,'+  ',4d15.5) 
 35      continue
c
         return
      end 
      subroutine pbneu (suba,dsolve,coef,jcoef,wfac,jwfac,
     a                  nd,wksp,nn,r,z) 
      implicit double precision (a-h, o-z)
c
c ... pbneu computes the nd-degree block neumann polynomial 
c ... approximation to the matrix inv(a).
c     if a = d - b, where d is a dense banded matrix
c     then the output vector is --
c
c             z = (i + p + p**2 + ... + p**nd)*inv(d) * r
c
c     where   p = inv(d)*b .
c
c ... parameters -- 
c
c         suba    matrix-vector multiplication routine
c         dsolve  routine for computing inv(d)*vector
c         nd      the degree of the polynomial desired
c         wksp    workspace of length 2*n
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      external  suba, dsolve
      dimension r(1), z(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(1)
c
      n = nn
      np1 = n + 1
      call dsolve (coef,jcoef,wfac,jwfac,n,r,z)
      if (nd .le. 0) return
c
      do 20 k = 1,nd
         call suba (coef,jcoef,wfac,jwfac,n,z,wksp)
         do 10 i = 1,n
 10      wksp(i) = r(i) - wksp(i)
         call dsolve (coef,jcoef,wfac,jwfac,n,wksp,wksp(np1))
         do 15 i = 1,n
 15      z(i) = z(i) + wksp(i+n)
 20   continue
      return
      end 
      subroutine pbpii (suba,dsolve,coef,jcoef,wfac,jwfac,
     a                  ainf,alpha,beta,nd,wksp,nn,r,z)
      implicit double precision (a-h, o-z)
c
c ... pbpii computes the block nd-degree least squares polynomial
c ... approximation to the matrix inv(a).  the output vector is --
c
c ...        z = inv(d)*p  (a*inv(d)) * r
c ...                    np
c
c ... parameters -- 
c
c         suba    matrix-vector multiplication routine
c         dsolve  routine to compute inv(d)*vector
c         ainf    the infinity norm of matrix inv(d)*a
c         alpha,  the least squares weighting factors
c          beta
c         nd      the degree of the polynomial desired
c         wksp    workspace of length 2*n
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      external  suba, dsolve
      dimension r(1), z(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(1)
c
c
      n = nn
      np1 = n + 1
      al = alpha
      be = beta
c
      c1 = ((al+be+2.0d0)*(al+be+3.0d0))/(ainf*(al+2.0d0)*(al+be+2.0d0))
      call dsolve (coef,jcoef,wfac,jwfac,n,r,z)
      do 10 i = 1,n 
 10   z(i) = c1*z(i)
      if (nd .le. 0) return
c
      do 15 i = 1,n 
 15   wksp(i) = r(i)
      do 35 k = 1,nd
         fk = dble (k)
         c1 = ((2.0d0*fk+al+be+2.0d0)*(2.0d0*fk+al+be+3.0d0))/
     a         (ainf*(fk+al+2.0d0)*(fk+al+be+2.0d0))
         c2 = (fk*(fk+be)*(2.0d0*fk+al+be))/
     a        ((fk+al+1.0d0)*(fk+al+be+1.0d0)*(2.0d0*fk+al+be+2.0d0))
         call suba (coef,jcoef,wfac,jwfac,n,z,wksp(np1))
         do 20 i = 1,n
 20      wksp(n+i) = r(i) - wksp(n+i)
         do 25 i = 1,n
 25      wksp(i) = wksp(i+n) + c2*wksp(i)
         call dsolve (coef,jcoef,wfac,jwfac,n,wksp,wksp(np1))
         do 30 i = 1,n
 30      z(i) = z(i) + c1*wksp(n+i)
 35   continue
      return
      end 
      subroutine pneu (suba,coef,jcoef,wfac,jwfac,d,nd,wksp,nn,r,z)
      implicit double precision (a-h, o-z)
c
c ... pneu computes the nd-degree point neumann polynomial
c ... approximation to the matrix inv(a).  the output vector is --
c ...        z = p  (a)*r
c ...             np
c
c ... parameters -- 
c
c         suba    matrix-vector multiplication routine
c         d       vector of length n giving the diagonal elements
c                  of the matrix
c         nd      the degree of the polynomial desired
c         wksp    workspace of length n 
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      external  suba
      dimension r(1), d(1), z(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(1)
c
      n = nn
      do 10 i = 1,n 
 10   z(i) = r(i)/d(i)
      if (nd .le. 0) return
c
      do 20 k = 1,nd
         call suba (coef,jcoef,wfac,jwfac,n,z,wksp)
         do 15 i = 1,n
 15      z(i) = z(i) + (r(i) - wksp(i))/d(i)
 20   continue
      return
      end 
      subroutine ppii (suba,coef,jcoef,wfac,jwfac,ainf,
     a                 alpha,beta,nd,wksp,nn,r,z) 
      implicit double precision (a-h, o-z)
c
c ... ppii computes the nd-degree least squares polynomial
c ... approximation to the matrix inv(a).  the output vector is --
c ...        z = p  (a)*r
c ...             np
c
c ... parameters -- 
c
c         suba    matrix-vector multiplication routine
c         ainf    the infinity norm of matrix a
c         alpha,  the least squares weighting factors
c          beta
c         nd      the degree of the polynomial desired
c         wksp    workspace of length 2*n
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      external  suba
      dimension r(1), z(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(1)
c
c
      n = nn
      np1 = n + 1
      al = alpha
      be = beta
c
      c1 = ((al+be+2.0d0)*(al+be+3.0d0))/(ainf*(al+2.0d0)*(al+be+2.0d0))
      do 10 i = 1,n 
 10   z(i) = c1*r(i)
      if (nd .le. 0) return
c
      do 15 i = 1,n 
 15   wksp(i) = r(i)
      do 35 k = 1,nd
         fk = dble (k)
         c1 = ((2.0d0*fk+al+be+2.0d0)*(2.0d0*fk+al+be+3.0d0))/
     a         (ainf*(fk+al+2.0d0)*(fk+al+be+2.0d0))
         c2 = (fk*(fk+be)*(2.0d0*fk+al+be))/
     a         ((fk+al+1.0d0)*(fk+al+be+1.0d0)*(2.0d0*fk+al+be+2.0d0))
         call suba (coef,jcoef,wfac,jwfac,n,z,wksp(np1))
         do 20 i = 1,n
 20      wksp(n+i) = r(i) - wksp(n+i)
         do 25 i = 1,n
 25      wksp(i) = wksp(i+n) + c2*wksp(i)
         do 30 i = 1,n
 30      z(i) = z(i) + c1*wksp(i)
 35   continue
      return
      end 
      subroutine pbs (n,t1,t2,x)
      implicit double precision (a-h, o-z)
c
c ... pbs does a penta-diagonal back substitution (i+t1+t2)*x = y
c     where t1 and t2 are the first and second super diagonals.
c
c ... parameters -- 
c
c          n      order of the system
c          t1     vector of length n-1 containing the first super-
c                  diagonal elements
c          t2     vector of length n-2 containing the second super-
c                  diagonal elements
c          x      on input, x contains y
c                 on output, x contains the solution to
c                  (i + t1 + t2)*x = y
c
c ... specifications for parameters
c
      dimension t1(1), t2(1), x(1)
c
      x(n-1) = x(n-1) - t1(n-1)*x(n)
      do 10 i = n-2,1,-1
 10   x(i) = x(i) - t1(i)*x(i+1) - t2(i)*x(i+2)
      return
      end 
      subroutine pbsm (nn,nsize,t1,t2,x)
      implicit double precision (a-h, o-z)
c
c ... pbsm does a penta-diagonal back substitution (i+t1+t2)*x = y
c     where t1 and t2 are superdiagonals of a system composed of
c     independent subsystems of size nsize.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  order of the individual subsystems
c          t1     linear array of length n-1 containing the first
c                  super-diagonal elements of the factorizations
c          t2     linear array of length n-2 containing the second
c                  super-diagonal elements of the factorizations
c          x      on input, x contains y
c                 the solution to (i + t1 + t2)*x = y
c
c ... specifications for parameters
c
      dimension t1(nsize,1), t2(nsize,1), x(nsize,1)
c
      n = nn
      nsys = n/nsize
      do 10 j = 1,nsys
 10   x(nsize-1,j) = x(nsize-1,j) - t1(nsize-1,j)*x(nsize,j)
      do 20 i = nsize-2,1,-1
         do 15 j = 1,nsys
 15      x(i,j) = x(i,j) - t1(i,j)*x(i+1,j) - t2(i,j)*x(i+2,j)
 20   continue
      return
      end 
      subroutine pfac (nn,d,t1,t2)
      implicit double precision (a-h, o-z)
c
c ... pfac computes a factorization of a single symmetric
c     pentadiagonal matrix contained in d, t1, and t2 and
c     replaces it.
c
c ... parameters -- 
c
c          n      order of system (= nn)
c          d      vector of length n containing the diagonal
c                  elements of the matrix
c          t1     vector of length n-1 containing the first 
c                  super-diagonal elements of the matrix
c          t2     vector of length n-2 containing the second
c                  super-diagonal elements of the matrix
c
c ... specifications for parameters
c
      dimension d(1), t1(1), t2(1)
c
      n = nn
      do 10 i = 1,n-2
         dii = 1.0d0/d(i)
         d(i+1) = d(i+1) - t1(i)*t1(i)*dii
         d(i+2) = d(i+2) - t2(i)*t2(i)*dii
         t1(i+1) = t1(i+1) - t1(i)*t2(i)*dii
 10   continue
      d(n) = d(n) - t1(n-1)*t1(n-1)/d(n-1)
      do 15 i = 1,n 
 15   d(i) = 1.0d0/d(i)
      do 20 i = 1,n-1
 20   t1(i) = d(i)*t1(i)
      do 25 i = 1,n-2
 25   t2(i) = d(i)*t2(i)
      return
      end 
      subroutine pfacm (nn,nsize,d,t1,t2)
      implicit double precision (a-h, o-z)
c
c ... pfacm computes factorizations of multiple independent 
c     symmetric pentadiagonal matrices contained in d, t1, and t2.
c
c ... parameters -- 
c
c          n      order of global system (= nn)
c          nsize  size of the individual subsystems
c          d      linear array of length n containing the
c                  diagonal elements of the systems
c          t1     linear array of length n-1 containing the 
c                  first super-diagonal elements of the systems
c          t2     linear array of length n-2 containing the 
c                  second super-diagonal elements of the systems
c
c ... specifications for parameters
c
      dimension d(nsize,1), t1(nsize,1), t2(nsize,1)
c
      n = nn
      nsys = n/nsize
      do 15 i = 1,nsize-2
         do 10 j = 1,nsys
            d(i+1,j) = d(i+1,j) - (t1(i,j)**2)/d(i,j)
            d(i+2,j) = d(i+2,j) - (t2(i,j)**2)/d(i,j)
            t1(i+1,j) = t1(i+1,j) - t1(i,j)*t2(i,j)/d(i,j)
 10      continue
 15   continue
      do 20 j = 1,nsys
 20   d(nsize,j) = d(nsize,j) - (t1(nsize-1,j)**2)/d(nsize-1,j)
      call vinv (n,d)
      call vexopy (n-1,t1,d,t1,3)
      call vexopy (n-2,t2,d,t2,3)
      return
      end 
      subroutine pfacn (nn,d,t1,t2,b1,b2)
      implicit double precision (a-h, o-z)
c
c ... pfacn computes a factorization of a single nonsymmetric
c     pentadiagonal matrix contained in d,t1,t2,b1, and b2
c     and replaces it.
c
c ... parameters -- 
c
c          n      order of system (= nn)
c          d      vector of length n containing the diagonal
c                  elements of the matrix
c          t1     vector of length n-1 containing the first 
c                  super-diagonal elements of the matrix
c          t2     vector of length n-2 containing the second
c                  super-diagonal elements of the matrix
c          b1     vector of length n-1 containing the first 
c                  sub-diagonal elements of the matrix
c          b2     vector of length n-2 containing the second
c                  sub-diagonal elements of the matrix
c
c ... specifications for parameters
c
      dimension d(1), t1(1), t2(1), b1(1), b2(1)
c
      n = nn
      do 10 i = 1,n-2
         dii = 1.0d0/d(i)
         d(i+1) = d(i+1) - b1(i)*t1(i)*dii
         d(i+2) = d(i+2) - b2(i)*t2(i)*dii
         t1(i+1) = t1(i+1) - b1(i)*t2(i)*dii
         b1(i+1) = b1(i+1) - b2(i)*t1(i)*dii
 10   continue
      d(n) = d(n) - b1(n-1)*t1(n-1)/d(n-1)
      do 15 i = 1,n 
 15   d(i) = 1.0d0/d(i)
      do 20 i = 1,n-1
         t1(i) = d(i)*t1(i)
         b1(i) = d(i)*b1(i)
 20   continue
      do 25 i = 1,n-2
         t2(i) = d(i)*t2(i)
         b2(i) = d(i)*b2(i)
 25   continue
      return
      end 
      subroutine pfacnm (nn,nsize,d,t1,t2,b1,b2)
      implicit double precision (a-h, o-z)
c
c ... pfacnm computes factorizations of multiple independent
c     nonsymmetric pentadiagonal matrices contained in
c     d,t1,t2,b1, and b2.
c
c ... parameters -- 
c
c          n      order of global system (= nn)
c          nsize  order of single subsystem
c          d      linear array of length n containing the
c                  diagonal elements of the systems
c          t1     linear array of length n-1 containing the first
c                  super-diagonal elements of the systems
c          t2     linear array of length n-2 containing the second
c                  super-diagonal elements of the systems
c          b1     linear array of length n-1 containing the first
c                  sub-diagonal elements of the systems
c          b2     linear array of length n-2 containing the second
c                  sub-diagonal elements of the systems
c
c ... specifications for parameters
c
      dimension d(nsize,1), t1(nsize,1), b1(nsize,1), t2(nsize,1),
     a            b2(nsize,1) 
c
      n = nn
      nsys = n/nsize
      do 15 i = 1,nsize-2
         do 10 j = 1,nsys
            d(i+1,j) = d(i+1,j) - b1(i,j)*t1(i,j)/d(i,j)
            d(i+2,j) = d(i+2,j) - b2(i,j)*t2(i,j)/d(i,j)
            t1(i+1,j) = t1(i+1,j) - b1(i,j)*t2(i,j)/d(i,j)
            b1(i+1,j) = b1(i+1,j) - b2(i,j)*t1(i,j)/d(i,j)
 10      continue
 15   continue
      do 20 j = 1,nsys
 20   d(nsize,j) = d(nsize,j) - b1(nsize-1,j)*t1(nsize-1,j)/
     a               d(nsize-1,j)
      call vinv (n,d)
      call vexopy (n-1,t1,d,t1,3)
      call vexopy (n-2,t2,d,t2,3)
      call vexopy (n-1,b1,d,b1,3)
      call vexopy (n-2,b2,d,b2,3)
      return
      end 
      subroutine pfs (n,b1,b2,x)
      implicit double precision (a-h, o-z)
c
c ... pfs does a penta-diagonal forward substitution  (i+b1+b2)*x = y 
c     where b1 and b2 are the first and second sub-diagonals.
c
c ... parameters -- 
c
c          n      order of system
c          b1     vector of length n-1 containing the first 
c                  sub-diagonal elements
c          b2     vector of length n-2 containing the second
c                  sub-diagonal elements
c          x      on input, x contains y
c                 on output, x contains the solution to
c                  (i + b1 + b2)*x = y
c
c ... specifications for parameters
c
      dimension b1(1), b2(1), x(2)
c
      x(2) = x(2) - b1(1)*x(1)
      do 10 i = 3,n 
 10   x(i) = x(i) - b1(i-1)*x(i-1) - b2(i-2)*x(i-2)
      return
      end 
      subroutine pfsm (nn,nsize,b1,b2,x)
      implicit double precision (a-h, o-z)
c
c ... pfsm does a penta-diagonal forward substitution (i+b1+b2)*x = y 
c     where b1 and b2 are subdiagonals of a system composed of
c     independent subsystems of size nsize.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  order of the individual subsystems
c          b1     linear array of length n-1 containing the first
c                  sub-diagonal elements of the factorizations
c          b2     linear array of length n-2 containing the second
c                  sub-diagonal elements of the factorizations
c          x      on input, x contains y
c                 on output, x contains 
c                 the solution to (i + b1 + b2)*x = y
c
c ... specifications for parameters
c
      dimension b1(nsize,1), b2(nsize,1), x(nsize,1)
c
      n = nn
      nsys = n/nsize
      do 10 j = 1,nsys
 10   x(2,j) = x(2,j) - b1(1,j)*x(1,j)
      do 20 i = 3,nsize
         do 15 j = 1,nsys
 15      x(i,j) = x(i,j) - b1(i-1,j)*x(i-1,j) - b2(i-2,j)*x(i-2,j)
 20   continue
      return
      end 
      subroutine psoln (nn,d,t1,t2,b1,b2,y,x)
      implicit double precision (a-h, o-z)
c
c ... psoln solves the system ax = y for x, where a is a single
c     pentadiagonal system.  d, t1, t2, b1, and b2 contain
c     the main, first and second super, and first and second sub
c     diagonals, respectively, of the factorization.
c
c ... parameters -- 
c
c          n      order of system
c          d      vector of length n containing the diagonal
c                  elements of the factorization matrix
c          t1     vector of length n-1 containing the first 
c                  super-diagonal elements of the factorization
c          t2     vector of length n-2 containing the second
c                  super-diagonal elements of the factorization
c          b1     vector of length n-1 containing the first 
c                  sub-diagonal elements of the factorization
c          b2     vector of length n-2 containing the second
c                  sub-diagonal elements of the factorization
c          y      the right-hand side
c          x      the solution to ax = y
c
c ... specifications for parameters
c
      dimension d(1), t1(1), t2(1), b1(1), b2(1), x(1), y(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      call pfs (n,b1,b2,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call pbs (n,t1,t2,x)
      return
      end 
      subroutine psolnm (nn,nsize,d,t1,t2,b1,b2,y,x)
      implicit double precision (a-h, o-z)
c
c ... psolnm solves the system ax = y for x, where a contains
c     multiple pentadiagonal systems.  d, t1, t2, b1, and b2 are
c     the main, first and second super, and the first and second
c     sub diagonals, respectively, of the factorization.
c
c ... parameters -- 
c
c          n      order of system
c          nsize  size of an individual subsystem 
c          d      vector of length n containing the diagonal
c                  elements of the factorization matrix
c          t1     vector of length n-1 containing the first 
c                  super-diagonal elements of the factorization
c          t2     vector of length n-2 containing the second
c                  super-diagonal elements of the factorization
c          b1     vector of length n-1 containing the first 
c                  sub-diagonal elements of the factorization
c          b2     vector of length n-2 containing the second
c                  sub-diagonal elements of the factorization
c          y      the right-hand side
c          x      the solution to ax = y
c
c ... specifications for parameters
c
      dimension d(1), t1(1), t2(1), b1(1), b2(1), x(1), y(1)
c
      n = nn
      do 10 i = 1,n 
 10   x(i) = y(i)
      call pfsm (n,nsize,b1,b2,x)
      do 15 i = 1,n 
 15   x(i) = d(i)*x(i)
      call pbsm (n,nsize,t1,t2,x)
      return
      end 
      subroutine parsi
      implicit double precision (a-h, o-z)
c
c ... parsi computes the iteration 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 / srelpr, keyzer, keygs
c
c *** end   -- package common 
c
      rhol = rho
      if (is - 1) 10,15,20
 10   rho = 1.0d0
      go to 25
 15   rho = 1.0d0/(1.0d0 - sigma*sigma/2.0d0) 
      go to 25
 20   rho = 1.0d0/(1.0d0 - sigma*sigma*rhol/4.0d0)
c
c ... compute alpha, beta.
c
 25   alpha = rho*gamma
      beta = rhol*(rho - 1.0d0)/rho
      return
      end 
      subroutine permas (isym,nn,nzz,ia,ja,a,wksp,p)
      implicit double precision (a-h, o-z)
c
c ... permas takes the sparse matrix representation of the
c     matrix and permutes both rows and columns, overwriting
c     the previous structure.  (sparse data structure)
c
c ... parameters -- 
c
c          isym      switch for symmetric storage 
c                     = 0   matrix is symmetric
c                     = 1   matrix is nonsymmetric
c          n         size of system
c          nz        length of ia, ja, and a vectors
c          ia        vector of i values 
c          ja        vector of j values 
c          a         vector of matrix coefficients
c          wksp      workspace vector of length n 
c          p         permutation vector 
c
c ... it is assumed that the i-th entry of the permutation vector
c     p indicates the row the i-th row gets mapped into.  (i.e.
c     if ( p(i) = j ) row i gets mapped into row j)
c
c ... specifications for parameters
c
      dimension a(1), wksp(1) 
      integer   ia(1), ja(1), p(1)
c
      n = nn
      nz = nzz
c
c ... explicit gathers.
c
      call vgathi (nz,p,ia,ia)
      call vgathi (nz,p,ja,ja)
      do 5 i = 1,n
 5    wksp(i) = a(i)
      call vscatr (n,wksp,p,a)
      do 10 i = 1,n 
         ia(i) = i
         ja(i) = i
 10   continue
c
c ... convert to upper triangular elements for symmetric storage
c
      if (isym .eq. 1) return 
      np1 = n + 1
      do 15 i = np1,nz
         if (ia(i) .le. ja(i)) go to 15 
         idum = ia(i)
         ia(i) = ja(i)
         ja(i) = idum
 15   continue
      return
      end 
      subroutine permat (ndim,maxnz,coef,jcoef,wksp,iwksp,nn,p)
      implicit double precision (a-h, o-z)
c
c ... permat takes the sparse matrix representation of the
c     matrix and permutes both rows and columns, overwriting
c     the previous structure.  (purdue data structure)
c
c ... parameters -- 
c
c          ndim      row dimension of coef array in defining routine
c          maxnz     number of columns in coef and jcoef arrays
c          coef      array of matrix coefficients 
c          jcoef     array of matrix columns numbers
c          wksp      workspace array of length n
c          iwksp     integer workspace array of length n
c          n         order of system (= nn)
c          p         permutation vector 
c
c ... it is assumed that the i-th entry of the permutation vector
c     p indicates the row the i-th row gets mapped into.  (i.e.
c     if ( p(i) = j ) row i gets mapped into row j)
c
c ... specifications for parameters
c
      dimension coef(ndim,1), wksp(1)
      integer   jcoef(ndim,1), iwksp(1), p(1)
c
      n = nn
      if (n .le. 0) return
      do 20 j = 1,maxnz
         do 10 i = 1,n
            wksp(i) = coef(i,j)
            iwksp(i) = jcoef(i,j)
 10      continue
         call vscatr (n,wksp,p,coef(1,j))
         call vscati (n,iwksp,p,jcoef(1,j))
         call vgathi (n,p,jcoef(1,j),jcoef(1,j))
 20   continue
      return
      end 
      subroutine perror1 (suba,coef,jcoef,wfac,jwfac,nn,u,rhs,
     a                   wksp,digtt1,digtt2,idgtts)
      implicit double precision (a-h, o-z)
c
c     perror1 computes the residual, r = rhs - a*u.  the user
c     also has the option of printing the residual and/or the
c     unknown vector depending on idgts.
c
c ... parameters -- 
c
c          suba   matrix-vector multiplication routine
c          n      dimension of matrix (= nn)
c          u      latest estimate of solution
c          rhs    right hand side of matrix problem
c          wksp   workspace vector of length n
c          digit1 output - measure of accuracy of stopping test (= digtt1
c          digit2 output - measure of accuracy of solution (= digtt2) 
c          idgts   parameter controlling level of output (= idgtts)
c                    if idgts < 1 or idgts > 4, then no output.
c                            = 1, then number of digits is printed, pro-
c                                 vided level .ge. 1
c                            = 2, then solution vector is printed, pro-
c                                 vided level .ge. 1
c                            = 3, then residual vector is printed, pro-
c                                 vided level .ge. 1
c                            = 4, then both vectors are printed, pro- 
c                                 vided level .ge. 1
c
c ... specifications for parameters
c
      external suba 
      dimension rhs(1), u(1), wksp(1), coef(1), jcoef(2),
     a          wfac(1), jwfac(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 / srelpr, keyzer, keygs
c
c *** end   -- package common 
c
      n = nn
      idgts = idgtts
      digit1 = 0.0d0
      digit2 = 0.0d0
c
      digit1 = -dlog10 (abs (srelpr))
      if (stptst .gt. 0.0d0) digit1 = -dlog10 (abs (stptst))
         call suba (coef,jcoef,wfac,jwfac,n,u,wksp)
         call vexopy (n,wksp,rhs,wksp,2)
         rnrm = sqrt ( vdot (n,wksp,wksp) )
         sum = vdot (n,rhs,rhs)
         bnorm = max ( sqrt(sum),srelpr )
         temp = rnrm/bnorm
         if (temp .eq. 0.0d0) go to 10
         digit2 = -dlog10 (abs (temp))
         go to 15
c
 10   digit2 = -dlog10 (abs (srelpr))
c
 15      if ((idgts .lt. 1) .or. (level .le. 0)) go to 25
         write (nout,20) digit1,digit2
 20      format (/10x,'approx. no. of digits in stopping test =',
     a          f5.1,2x,'(digit1)'
     b          /10x,'approx. no. of digits in ratio test    =',
     c          f5.1,2x,'(digit2)')
c
         if (idgts .le. 1 .or. idgts .gt. 4) go to 25
         if (idgts .ge. 3) call out (n,wksp,1,nout)
         if (idgts .ne. 3) call out (n,u,2,nout)
c
 25   continue
      digtt1 = digit1
      digtt2 = digit2
      return
      end 
      subroutine pervec (nn,p,v,wksp)
      implicit double precision (a-h, o-z)
c
c ... pervec permutes a vector as dictated by the permutation
c ... vector p.  if p(i) = j, then v(j) gets v(i).
c
c ... parameters -- 
c
c          n       length of vectors p, v, and wksp  (= nn) 
c          p       integer permutation vector
c          v       vector to be permuted
c          wksp    workspace vector of length n
c
c ... specifications for parameters
c
      integer   p(1)
      dimension v(1), wksp(1) 
c
      n = nn
      if (n .le. 0) return
      do 10 i = 1,n 
         wksp(i) = v(i)
 10   continue
      call vscatr (n,wksp,p,v)
      return
      end 
      subroutine pgen (nn,p,ip,nc,ncolor)
      implicit double precision (a-h, o-z)
c
c ... pgen constructs the permutation vector p and its inverse
c ... ip for a coloring given by p.
c
c ... parameters -- 
c
c         n         order of system (= nn)
c         p         vector from prbndx upon input 
c                   permutation vector upon output
c         ip        integer workspace vector upon input
c                   inverse permutation vector upon output
c         nc        number of points for each color (output)
c         ncolor    number of colors
c
c ... specifications for parameters
c
      integer p(1), ip(1), nc(1)
c
      n = nn
c
c ... determine number of colors and number of elements for each
c     color.
c
      ncolor = 0
      do 5 i = 1,n
 5    nc(i) = 0
      do 10 i = 1,n 
         ic = p(i)
         if (ncolor .lt. ic) ncolor = ic
         nc(ic) = nc(ic) + 1
 10   continue
c
c ... construct permutation vector.
c
      ip(1) = 1
      do 15 i = 2,ncolor
         ip(i) = ip(i-1) + nc(i-1)
 15   continue
      do 20 i = 1,n 
         ic = p(i)
         p(i) = ip(ic)
         ip(ic) = ip(ic) + 1
 20   continue
c
c ... construct inverse permutation vector.
c
      do 25 i = 1,n 
         j = p(i)
         ip(j) = i
 25   continue
      return
      end 
      subroutine pjac (diag,nn,r,z)
      implicit double precision (a-h, o-z)
c
c ... pjac does the point jacobi preconditioning. 
c
c ... parameters -- 
c
c         diag    vector of length n containing the diagonal
c                  elements of the coefficient matrix
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      dimension r(1), z(1), diag(1)
c
      n = nn
      do 10 i = 1,n 
 10   z(i) = r(i)/diag(i)
      return
      end 
      subroutine pmdg (ndim,mdim,nn,maxnz,jcoef,coef,ncol,nc,p,ip,
     a                 maxd,maxnew,jcnew,wksp,iwksp,isym,ier)
      implicit double precision (a-h, o-z)
c
c ... pmdg permutes the matrix according to index vector p, 
c     and, if room allows, stores the permuted matrix in
c     diagonal format.  there will be enough room if the number
c     of diagonals needed does not exceed mdim.
c
c ... parameters -- 
c
c        ndim      row dimension of coef and jcoef arrays
c                   in defining routine 
c        mdim      column dimension of coef and jcoef arrays in
c                   defining routine
c        n         order of system (active row size of coef and jcoef)
c        maxnz     active column size of coef and jcoef
c        jcoef     integer array of column numbers
c        coef      floating point array of coefficients
c        ncolor    number of colors in the permutation (= ncol)
c        nc        integer vector of length ncolor giving the
c                   number of nodes for each color
c        p         permutation vector
c        ip        inverse permuation vector
c        maxd      active columns in permuted matrix
c        maxnew    integer vector giving the number of diagonals
c                   created for each color
c        jcnew     integer array of size ncolor*max(maxnew(i))
c                   giving the diagonal numbers for each color
c        wksp      floating point workspace of length n
c        iwksp     integer workspace of length 2*n
c        isym      symmetric storage switch
c                   = 0    symmetric storage
c                   = 1    nonsymmetric storage
c        ier       error flag 
c                  =  0   no errors detected
c                  = -9   mdim is less than the number of columns
c                          needed in coef to store the permuted
c                          matrix in diagonal format
c
c ... specifications for parameters
c
      integer jcoef(2), nc(1), p(1), maxnew(1), jcnew(ncol,1),
     a        iwksp(1), ip(1) 
      dimension coef(ndim,1), wksp(1)
c
c
      n = nn
      ncolor = ncol 
c
c ... fill out rest of matrix if symmetric storage is used. 
c
      if (isym .ne. 0) go to 2
      maxd = 2*maxnz - 1
      if (mdim .lt. maxd) ier = -9
      if (ier .lt. 0) return
c
      do 1 j = 2,maxnz
         ind = jcoef(j)
         len = n - ind
         jcol = maxnz + j - 1 
         jcoef(jcol) = -ind
         call vfill (ind,coef(1,jcol),0.0d0)
         call vcopy (len,coef(1,j),coef(ind+1,jcol))
 1    continue
      maxnz = maxd
c
c ... determine the number of created diagonals.
c
 2    do 5 i = 1,ncolor
         maxnew(i) = 1
         jcnew(i,1) = 0
 5    continue
      do 35 j = 2,maxnz
         ind = jcoef(j)
         do 10 i = 1,n
            iwksp(n+i) = i + ind
            if (coef(i,j) .eq. 0.0d0) iwksp(n+i) = i
 10      continue
         call vscati (n,iwksp(n+1),p,iwksp)
         call vgathi (n,p,iwksp,iwksp)
         do 15 i = 1,n
 15      iwksp(i) = iwksp(i) - i
         ist = 1
         do 30 k = 1,ncolor
            ncc = nc(k)
            ied = ist + ncc - 1
            lim = maxnew(k)
            do 25 i = ist,ied 
               id = iwksp(i)
               do 20 jj = 1,lim
                  if (jcnew(k,jj) .eq. id) go to 25
 20            continue
               lim = lim + 1
               maxnew(k) = lim
               if (lim .gt. mdim) go to 40
               jcnew(k,lim) = id
 25         continue
            ist = ist + ncc
 30      continue
 35   continue
c
c ... determine maxd.
c
 40   maxd = -1
      do 45 k = 1,ncolor
 45   maxd = max (maxd,maxnew(k))
      if (mdim .lt. maxd) ier = -9
      if (ier .lt. 0) return
c
c ... permute matrix.
c
      do 55 j = 1,maxnz
         do 50 i = 1,n
 50      wksp(i) = coef(i,j)
         call vscatr (n,wksp,p,coef(1,j))
 55   continue
c
c ... rearrange rows.
c
      ist = 1
      do 85 k = 1,ncolor
         ncc = nc(k)
         ied = ist + ncc - 1
         lim = maxnew(k)
         do 62 l = 1,lim
            jcol = jcnew(k,l) 
            iwksp(n+jcol) = l 
 62      continue
         do 80 i = ist,ied
            iip = ip(i)
            do 60 j = 2,maxnz 
 60         wksp(j) = coef(i,j)
            do 63 j = 2,maxd
 63         coef(i,j) = 0.0d0
            do 75 j = 2,maxnz 
               if (wksp(j) .eq. 0.0d0) go to 75
               icol = p(iip + jcoef(j)) - i
               l = iwksp(n+icol)
               coef(i,l) = wksp(j)
 75         continue
 80      continue
         ist = ist + ncc
 85   continue
      return
      end 
      subroutine prbndx (nn,ndim,maxnzz,jcoef,coef,p,ip,propa,nstore) 
      implicit double precision (a-h, o-z)
c
c**************************************************************
c
c     (purdue, diagonal data structures)
c     prbndx determines if the matrix has property a.
c     this algorithm assumes all neighbors of a particular node
c     are known.
c
c     the algorithm is to mark the first node as red (arbitrary).
c     all of its adjacent nodes are marked black and placed in
c     a stack.  the remainder of the code pulls the first node
c     off the top of the stack and tries to type its adjacent nodes.
c     the typing of the adjacent point is a five way case statement
c     which is well commented below (see do loop 50).
c
c     the array p is used both to keep track of the color of a node
c     (red node is positive, black is negative) but also the father
c     node that caused the color marking of that point.  since
c     complete information on the adjacency structure is hard to come 
c     by, this forms a link to enable the color change of a partial
c     tree when a recoverable color conflict occurs.
c
c     the array ip is used as a stack to point to the set of nodes
c     left to be typed that are known to be adjacent to the current
c     father node.
c
c
c*********************************************************************
c
c ... input parameters --
c
c        n      number of nodes.  (integer, scalar) (= nn)
c        ndim   row dimension of coef array
c        maxnz  maximum number of nonzeros per row
c        jcoef  integer data array
c        coef   floating point data array
c        p,ip   integer workspace vectors of length n
c        nstore data structure switch
c                = 1  purdue
c                = 2  diagonal (symmetric or nonsymmetric)
c
c ... output parameters --
c
c        p      contains information for constructing the permutation 
c               array upon output
c        propa  a logical variable which is set to .true. if the
c               matrix has property a and .false. otherwise 
c
c
c******************************************************************** 
c
c ... specifications for parameters
c
      integer   p(1), ip(1), jcoef(ndim,1)
      dimension coef(ndim,1)
      logical   propa
c
c
c ... specifications for local variables
c
      integer   first, old, young, curtyp, type
c
c-----------------------------------------------------------------------
c
         n = nn
         maxnz = maxnzz
         do 5 i = 1,n
            p(i) = 0
            ip(i) = 0
 5       continue
c
c ... handle the first set of points until some adjacent points
c ... are found
c
         first = 1
c
 10      p(first) = first
         if (maxnz .gt. 1) go to 20
c
c ... search for next entry that has not been marked
c
         if (first .eq. n) go to 65
         ibgn = first + 1
         do 15 i = ibgn,n
            if (p(i) .ne. 0) go to 15
            first = i
            go to 10
 15      continue
         go to 65
c
c ... first set of adjacent points found
c
 20      next = 1
         last = 1
         ip(1) = first
c
c ... loop over labeled points indicated in the stack stored in
c ... the array ip
c
 25      k = ip(next)
         curtyp = p(k)
         nxttyp = -curtyp
         if (maxnz .le. 0) go to 55
         do 50 j = 1,maxnz
            if (nstore .eq. 1) jcol = jcoef(k,j)
            if (nstore .ge. 2) jcol = k + jcoef(j,1)
c
c ... determine if element (k,j) is a diagonal element or zero.
c
            if (jcol .lt. 1 .or. jcol .gt. n .or. jcol .eq. k)
     a                                    go to 50
            if (coef(k,j) .eq. 0.0d0) go to 50
c
            type = p(jcol)
c
c==================================================================
c
c     the following is a five way case statement dealing with the
c     labeling of the adjacent node.
c
c ... case i.  if the adjacent node has already been labeled with
c              label equal to nxttyp, then skip to the next adjacent
c              node.
c
         if (type .eq. nxttyp) go to 50 
c
c ... case ii.  if the adjacent node has not been labeled yet label
c               it with nxttyp and enter it in the stack
c
         if (type .ne. 0) go to 30
            last = last + 1
            ip(last) = jcol
            p(jcol) = nxttyp
            go to 50
c
c ... case iii.  if the adjacent node has already been labeled with
c                opposite color and the same father seed, then there
c                is an irrecoverable color conflict.
c
 30      if (type .eq. curtyp) go to 999
c
c ... case iv.  if the adjacent node has the right color and a different
c               father node, then change all nodes of the youngest fathe
c               node to point to the oldest father seed and retain the
c               same colors.
c
         if (type * nxttyp .lt. 1) go to  40
            old   = min ( iabs(type), iabs(nxttyp) )
            young = max ( iabs(type), iabs(nxttyp) )
            do  35 l = young,n
               if (iabs(p(l)) .eq. young) p(l) = isign(old, p(l))
 35         continue
            curtyp = p(k)
            nxttyp = -curtyp
            go to 50
c
c ... case v.  if the adjacent node has the wrong color and a different
c              father node, then change all nodes of the youngest father
c              node to point to the oldest father node along with
c              changing their colors.  since until this time the
c              youngest father node tree has been independent no other
c              color conflicts will arise from this change. 
c
 40      old   = min ( iabs(type), iabs(nxttyp) )
         young = max ( iabs(type), iabs(nxttyp) )
         do  45 l = young,n
            if (iabs(p(l)) .eq. young) p(l) = isign(old, -p(l))
 45      continue
         curtyp = p(k)
         nxttyp = -curtyp
c
c
c ... end of case statement
c
c==================================================================
 50      continue
c
c ... advance to next node in the stack 
c
 55      next = next + 1
         if (next .le. last) go to 25
c
c ... all nodes in the stack have been removed
c
c ... check for nodes not labeled.  if any are found
c ... start the labeling process again at the first
c ... node found that is not labeled.
c
         ibgn = first + 1
         do 60 i = ibgn,n
            if (p(i) .ne. 0) go to 60
               first = i
               go to 10
 60      continue
c
c
c===================================================================
c
c
c ... all nodes are now typed either red or black.
c ... red-black ordering possible.
c
 65      propa = .true.
         do 70 i = 1,n
            if (p(i) .ge. 0) p(i) = 1
            if (p(i) .le. 0) p(i) = 2
 70      continue
         return
c
c ...... type conflict
c
 999  propa = .false.
      return
      end 
      subroutine prbblk (ncol,ndis,iblock,lbhb,p,ip,propa)
      implicit double precision (a-h, o-z)
c
c**************************************************************
c
c     (block structure)
c     prbblk determines if the matrix has block property a. 
c     see routine prbndx for an explanation of the algorithm
c
c**************************************************************
c
c ... input parameters --
c
c         ncolor   number of diagonal blocks
c         ndis     number of distinct diagonal blocks
c         iblock   integer array of size 3 by ndis by max(lbhb(i))
c                   giving block constants
c         lbhb     integer vector of size ndis giving the number
c                   of diagonal blocks for each distinct block size.
c         p,ip     integer workspace vectors of length ncolor
c
c ... output parameters --
c
c        p      contains information for constructing the permutation 
c               array upon output
c        propa  a logical variable which is set to .true. if the
c               matrix has block property a and .false. otherwise
c
c
c******************************************************************** 
c
c ... specifications for parameters
c
      integer   p(1), ip(1), iblock(3,ndis,1), lbhb(1)
      logical   propa
c
c ... specifications for local variables
c
      integer   first, old, young, curtyp, type
c
c
         ncolor = ncol
         ndist = ndis
         index = 1
         do 5 i = 1,ncolor
            p(i) = 0
            ip(i) = 0
 5       continue
c
c ... handle the first set of points until some adjacent points
c ... are found
c
         first = 1
c
 10      p(first) = first
         if (ndist .gt. 1) index = first
         maxnz = lbhb(index)
         if (maxnz .gt. 1) go to 20
c
c ... search for next entry that has not been marked
c
         if (first .eq. ncolor) go to 65
         do 15 i = first+1,ncolor
            if (p(i) .ne. 0) go to 15
            first = i
            go to 10
 15      continue
         go to 65
c
c ... first set of adjacent points found
c
 20      next = 1
         last = 1
         ip(1) = first
c
c ... loop over labeled points indicated in the stack stored in
c ... the array ip
c
 25      k = ip(next)
         curtyp = p(k)
         nxttyp = -curtyp
         if (ndist .gt. 1) index = k
         maxnz = lbhb(index)
         if (maxnz .le. 0) go to 55
         do 50 j = 1,maxnz
            jcol = k + iblock(1,index,j)
c
c ... determine if element (k,j) is a diagonal element or zero.
c
            if (jcol .lt. 1 .or. jcol .gt. ncolor .or.
     a          jcol .eq. k) go to 50
            if (iblock(3,index,j) .eq. 0) go to 50
c
            type = p(jcol)
c
c==================================================================
c
c     the following is a five way case statement dealing with the
c     labeling of the adjacent node.
c
c ... case i.  if the adjacent node has already been labeled with
c              label equal to nxttyp, then skip to the next adjacent
c              node.
c
         if (type .eq. nxttyp) go to 50 
c
c ... case ii.  if the adjacent node has not been labeled yet label
c               it with nxttyp and enter it in the stack
c
         if (type .ne. 0) go to 30
            last = last + 1
            ip(last) = jcol
            p(jcol) = nxttyp
            go to 50
c
c ... case iii.  if the adjacent node has already been labeled with
c                opposite color and the same father seed, then there
c                is an irrecoverable color conflict.
c
 30      if (type .eq. curtyp) go to 999
c
c ... case iv.  if the adjacent node has the right color and a different
c               father node, then change all nodes of the youngest fathe
c               node to point to the oldest father seed and retain the
c               same colors.
c
         if (type * nxttyp .lt. 1) go to  40
            old   = min ( iabs(type), iabs(nxttyp) )
            young = max ( iabs(type), iabs(nxttyp) )
            do  35 l = young,ncolor
               if (iabs(p(l)) .eq. young) p(l) = isign(old, p(l))
 35         continue
            curtyp = p(k)
            nxttyp = -curtyp
            go to 50
c
c ... case v.  if the adjacent node has the wrong color and a different
c              father node, then change all nodes of the youngest father
c              node to point to the oldest father node along with
c              changing their colors.  since until this time the
c              youngest father node tree has been independent no other
c              color conflicts will arise from this change. 
c
 40      old   = min ( iabs(type), iabs(nxttyp) )
         young = max ( iabs(type), iabs(nxttyp) )
         do  45 l = young,ncolor
            if (iabs(p(l)) .eq. young) p(l) = isign(old, -p(l))
 45      continue
         curtyp = p(k)
         nxttyp = -curtyp
c
c
c ... end of case statement
c
c==================================================================
 50      continue
c
c ... advance to next node in the stack 
c
 55      next = next + 1
         if (next .le. last) go to 25
c
c ... all nodes in the stack have been removed
c
c ... check for nodes not labeled.  if any are found
c ... start the labeling process again at the first
c ... node found that is not labeled.
c
         do 60 i = first+1,ncolor
            if (p(i) .ne. 0) go to 60
               first = i
               go to 10
 60      continue
c
c
c===================================================================
c
c
c ... all nodes are now typed either red or black.
c ... red-black ordering possible.
c
 65      propa = .true.
         do 70 i = 1,ncolor
            if (p(i) .ge. 0) p(i) = 1
            if (p(i) .le. 0) p(i) = 2
 70      continue
         return
c
c ...... type conflict
c
 999  propa = .false.
      return
      end 
      subroutine prep1 (nn,ndim,maxnzz,jcoef,coef,ier)
      implicit double precision (a-h, o-z)
c
c ... prep1 puts the diagonal elements of the matrix in column one
c     of coef  (purdue data structure)
c
c ... parameters -- 
c
c         n       dimension of matrix ( = 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         ier     error flag -- on return, values mean
c                      0 -- no errors detected
c                     -5 -- nonexistent diagonal element
c
c ... specifications for parameters
c
      integer   jcoef(ndim,1) 
      dimension coef(ndim,1)
c
      n = nn
      maxnz = maxnzz
c
      do 20 i = 1,n 
         do 10 j = 1,maxnz
            if (jcoef(i,j) .eq. i) go to 15
 10      continue
c
c ... no diagonal entry for row i.
c
         ier = -5
         return
c
c ... switch entries so that diagonal element is in column 1.
c
 15      if (j .eq. 1) go to 20
         save = coef(i,j)
         coef(i,j) = coef(i,1)
         jcoef(i,j) = jcoef(i,1)
         coef(i,1) = save
         jcoef(i,1) = i
 20   continue
      return
      end 
      subroutine prep2 (nn,ndim,maxnzz,jcoef,coef,wksp,ier) 
      implicit double precision (a-h, o-z)
c
c ... prep2 puts the diagonal entries of the matrix into column
c     one of coef.  (diagonal data structure)
c
c ... parameters -- 
c
c         n       dimension of matrix ( = 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         wksp    workspace array of size n
c         ier     error flag -- on return, values mean
c                      0 -- no errors detected
c                     -5 -- nonexistent diagonal element
c
c ... specifications for parameters
c
      integer   jcoef(2)
      dimension coef(ndim,1), wksp(1)
c
      n = nn
      maxnz = maxnzz
c
      do 10 j = 1,maxnz
         if (jcoef(j) .eq. 0) go to 15
 10   continue
c
c ... no main diagonal.
c
      ier = -5
      return
c
c ... switch diagonals so that main diagonal is in column 1.
c
 15   if (j .eq. 1) return
      do 20 i = 1,n 
         wksp(i) = coef(i,1)
         coef(i,1) = coef(i,j)
         coef(i,j) = wksp(i)
 20   continue
      jcoef(j) = jcoef(1)
      jcoef(1) = 0
      return
      end 
      subroutine prep3 (n,nz,ia,ja,a,m,np,iwksp)
      implicit double precision (a-h, o-z)
c
c ... prep3 puts the diagonal elements of the matrix in the 
c     first n locations of the data structure, adds duplicate
c     triples, and defines the partition for matrix-vector
c     products.
c
c ... parameters -- 
c
c         n       number of equations
c         nz      length of ia, ja, and a vectors 
c         ia      vector of i values
c         ja      vector of j values
c         a       vector of matrix coefficients
c         m       number of partitions (output)
c         np      on output, np contains the partition pointers.
c                  it must be at least m+1 in length.
c         iwksp   integer workspace vector of length n
c
c ... specifications for parameters
c
      integer ia(1), ja(1), iwksp(1), np(1)
      dimension a(1)
c
c ... eliminate duplicates from the vectors by adding their 
c     values in the a vector.  first, sort the vectors by
c     rows first and then by columns within each row.
c
      call vsrta1 (nz,ia,ja,a)
c
c ... add duplicates.
c
      l = 1
      do 10 k = 2,nz
         i = ia(k)
         j = ja(k)
         aval = a(k)
         if (i .eq. ia(l) .and. j .eq. ja(l)) go to 5
         l = l + 1
         ia(l) = i
         ja(l) = j
         a(l) = aval
         go to 10
 5       a(l) = a(l) + aval
 10   continue
      nz = l
c
c ... put main diagonal elements first. 
c
      do 20 k = 1,nz
 15      i = ia(k)
         j = ja(k)
         if (i .ne. j) go to 20
         if (i .eq. k) go to 20
         val = a(k) 
         ia(k) = ia(i)
         ja(k) = ja(i)
         a(k) = a(i)
         ia(i) = i
         ja(i) = i
         a(i) = val 
         go to 15
 20   continue
c
c ... define partitions.
c
      kbgn = n + 1
      krep = kbgn
      mm = 1
      np(1) = 1
 25   mm = mm + 1
      np(mm) = kbgn 
      do 30 i = 1,n 
 30   iwksp(i) = 0
      nval = 0
      if (kbgn .gt. nz) go to 50
      do 40 k = kbgn,nz
         i = ia(k)
         j = ja(k)
         if (iwksp(i) .eq. 1 .or. iwksp(i) .eq. 3 .or.
     a       iwksp(j) .ge. 2) go to 40
         nval = nval + 1
         iwksp(i) = iwksp(i) + 1
         iwksp(j) = iwksp(j) + 2
         if (k .eq. krep) go to 35
         at = a(krep)
         it = ia(krep)
         jt = ja(krep)
         a(krep) = a(k)
         ia(krep) = i
         ja(krep) = j
         a(k) = at
         ia(k) = it 
         ja(k) = jt 
 35      krep = krep + 1
         if (nval .ge. n) go to 45
 40   continue
 45   kbgn = krep
      go to 25
 50   m = mm - 1
      return
      end 
      subroutine prich (nn,r,z)
      implicit double precision (a-h, o-z)
c
c ... prich does the richardson preconditioning.
c
c ... parameters -- 
c
c         n       order of system (= nn)
c         r       residual
c         z       output vector
c
c ... specifications for parameters
c
      dimension r(1), z(1)
c
      n = nn
      do 10 i = 1,n 
 10   z(i) = r(i)
      return
      end 
      subroutine pstops (nn,r,z,u,ubar,ier)
      implicit double precision (a-h, o-z)
c
c ... pstops performs a test to see if the iterative method has
c     converged to a solution inside the error tolerance, zeta.
c     (cg and si routines)
c
c     the stopping tests are --
c
c    (1)  (emax/emin) * sqrt ( (r ,zt)/(rhs,inv(q)*rhs) )
c    (2)  ( 1.0/emin) * sqrt ( (zt,zt)/(u,u) )
c    (3)  (emax/emin) * sqrt ( (zt,zt)/(inv(q)*rhs,inv(q)*rhs) )
c    (4)                sqrt ( (zt,zt)/(inv(q)*rhs,inv(q)*rhs) )
c    (5)                sqrt ( (r ,r )/(rhs,rhs) )
c    (6)                sqrt ( (u-ubar,u-ubar)/(ubar,ubar) )
c    (7)  (emax/emin) * sqrt ( (r,z)/(rhs,inv(ql)*rhs) )
c    (8)  ( 1.0/emin) * sqrt ( (z,z)/(u,u) )
c    (9)  (emax/emin) * sqrt ( (z,z)/(inv(ql)*rhs,inv(ql)*rhs) )
c   (10)                sqrt ( (z,z)/(inv(ql)*rhs,inv(ql)*rhs) )
c
c
c ... parameters -- 
c
c         n       order of system
c         r       residual vector
c         z       pseudo-residual vector
c         u       solution estimate
c         ier     error flag
c                  =  0   no errors detected
c                  = -7   splitting matrix is not positive definite
c
c ... specifications for parameters
c
      dimension r(1), z(1), u(1), ubar(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 / srelpr, keyzer, keygs
      common / itcom9 / rdot, rzdot, rztdot, zdot, zztdot, ztdot,
     a           rhave, zhave, zthave, rcalp, zcalp, ztcalp,
     a           udhav, rdhav, rzhav, rzthav, zdhav, zzthav, ztdhav
      logical rhave, zhave, zthave, rcalp, zcalp, ztcalp,
     a        udhav, rdhav, rzhav, rzthav, zdhav, zzthav, ztdhav
c
c *** end   -- package common 
c
      logical q1
      save    q1
      n = nn
      halt = .false.
      tiny = 500.0d0*srelpr
      nteste = ntest
      if (ntest .gt. 6) nteste = nteste - 6
c
      go to (10,20,30,40,50,60), nteste 
c
c ... test 1
c
 10   if (rzdot .ge. 0.0d0) go to 15
      ier = -7
      call ershow (ier,'pstops')
      return
 15   emaxl = emax
      eminl = emin
      if (eminl .lt. tiny) eminl = tiny 
      tl = emaxl*sqrt (rzdot) 
      tr = eminl*bnorm1
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
c
c ... test 2
c
c ... special procedure for zeroth iteration
c
 20   if (in .ge. 1) go to 25 
      q1 = .false.
      udnm = 1.0d0
      stptst = sqrt (rzdot)
      if (stptst .lt. tiny) halt = .true.
      return
c
c ... in .ge. 1
c
c ... test if udnm needs to be recomputed.
c
 25   if (q1) go to 28
      if ((in .gt. 5)  .and.  (mod(in,5) .ne. 0)) go to 28
         uold = udnm
         udnm = sqrt ( vdot (n,u,u) )
         if (udnm .lt. tiny) udnm = 1.0d0 
         if ((in .gt. 5)  .and.
     a       (abs (udnm-uold) .lt. udnm*zeta)) q1 = .true.
c
c ... compute stopping test.
c
 28   eminl = emin
      if (eminl .lt. tiny) eminl = tiny 
      tl = sqrt ( vdot (n,z,z) )
      tr = udnm*eminl
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
c
c ... test 3.
c
 30   emaxl = emax
      eminl = emin
      if (eminl .lt. tiny) eminl = tiny 
      tl = emaxl*sqrt ( vdot (n,z,z) )
      tr = eminl*bnorm1
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
c
c ... test 4.
c
 40   tl = sqrt ( vdot (n,z,z) )
      tr = bnorm1
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
c
c ... test 5.
c
 50   tl = sqrt ( vdot (n,r,r) )
      tr = bnorm
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
c
c ... test 6.
c
 60   sum = 0.0d0
      do 65 i = 1,n 
 65   sum = sum + (u(i) - ubar(i))**2
      tl = sqrt (sum)
      tr = ubarnm
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
      end 
      subroutine rowise (maxnz,jcoef,irwise)
      implicit double precision (a-h, o-z)
c
c ... rowise determines whether a row-wise or diagonal-wise 
c     algorithm should be used for ic and ssor splittings with
c     diagonal storage.  this routine should be called after
c     final factorization is computed.
c
c ... parameters -- 
c
c          maxnz  number of number of diagonals stored
c          jcoef  vector of diagonal numbers for factorization
c                  array or matrix
c          irwise has a value upon output of
c                  0   if diagonal-wise algorithm should be used
c                  1   if row-wise algorithm should be used 
c
c ... specifications for parameters
c
      integer   jcoef(2)
c
c ... use a rowwise algorithm if  2 .le. /jcoef(j)/ .le. maxd
c     some j.
c
      maxd = 10
c
      irwise = 0
      do 15 j = 1,maxnz
         jcol = iabs(jcoef(j))
         if (jcol .le. 1 .or. jcol .gt. maxd) go to 15
         irwise = 1 
         return
 15   continue
      return
      end 
      subroutine rowsum (lda,n,maxnzz,a,x,isym)
      implicit double precision (a-h, o-z)
c
c ... rowsum computes the row sum of the matrix a.
c
c ... parameters -- 
c
c        lda     leading dimension of array a
c        n       active size of array a 
c        maxnz   number of columns in array a
c        a       array of size n by maxnz
c        x       vector of length n containing the row
c                 sum of a upon output
c        isym    symmetry switch
c                 = 0  matrix is a banded symmetric matrix
c                       with the diagonal in column one
c                 = 1  matrix is nonsymmetric
c
c ... specifications for parameters
c
      dimension a(lda,1), x(1)
c
      maxnz = maxnzz
      do 10 i = 1,n 
 10   x(i) = 0.0d0
      do 20 j = 1,maxnz
         do 15 i = 1,n
 15      x(i) = x(i) + a(i,j) 
 20   continue
      if (isym .eq. 1 .or. maxnz .le. 1) return
      do 30 j = 2,maxnz
         do 25 i = j,n
 25      x(i) = x(i) + a(i-j+1,j)
 30   continue
      return
      end 
      subroutine rsad (nn,nsize,nrr,ndim,maxnew,ndtt,ndbb,jcnew,
     a                 coef,c,b,dfac,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsad computes  c = (dr - t*inv(db)*b)*b
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c        nsize      size of an individual subsystem (if multiple
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        ndt        number of upper diagonals in diagonal block
c        ndb        number of lower diagonals in diagonal block
c        coef       floating point data structure 
c        b          vector of length n containing bb behind br
c        c          vector of length nr containing cr
c        dfac       vector of length (1+nt+nb)*n to contain 
c                    factorization of diagonal block upon output
c        wksp       workspace vector of length nb 
c
c ... specifications for parameters
c
      integer   jcnew(2,1), maxnew(2)
      dimension coef(ndim,2), b(1), c(1), dfac(1), wksp(1)
c
      n = nn
      nr = nrr
      ndt = ndtt
      ndb = ndbb
      nrp1 = nr + 1 
      nb = n - nr
      maxd = 1 + ndt + ndb
      maxz = maxnew(1) - maxd 
      max2 = maxnew(2) - maxd 
c
c ... cr = dr*br.
c
      if (ndt+ndb .gt. 0) go to 15
      do 10 i = 1,nr
 10   c(i) = coef(i,1)*b(i)
      go to 20
 15   call bmuln (ndim,nr,ndt,ndb,coef,coef(1,2),coef(1,ndt+2),b,c)
c
c ... wksp = b*br
c
 20   if (maxz*max2 .eq. 0) return
      do 25 i = 1,nb
 25   wksp(i) = 0.0d0 
      call vaddd (ndim,2,nb,nr,max2,coef(nrp1,maxd+1),
     a            jcnew(2,maxd+1),wksp,b,-nr)
c
c ... wksp = inv(db)*wksp
c
      if (ndt+ndb .gt. 0) go to 35
      do 30 i = 1,nb
 30   wksp(i) = wksp(i)*dfac(i+nr)
      go to 40
 35   call bdsol (n,nb,nsize,ndt,ndb,dfac(nrp1),wksp,wksp,1)
c
c ... cr = cr - t*wksp
c
 40   call vsubd (ndim,2,nr,nb,maxz,coef(1,maxd+1),jcnew(1,maxd+1),
     a            c,wksp,nr)
      return
      end 
      subroutine rsap (ndimm,n,nr,maxnz,jcoef,coef,b,c,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsap computes  c = (dr - t*inv(db)*b)*b
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     purdue format 
c
c ... parameters -- 
c
c        ndim       row dimension of coef,jcoef arrays
c        n          order of total system
c        nr         order of red subsystem
c        maxnz      number of columns in coef,jcoef arrays
c        jcoef      integer array of matrix column numbers
c        coef       floating point array of matrix coefficients
c        b,c        vectors of length nr
c        wksp       workspace array of length n + nb
c
c ... specifications for parameters
c
      integer   jcoef(ndimm,2)
      dimension coef(ndimm,2), b(1), c(1), wksp(1)
c
      ndim = ndimm
      do 10 i = 1,nr
 10   c(i) = coef(i,1)*b(i)
      if (maxnz .le. 1) return
      np1 = n + 1
      nb = n - nr
      nrp1 = nr + 1 
      maxm1 = maxnz - 1
      do 15 i = 1,n 
 15   wksp(i) = 0.0d0 
      call vaddp (ndim,ndim,nb,maxm1,coef(nrp1,2),jcoef(nrp1,2),
     a            wksp(nrp1),b,wksp(np1))
      do 20 i = nrp1,n
 20   wksp(i) = wksp(i)/coef(i,1)
      call vsubp (ndim,ndim,nr,maxm1,coef(1,2),jcoef(1,2),c,wksp,wksp)
      return
      end 
      subroutine rsatd (nn,nsize,nrr,ndim,maxnew,ndtt,ndbb,jcnew,
     a                  coef,c,b,dfac,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsatd computes  c = ((dr**t) - (b**t)*(db**(-t))*(t**t))*b
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c        nsize      size of an individual subsystem (if multiple
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        ndt        number of upper diagonals in diagonal block
c        ndb        number of lower diagonals in diagonal block
c        coef       floating point data structure 
c        b          vector of length n containing bb behind br
c        c          vector of length nr containing cr
c        dfac       vector of length (1+nt+nb)*n to contain 
c                    factorization of diagonal block upon output
c        wksp       workspace vector of length nb 
c
c ... specifications for parameters
c
      integer   jcnew(2,1), maxnew(2)
      dimension coef(ndim,2), b(1), c(1), dfac(1), wksp(1)
c
      n = nn
      nr = nrr
      ndt = ndtt
      ndb = ndbb
      nrp1 = nr + 1 
      nb = n - nr
      maxd = 1 + ndt + ndb
      maxz = maxnew(1) - maxd 
      max2 = maxnew(2) - maxd 
c
c ... cr = (dr**t)*br.
c
      if (ndt+ndb .gt. 0) go to 15
      do 10 i = 1,nr
 10   c(i) = coef(i,1)*b(i)
      go to 20
 15   call bmulnt (ndim,nr,ndt,ndb,coef,coef(1,2),coef(1,ndt+2),b,c)
c
c ... wksp = (t**t)*br
c
 20   if (maxz*max2 .eq. 0) return
      do 25 i = 1,nb
 25   wksp(i) = 0.0d0 
      call vadddt (ndim,2,nr,nb,maxz,coef(1,maxd+1),
     a            jcnew(1,maxd+1),wksp,b,nr)
c
c ... wksp = (db**(-t))*wksp
c
      if (ndt+ndb .gt. 0) go to 35
      do 30 i = 1,nb
 30   wksp(i) = wksp(i)*dfac(i+nr)
      go to 40
 35   call bdsolt (n,nb,nsize,ndt,ndb,dfac(nrp1),wksp,wksp) 
c
c ... cr = cr - (b**t)*wksp
c
 40   call vsubdt (ndim,2,nb,nr,max2,coef(nrp1,maxd+1),
     a            jcnew(2,maxd+1),c,wksp,-nr)
      return
      end 
      subroutine rsatp (ndimm,n,nr,maxnz,jcoef,coef,b,c,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsatp computes  c = (dr - (b**t)*inv(db)*(t**t))*b
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     purdue format 
c
c ... parameters -- 
c
c        ndim       row dimension of coef,jcoef arrays
c        n          order of total system
c        nr         order of red subsystem
c        maxnz      number of columns in coef,jcoef arrays
c        jcoef      integer array of matrix column numbers
c        coef       floating point array of matrix coefficients
c        b,c        vectors of length nr
c        wksp       workspace array of length n + nb
c
c ... specifications for parameters
c
      integer   jcoef(ndimm,2)
      dimension coef(ndimm,2), b(1), c(1), wksp(1)
c
      ndim = ndimm
      do 10 i = 1,nr
 10   c(i) = coef(i,1)*b(i)
      if (maxnz .le. 1) return
      np1 = n + 1
      nb = n - nr
      nrp1 = nr + 1 
      maxm1 = maxnz - 1
      do 15 i = 1,n 
 15   wksp(i) = 0.0d0 
      call vaddpt (ndim,ndim,nr,maxm1,coef(1,2),jcoef(1,2),wksp,b,
     a             wksp)
      do 20 i = nrp1,n
 20   wksp(i) = wksp(i)/coef(i,1)
      call vsubpt (ndim,ndim,nb,maxm1,coef(nrp1,2),jcoef(nrp1,2),c,
     a             wksp(nrp1),wksp(np1))
      return
      end 
      subroutine rsbegd (nn,nsize,nrr,ndim,maxnew,ndtt,ndbb,jcnew,
     a                   coef,c,b,dfac,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsbegd computes  cr = br - t*inv(db)*bb.
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c        nsize      size of an individual subsystem (if multiple
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        ndt        number of upper diagonals in diagonal block
c        ndb        number of lower diagonals in diagonal block
c        coef       floating point data structure 
c        b          vector of length n containing bb behind br
c        c          vector of length nr containing cr
c        dfac       vector of length (1+nt+nb)*n containing 
c                    factorization of diagonal block upon input
c        wksp       workspace vector of length nb 
c
c ... specifications for parameters
c
      integer   jcnew(2,1), maxnew(2)
      dimension coef(ndim,2), b(1), c(1), dfac(1), wksp(1)
c
      n = nn
      nr = nrr
      ndt = ndtt
      ndb = ndbb
      nrp1 = nr + 1 
      nb = n - nr
      maxd = 1 + ndt + ndb
c
c ... compute cr.
c
      do 10 i = 1,nr
 10   c(i) = b(i)
      call bdsol (n,nb,nsize,ndt,ndb,dfac(nrp1),b(nrp1),wksp,1)
      maxm1 = maxnew(1) - maxd
      call vsubd (ndim,2,nr,nb,maxm1,coef(1,maxd+1),jcnew(1,maxd+1),
     a            c,wksp,nr)
      return
      end 
      subroutine rsbegp (n,nr,ndim,maxnz,jcoef,coef,c,b,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsbegp computes  cr = br - t*inv(db)*bb.
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     purdue storage
c
c ... parameters -- 
c
c        n          order of system
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnz      number of columns in coef array
c        jcoef      integer data structure
c        coef       floating point data structure 
c        b          vector of length n containing bb behind br
c        c          vector of length nr containing cr
c        wksp       workspace vector of length n
c
c ... specifications for parameters
c
      integer   jcoef(ndim,2) 
      dimension coef(ndim,2), b(1), c(1), wksp(1) 
c
      nrp1 = nr + 1 
      do 10 i = 1,nr
 10   c(i) = b(i)
      if (maxnz .le. 1) return
      do 15 i = nrp1,n
 15   wksp(i) = b(i)/coef(i,1)
      maxm1 = maxnz - 1
      call vsubp (ndim,ndim,nr,maxm1,coef(1,2),jcoef(1,2),c,
     a            wksp,wksp)
      return
      end 
      subroutine rsendd (nn,nsize,nrr,ndim,maxnew,ndtt,ndbb,jcnew,
     a                   coef,x,b,dfac) 
      implicit double precision (a-h, o-z)
c
c ... rsendd computes  xb = inv(db)*(bb - b*xr)
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c        nsize      size of an individual subsystem (if multiple
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        ndt        number of upper diagonals in diagonal block
c        ndb        number of lower diagonals in diagonal block
c        coef       floating point data structure 
c        x          vector of length n containing  xr, xb
c        b          vector of length n containing bb in the last
c                    nb locations
c        dfac       vector of length (1+nt+nb)*n containing 
c                    factorization of diagonal block upon input
c
c ... specifications for parameters
c
      integer   jcnew(2,1), maxnew(2)
      dimension coef(ndim,2), x(1), b(1), dfac(1) 
c
      n = nn
      nr = nrr
      ndt = ndtt
      ndb = ndbb
      nrp1 = nr + 1 
      nb = n - nr
      maxd = 1 + ndt + ndb
c
c ... compute xb.
c
      do 10 i = nrp1,n
 10   x(i) = b(i)
      max2 = maxnew(2) - maxd 
      call vsubd (ndim,2,nb,nr,max2,coef(nrp1,maxd+1),
     a            jcnew(2,maxd+1),x(nrp1),x,-nr)
      call bdsol (n,nb,nsize,ndt,ndb,dfac(nrp1),x(nrp1),x(nrp1),1)
      return
      end 
      subroutine rsendp (n,nr,ndim,maxnz,jcoef,coef,x,b,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsendp computes  xb = inv(db)*(bb - b*xr)
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     purdue format 
c
c ... parameters -- 
c
c        n          order of matrix
c        nr         order of red subsystem
c        ndim       row dimension of ah and jah arrays
c        maxnz      number of columns in coef and jcoef arrays
c        jcoef      integer array of column numbers
c        coef       floating point array of matrix coefficients
c        x          vector of length n containing  xr, xb
c        b          vector of length n containing bb in the last
c                    nb locations
c        wksp       workspace array of length nb
c
c ... specifications for parameters
c
      integer   jcoef(ndim,2) 
      dimension coef(ndim,2), x(1), b(1), wksp(1) 
c
      nrp1 = nr + 1 
      nb = n - nr
      do 10 i = nrp1,n
 10   x(i) = b(i)
      if (maxnz .le. 1) go to 15
      maxm1 = maxnz - 1
      call vsubp (ndim,ndim,nb,maxm1,coef(nrp1,2),jcoef(nrp1,2),
     a            x(nrp1),x,wksp)
 15   do 20 i = nrp1,n
 20   x(i) = x(i)/coef(i,1)
      return
      end 
      subroutine rsmatd (ndim,nrr,nb,maxnew,jcnew,dr,ah,ak,db,
     a                   maxrss,jcrs,rs,maxlim,isym,ier)
      implicit double precision (a-h, o-z)
c
c ... rsmatd computes  rs = dr - ah*inv(db)*ak  where a has been
c     permuted to red-black form --
c
c                   * dr  ah *
c             a =   *        *
c                   * ak  db *
c
c     (diagonal storage)
c
c      dr is nr x nr        ah is nr x nb
c      ak is nb x nr        db is nb x nb
c
c ... definition of parameters --
c
c         ndim          row dimension of ah and ak arrays
c         nr            number of red points
c         nb            number of black points
c         maxnew        integer vector of length 2 indicating number
c                        of diagonals stored in ah and ak,
c                        respectively.
c         jcnew         integer array of diagonal numbers
c         dr            vector of length nr
c         ah            array of size nr by (maxnew(1)-1)
c         ak            array of size nb by (maxnew(2)-1)
c         db            vector of length nb
c         maxrs         number of columns needed to store reduced
c                        system (output)
c         jcrs          diagonal numbers for rs (output)
c         rs            array to contain reduced system
c         maxlim        maximum column width to be allowed for rs
c         isym          symmetry switch for rs matrix
c                        = 0   store only upper half of rs
c                        = 1   store all of rs
c         ier           error code
c                        =  0     no errors detected
c                        = -2     maxlim .lt. maxrs
c
c ... specifications for parameters
c
      integer maxnew(2), jcnew(2,1), jcrs(1)
      dimension db(1), ak(ndim,1), ah(ndim,1), dr(1), rs(nrr,1)
c
      nr = nrr
      maxrs = 1
      jcrs(1) = 0
      do 5 i = 1,nr 
 5    rs(i,1) = dr(i)
      maxh = maxnew(1) - 1
      maxk = maxnew(2) - 1
      do 35 lh = 1,maxh
         i = jcnew(1,lh+1) - nr
         ia1 = max (1,1-i)
         ib1 = min (nr,nb-i) 
         do 30 lk = 1,maxk
            k = jcnew(2,lk+1) + nr
            l = i + k
            if (l .lt. 0 .and. isym .eq. 0) go to 30
            do 10 ld = 1,maxrs
               if (jcrs(ld) .eq. l) go to 20
 10         continue
            if (maxrs .eq. maxlim) go to 999
            maxrs = maxrs + 1 
            ld = maxrs
            jcrs(maxrs) = l
            do 15 ii = 1,nr
 15         rs(ii,maxrs) = 0.0d0
 20         ist = max (ia1,1-l)
            ied = min (ib1,nr-l)
            do 25 m = ist,ied 
 25         rs(m,ld) = rs(m,ld) - ah(m,lh)*ak(m+i,lk)/db(m+i)
 30      continue
 35   continue
      maxrss = maxrs
      return
c
c ... error exit -- maxlim too small.
c
 999  ier = -2
      return
      end 
      subroutine rsmatp (ndim,nrr,maxnzz,jcoef,coef,maxrss,jcrs,
     a                   rs,maxlim,wksp,iwksp,ier)
      implicit double precision (a-h, o-z)
c
c ... rsmatp computes  rs = dr - ah*inv(db)*ak  where a has been
c     permuted to red-black form --
c
c                   * dr  ah *
c             a =   *        *
c                   * ak  db *
c
c     (purdue storage)
c
c      dr is nr x nr        ah is nr x nb
c      ak is nb x nr        db is nb x nb
c
c ... definition of parameters --
c
c         ndim          row dimension of coef and jcoef arrays
c         nr            number of red points
c         maxnz         number of columns in coef and jcoef 
c         jcoef         array of column indices
c         coef          array of matrix coefficients
c         maxrs         number of columns needed to store reduced
c                        system (output)
c         jcrs          column numbers for rs (output)
c         rs            array to contain reduced system
c         maxlim        maximum column width to be allowed for rs
c         wksp          workspace of length 2*nr
c         iwksp         integer workspace of length nr
c         ier           error code
c                        =  0     no errors detected
c                        = -2     maxlim .lt. maxrs
c
c ... specifications for parameters
c
      integer jcoef(ndim,1), jcrs(nrr,1), iwksp(1)
      dimension coef(ndim,1), rs(nrr,1), wksp(1)
c
      nr = nrr
      maxnz = maxnzz
      maxrs = 1
      do 5 i = 1,nr 
         rs(i,1) = coef(i,1)
         jcrs(i,1) = i
 5    continue
      do 50 j = 2,maxnz
         call vgathr (nr,coef,jcoef(1,j),wksp)
         do 10 i = 1,nr
 10      wksp(i) = coef(i,j)/wksp(i)
         do 45 jj = 2,maxnz
            call vgathr (nr,coef(1,jj),jcoef(1,j),wksp(nr+1))
            call vgathi (nr,jcoef(1,jj),jcoef(1,j),iwksp)
            do 15 i = 1,nr
 15         wksp(nr+i) = wksp(i)*wksp(nr+i)
            do 40 i = 1,nr
               jcol = iwksp(i)
               term = wksp(nr+i)
               if (jcol .gt. nr) go to 40
               do 20 jjj = 1,maxrs
                  if (jcrs(i,jjj) .ne. jcol) go to 20
                  rs(i,jjj) = rs(i,jjj) - term
                  go to 40
 20            continue
               if (maxrs .eq. 1) go to 30
               do 25 jjj = 2,maxrs
                  if (jcrs(i,jjj) .ne. i) go to 25
                  rs(i,jjj) = rs(i,jjj) - term
                  jcrs(i,jjj) = jcol
                  go to 40
 25            continue
 30            if (maxrs .eq. maxlim) go to 999
               maxrs = maxrs + 1
               do 35 ii = 1,nr
                  jcrs(ii,maxrs) = ii
                  rs(ii,maxrs) = 0.0d0
 35            continue
               rs(i,maxrs) = -term
               jcrs(i,maxrs) = jcol
 40         continue
 45      continue
 50   continue
      maxrss = maxrs
      return
c
c ... error exit -- maxlim too small.
c
 999  ier = -2 
      return
      end 
      subroutine rsrhsd (nn,nrr,ndim,maxnew,jcnew,coef,c,b,wksp)
      implicit double precision (a-h, o-z)
c
c ... rsrhsd computes  cr = br - t*inv(db)*bb.
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        coef       floating point data structure 
c        b          vector of length n containing bb behind br
c        c          vector of length nr containing cr
c        wksp       workspace vector of length nb 
c
c ... specifications for parameters
c
      integer   jcnew(2,2), maxnew(2)
      dimension coef(ndim,2), b(1), c(1), wksp(1) 
c
      n = nn
      nr = nrr
      nb = n - nr
c
c ... compute cr.
c
      do 10 i = 1,nr
 10   c(i) = b(i)
      do 15 i = 1,nb
 15   wksp(i) = b(nr+i)/coef(nr+i,1)
      maxm1 = maxnew(1) - 1
      call vsubd (ndim,2,nr,nb,maxm1,coef(1,2),jcnew(1,2),
     a            c,wksp,nr)
      return
      end 
      subroutine rsxbd (nn,nrr,ndim,maxnew,jcnew,coef,x,b)
      implicit double precision (a-h, o-z)
c
c ... rsxbd computes  xb = inv(db)*(bb - b*xr)
c
c       where      a  = ( dr   t )
c                       ( b   db )
c
c     diagonal storage
c
c ... parameters -- 
c
c        n          order of system
c                    systems) 
c        nr         order of the red subsystem
c        ndim       row dimension of coef array
c        maxnew     number of columns in coef array
c        coef       floating point data structure 
c        x          vector of length n containing  xr, xb
c        b          vector of length n containing bb in the last
c                    nb locations
c
c ... specifications for parameters
c
      integer   jcnew(2,2), maxnew(2)
      dimension coef(ndim,2), x(1), b(1)
c
      n = nn
      nr = nrr
      nrp1 = nr + 1 
      nb = n - nr
c
c ... compute xb.
c
      do 10 i = nrp1,n
 10   x(i) = b(i)
      max2 = maxnew(2) - 1
      call vsubd (ndim,2,nb,nr,max2,coef(nrp1,2), 
     a            jcnew(2,2),x(nrp1),x,-nr)
      do 15 i = nrp1,n
 15   x(i) = x(i)/coef(i,1)
      return
      end 
      subroutine sbbs (ldd,ldt,n,kblszz,nsize,lbhb,iblock,d,t,
     a                 jt,x,omega)
      implicit double precision (a-h, o-z)
c
c ... sbbs does an block ssor backward pass.
c     symmetric diagonal data structure, natural ordering.
c     block ssor 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         kblsz    block size 
c         nsize    size of an individual subsystem within a 
c                   diagonal block
c         lbhb     number of blocks per block row 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         x        input/output vector of length n
c         omega    over-relaxation factor
c
c ... specifications for parameters
c
      integer   jt(1), iblock(3,1)
      dimension d(ldd,2), t(ldt,1), x(1)
c
      kblsz = kblszz
      l = n/kblsz
      nt = iblock(3,1) - 1
      do 35 k = l,1,-1
         ist = (k - 1)*kblsz + 1
         ied = k*kblsz
         if (k .eq. l) go to 15
         jjlim = min (lbhb,l-k+2)
         do 10 jj = 3,jjlim
            jblk = iblock(1,jj)
            jst = iblock(2,jj)
            mjj = iblock(3,jj)
            inc = jblk*kblsz
            istf = ist + inc
            if (istf .gt. n) go to 10
            call vsubd (ldt,1,kblsz,kblsz,mjj,t(ist,jst),jt(jst),
     a                   x(ist),x(istf),inc)
 10      continue
 15      if (nt .ge. 1) go to 25
         do 20 i = ist,ied
 20      x(i) = omega*d(i,1)*x(i)
         go to 35
 25      call bdsol (ldd,kblsz,nsize,nt,0,d(ist,1),x(ist),x(ist),
     a               0)
         do 30 i = ist,ied
 30      x(i) = omega*x(i)
 35   continue
      return
      end 
      subroutine sbbsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif,wksp)
      implicit double precision (a-h, o-z)
c
c ... sbbsn does an block ssor backward solve.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         x        input/output vector of length n
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,2), t(ldt,1), wksp(1), x(1) 
      logical   unif
c
      unif = iunif .eq. 1
c
      l = ncolor
      if (.not. unif) go to 10
      na = nci(1)
      nb = na
      jlim = lbhb(1)
      l = n/na
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      kk = 1
c
c ... do backward solution.
c
 10   lm1 = l - 1
      do 50 k = lm1,1,-1
         if (unif) go to 15
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         go to 20
 15      ist = (k - 1)*na + 1 
 20      ied = ist + na - 1
         do 25 i = 1,na
 25      wksp(i) = 0.0d0
         do 30 j = 3,jlim
            jcol = k + iblock(1,kk,j)
            if (jcol .le. k) go to 30
            jstb = iblock(2,kk,j)
            mb = iblock(3,kk,j)
            if (unif) inc = (jcol - k)*na
            if (.not. unif) inc = ipt(jcol) - ipt(k)
            if (.not. unif) nb = nci(jcol)
            istb = ist + inc
            call vaddd (ldt,ncolor,na,nb,mb,t(ist,jstb),jt(kk,jstb),
     a                  wksp,x(istb),inc)
 30      continue
         if (ndt + ndb .ge. 1) go to 40 
         do 35 i = ist,ied
 35      x(i) = x(i) - omega*d(i,1)*wksp(i-ist+1) 
         go to 50
 40      call bdsol (ldd,na,nsize,ndt,ndb,d(ist,1),wksp,wksp,1)
         do 45 i = ist,ied
 45      x(i) = x(i) - omega*wksp(i-ist+1)
 50   continue
      return
      end 
      subroutine sbbsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                   iblock,d,t,jt,x,omega,iunif)
      implicit double precision (a-h, o-z)
c
c ... sbbsnt does an block ssor transpose backward solve.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         x        input/output vector of length n
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,2), t(ldt,1), x(1)
      logical   unif
c
      unif = iunif .eq. 1
c
      l = ncolor
      if (.not. unif) go to 10
      na = nci(1)
      nb = na
      jlim = lbhb(1)
      l = n/na
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      kk = 1
c
c ... do backward solution.
c
 10   do 50 k = l,1,-1
         if (unif) go to 15
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         go to 20
 15      ist = (k - 1)*na + 1 
 20      ied = ist + na - 1
         if (ndt + ndb .ge. 1) go to 30 
         do 25 i = ist,ied
 25      x(i) = omega*d(i,1)*x(i)
         go to 40
 30      call bdsolt (ldd,na,nsize,ndt,ndb,d(ist,1),x(ist),x(ist))
         do 35 i = ist,ied
 35      x(i) = omega*x(i)
 40      do 45 j = 3,jlim
            jcol = k + iblock(1,kk,j)
            if (jcol .ge. k) go to 45
            jstb = iblock(2,kk,j)
            mb = iblock(3,kk,j)
            if (unif) inc = (jcol - k)*na
            if (.not. unif) inc = ipt(jcol) - ipt(k)
            if (.not. unif) nb = nci(jcol)
            istb = ist + inc
            call vsubdt (ldt,ncolor,na,nb,mb,t(ist,jstb),jt(kk,jstb), 
     a                  x(istb),x(ist),inc)
 45      continue
 50   continue
      return
      end 
      subroutine sbfs (ldd,ldt,n,kblszz,nsize,lbhb,iblock,d,t,
     a                 jt,x,omega,wksp) 
      implicit double precision (a-h, o-z)
c
c ... sbfs does an block ssor forward pass.
c     symmetric diagonal data structure, natural ordering.
c     block ssor 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         kblsz    block size 
c         nsize    size of an individual subsystem within a 
c                   diagonal block
c         lbhb     number of blocks per block row 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         x        input/output vector of length n
c         omega    over-relaxation factor
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   jt(1), iblock(3,1)
      dimension d(ldd,2), t(ldt,1), wksp(1), x(1) 
c
      kblsz = kblszz
      l = n/kblsz
      lm1 = l - 1
      nt = iblock(3,1) - 1
      do 35 k = 1,lm1
         ist = (k - 1)*kblsz + 1
         ied = k*kblsz
         if (nt .ge. 1) go to 15
         do 10 i = ist,ied
 10      wksp(i-ist+1) = omega*d(i,1)*x(i)
         go to 25
 15      call bdsol (ldd,kblsz,nsize,nt,0,d(ist,1),
     a                          x(ist),wksp,0)
         do 20 i = 1,kblsz
 20      wksp(i) = omega*wksp(i)
 25      jjlim = min (lbhb,l-k+2)
         do 30 jj = 3,jjlim
            jblk = iblock(1,jj)
            jst = iblock(2,jj)
            mjj = iblock(3,jj)
            inc = jblk*kblsz
            istf = ist + inc
            if (istf .gt. n) go to 30
            call vsubdt (ldt,1,kblsz,kblsz,mjj,t(ist,jst),jt(jst),
     a                   x(istf),wksp,inc)
 30      continue
 35   continue
      return
      end 
      subroutine sbfsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif)
      implicit double precision (a-h, o-z)
c
c ... sbfsn does an block ssor forward solve.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         x        input/output vector of length n
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,2), t(ldt,1), x(1)
      logical   unif
c
      unif = iunif .eq. 1
c
      l = ncolor
      if (.not. unif) go to 10
      na = nci(1)
      nb = na
      jlim = lbhb(1)
      l = n/na
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      kk = 1
c
c ... do forward solution.
c
 10   do 45 k = 1,l 
         if (unif) go to 15
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         go to 20
 15      ist = (k - 1)*na + 1 
 20      ied = ist + na - 1
         do 25 j = 3,jlim
            jcol = k + iblock(1,kk,j)
            if (jcol .ge. k) go to 25
            jstb = iblock(2,kk,j)
            mb = iblock(3,kk,j)
            if (unif) inc = (jcol - k)*na
            if (.not. unif) inc = ipt(jcol) - ipt(k)
            if (.not. unif) nb = nci(jcol)
            istb = ist + inc
            call vsubd (ldt,ncolor,na,nb,mb,t(ist,jstb),jt(kk,jstb),
     a                  x(ist),x(istb),inc)
 25      continue
         if (ndt + ndb .ge. 1) go to 35 
         do 30 i = ist,ied
 30      x(i) = omega*d(i,1)*x(i)
         go to 45
 35      call bdsol (ldd,na,nsize,ndt,ndb,d(ist,1),x(ist),x(ist),1)
         do 40 i = ist,ied
 40      x(i) = omega*x(i)
 45   continue
      return
      end 
      subroutine sbfsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                   iblock,d,t,jt,x,omega,iunif,wksp)
      implicit double precision (a-h, o-z)
c
c ... sbfsnt does an block ssor transpose forward solve.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         x        input/output vector of length n
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,2), t(ldt,1), wksp(1), x(1) 
      logical   unif
c
      unif = iunif .eq. 1
c
      l = ncolor
      if (.not. unif) go to 10
      na = nci(1)
      nb = na
      jlim = lbhb(1)
      l = n/na
      ndt = iblock(3,1,1) - 1 
      ndb = iblock(3,1,2)
      kk = 1
c
c ... do forward solution.
c
 10   lm1 = l - 1
      do 50 k = 1,lm1
         if (unif) go to 15
         kk = k
         ist = ipt(k) + 1
         jlim = lbhb(k)
         na = nci(k)
         ndt = iblock(3,k,1) - 1
         ndb = iblock(3,k,2)
         go to 20
 15      ist = (k - 1)*na + 1 
 20      ied = ist + na - 1
         if (ndt + ndb .ge. 1) go to 30 
         do 25 i = ist,ied
 25      wksp(i-ist+1) = omega*d(i,1)*x(i)
         go to 40
 30      call bdsolt (ldd,na,nsize,ndt,ndb,d(ist,1),x(ist),wksp)
         do 35 i = 1,na
 35      wksp(i) = omega*wksp(i)
 40      do 45 j = 3,jlim
            jcol = k + iblock(1,kk,j)
            if (jcol .le. k) go to 45
            jstb = iblock(2,kk,j)
            mb = iblock(3,kk,j)
            if (unif) inc = (jcol - k)*na
            if (.not. unif) inc = ipt(jcol) - ipt(k)
            if (.not. unif) nb = nci(jcol)
            istb = ist + inc
            call vsubdt (ldt,ncolor,na,nb,mb,t(ist,jstb),jt(kk,jstb), 
     a                  x(istb),wksp,inc)
 45      continue
 50   continue
      return
      end 
      subroutine sbsl (ldd,ldt,n,kblsz,nsize,lbhb,iblock,d,t,
     a                 jt,y,x,omega,wksp)
      implicit double precision (a-h, o-z)
c
c ... sbsl does an block ssor solution. 
c     symmetric diagonal data structure, natural ordering.
c     block ssor 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         kblsz    block size 
c         nsize    size of an individual subsystem within a 
c                   diagonal block
c         lbhb     number of blocks per block row 
c         iblock   integer array of size 3 by lbhb
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
c         jt       integer vector giving the diagonal numbers
c                   for the off-diagonal blocks
c         y        input vector for the right-hand-side
c         x        output vector for the solution to q*x = y
c         omega    over-relaxation factor
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   jt(1), iblock(3,1)
      dimension d(ldd,1), t(ldt,1), wksp(1), x(1), y(1)
c
      const = 2.0d0 - omega
      do 10 i = 1,n 
 10   x(i) = const*y(i)
      call sbfs (ldd,ldt,n,kblsz,nsize,lbhb,iblock,d,t,
     a           jt,x,omega,wksp)
      call sbbs (ldd,ldt,n,kblsz,nsize,lbhb,iblock,d,t,
     a           jt,x,omega)
      return
      end 
      subroutine sbsln (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,y,x,omega,iunif,wksp) 
      implicit double precision (a-h, o-z)
c
c ... sbsln does an block ssor solution.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), wksp(1), x(1), y(1)
c
      const = 2.0d0 - omega
      do 10 i = 1,n 
 10   x(i) = const*y(i)
      call sbfsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif)
      call sbbsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif,wksp)
      return
      end 
      subroutine sbslnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                   iblock,d,t,jt,y,x,omega,iunif,wksp)
      implicit double precision (a-h, o-z)
c
c ... sbslnt does an block ssor transpose solution.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), wksp(1), x(1), y(1)
c
      const = 2.0d0 - omega
      do 10 i = 1,n 
 10   x(i) = const*y(i)
      call sbfsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif,wksp)
      call sbbsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif)
      return
      end 
      subroutine sbsln1 (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,y,x,omega,iunif)
      implicit double precision (a-h, o-z)
c
c ... sbsln1 does an block ssor forward solution. 
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), x(1), y(1)
c
      const = 2.0d0 - omega
      do 10 i = 1,n 
 10   x(i) = const*y(i)
      call sbfsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif)
      return
      end 
      subroutine sbsln2 (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,y,x,omega,iunif,wksp) 
      implicit double precision (a-h, o-z)
c
c ... sbsln2 does an block ssor back solution.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), wksp(1), x(1), y(1)
c
      do 10 i = 1,n 
 10   x(i) = y(i)
      call sbbsn (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif,wksp)
      return
      end 
      subroutine sbsln3 (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                   iblock,d,t,jt,y,x,omega,iunif)
      implicit double precision (a-h, o-z)
c
c ... sbsln3 does an block ssor transpose forward solution. 
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), x(1), y(1)
c
      const = 2.0d0 - omega
      do 10 i = 1,n 
 10   x(i) = const*y(i)
      call sbbsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif)
      return
      end 
      subroutine sbsln4 (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                   iblock,d,t,jt,y,x,omega,iunif,wksp)
      implicit double precision (a-h, o-z)
c
c ... sbsln4 does an block ssor transpose back solution.
c     nonsymmetric diagonal data structure, natural or multi-color
c     orderings, block ssor 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         nsize    size of an individual subsystem within a 
c                   diagonal block
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         ipt      integer pointer vector of length ncolor+1 if
c                   iunif = 0.  formed in the factorization routine.
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         iblock   integer array of size 3 by ncolor by max(lbhb(i))
c                   giving block constants
c         d        array for diagonal block
c         t        array for off-diagonal blocks
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         y        input vector of length n containing right-hand-side
c         x        output vector containing the solution to q*x = y
c         omega    over-relaxation factor
c         iunif    uniform block size switch
c                   = 0   diagonal blocks are not of uniform size
c                   = 1   diagonal blocks are of uniform size
c         wksp     floating point workspace vector
c
c ... specifications for parameters
c
      integer   ipt(1), jt(ncolor,1), nci(1), lbhb(1),
     a          iblock(3,ncolor,2)
      dimension d(ldd,1), t(ldt,1), wksp(1), x(1), y(1)
c
      do 10 i = 1,n 
 10   x(i) = y(i)
      call sbfsnt (ldd,ldt,n,nsize,ncolor,nci,ipt,lbhb,
     a                  iblock,d,t,jt,x,omega,iunif,wksp)
      return
      end 
      subroutine scal1 (nn,ndim,maxnzz,jcoef,coef,rhs,u,ubar,
     a                  diag,work,iflag,ier)
      implicit double precision (a-h, o-z)
c
c ... scal1 scales the original matrix to a unit diagonal matrix.
c     (purdue data structure) 
c     rhs and u vectors are scaled accordingly.  upon output, diag
c     contains the reciprocal square roots of the diagonal elements.
c     it is assumed that the diagonal of the matrix is in column one
c     of coef.
c
c ... parameters -- 
c
c         n       dimension of matrix
c         ndim    row dimension of coef array in defining routine
c         maxnz   number of columns in coef array 
c         jcoef   integer matrix representation array
c         coef    matrix representation array
c         rhs     right hand side of matrix problem
c         u       latest estimate of solution
c         ubar    exact solution (optional)
c         diag    work array of length n (nonvolatile)
c         work    work array of length n (volatile)
c         iflag   flag for ubar
c                  = 0  do not scale ubar
c                  = 1  scale ubar
c         ier     error flag -- on return, values mean
c                      0 -- no errors detected
c                     -4 -- nonpositive diagonal element
c
c ... specifications for parameters
c
      integer   jcoef(ndim,1) 
      dimension coef(ndim,1), rhs(1), u(1), diag(1), work(1),
     a          ubar(1)
c
c *** begin -- package common 
c
      common / itcom4 / srelpr, keyzer, keygs
c
c *** end   -- package common 
c
      n = nn
      maxnz = maxnzz
c
c ... check for positive diagonal entries for each row.
c
      cmin = vmin (n,coef)
      if (cmin .gt. 0.0d0) go to 10
c
c ... fatal error -- nonpositive diagonal element.
c
      ier = -4
      return
c
c ... scale matrix.  store reciprocal square roots
c ... of diagonal entries in diag.
c
 10   do 15 i = 1,n 
 15   diag(i) = sqrt (coef(i,1))
c
c ... scale rhs, u, and ubar. 
c
      do 20 i = 1,n 
 20   u(i) = diag(i)*u(i)
      if (iflag .eq. 0) go to 30
      do 25 i = 1,n 
 25   ubar(i) = diag(i)*ubar(i)
 30   do 35 i = 1,n 
 35   diag(i) = 1.0d0/diag(i)
      do 40 i = 1,n 
 40   rhs(i) = diag(i)*rhs(i) 
      if (keygs .eq. 2) go to 55
c
c ... using gathers.
c
      do 50 j = 1,maxnz
         call vgathr (n,diag,jcoef(1,j),work)
         do 45 i = 1,n
 45      coef(i,j) = diag(i)*coef(i,j)*work(i)
 50   continue
      return
c
c ... not using gathers.
c
 55   do 65 j = 1,maxnz
         do 60 i = 1,n
 60      coef(i,j) = diag(i)*coef(i,j)*diag(jcoef(i,j))
 65   continue
      return
      end 
      subroutine scal2 (nn,ndim,maxnz,jcoef,coef,rhs,u,ubar,
     a                  diag,iflag,ier) 
      implicit double precision (a-h, o-z)
c
c ... scal2 scales the original matrix to a unit diagonal matrix.
c     (diagonal data structure)
c     rhs and u vectors are scaled accordingly.  upon output, diag
c     contains the reciprocal square roots of the diagonal elements.
c     it is assumed that the diagonal of the matrix is in column one
c     of coef.
c
c ... parameters -- 
c
c         n       dimension of matrix
c         ndim    row dimension of coef array in defining routine
c         maxnz   number of columns in coef array 
c         jcoef   integer matrix representation array
c         coef    matrix representation array
c         rhs     right hand side of matrix problem
c         u       latest estimate of solution
c         ubar    exact solution (optional)
c         diag    work array of length n (nonvolatile)
c         iflag   flag for ubar
c                  = 0  do not scale ubar
c                  = 1  scale ubar
c         ier     error flag -- on return, values mean
c                      0 -- no errors detected
c                     -4 -- nonpositive diagonal element
c
c ... specifications for parameters
c
      integer   jcoef(2)
      dimension coef(ndim,1), rhs(1), u(1), diag(1), ubar(1)
c
c
      n = nn
c
c ... check for positive diagonal entries for each row.
c
      cmin = vmin (n,coef)
      if (cmin .gt. 0.0d0) go to 10
c
c ... fatal error -- nonpositive diagonal element.
c
      ier = -4
      return
c
c ... scale matrix.  store reciprocal square roots
c ... of diagonal entries in diag.
c
 10   do 15 i = 1,n 
 15   diag(i) = sqrt (coef(i,1))
c
c ... scale rhs, u, and ubar. 
c
      do 20 i = 1,n 
 20   u(i) = diag(i)*u(i)
      if (iflag .eq. 0) go to 30
      do 25 i = 1,n 
 25   ubar(i) = diag(i)*ubar(i)
 30   do 35 i = 1,n 
 35   diag(i) = 1.0d0/diag(i)
      do 40 i = 1,n 
 40   rhs(i) = diag(i)*rhs(i) 
c
c ... scale matrix. 
c
      do 60 j = 1,maxnz
         ind = jcoef(j)
         len = n - iabs(ind)
         if (ind .lt. 0) go to 50
         do 45 i = 1,len
 45      coef(i,j) = diag(i)*coef(i,j)*diag(i+ind)
         go to 60
 50      do 55 i = 1,len
 55      coef(i-ind,j) = diag(i-ind)*coef(i-ind,j)*diag(i)
 60   continue
      return
      end 
      subroutine scal3 (nn,nz,ia,ja,a,rhs,u,ubar,diag,
     a                  work,iflag,ier) 
      implicit double precision (a-h, o-z)
c
c ... scal3 scales the original matrix to a unit diagonal matrix.
c     (sparse data structure) 
c     rhs and u vectors are scaled accordingly.  upon output, diag
c     contains the reciprocal square roots of the diagonal elements.
c     it is assumed that the diagonal of the matrix is in the
c     n first locations of a. 
c
c ... parameters -- 
c
c         n       dimension of matrix
c         nz      length of ia, ja, and a vectors 
c         a       vector containing matrix coefficients
c         ia      vector of i values
c         ja      vector of j values
c         rhs     right hand side of matrix problem
c         u       latest estimate of solution
c         ubar    exact solution (optional)
c         diag    vector of length n containing the reciprocal
c                  square roots of the diagonal elements upon
c                  output
c         work    workspace vector of length n
c         iflag   flag for ubar
c                  = 0  do not scale ubar
c                  = 1  scale ubar
c         ier     error flag -- on return, values mean
c                      0 -- no errors detected
c                     -4 -- nonpositive diagonal element
c
c ... specifications for parameters
c
      integer   ia(1), ja(1)
      dimension a(1), rhs(1), u(1), diag(1), work(1),
     a          ubar(1)
c
c *** begin -- package common 
c
      common / itcom4 / srelpr, keyzer, keygs
c
c *** end   -- package common 
c
      n = nn
c
c ... check for positive diagonal entries for each row.
c
      cmin = vmin (n,a)
      if (cmin .gt. 0.0d0) go to 10
c
c ... fatal error -- nonpositive diagonal element.
c
      ier = -4
      return
c
c ... scale matrix.  store reciprocal square roots
c ... of diagonal entries in diag.
c
 10   do 15 i = 1,n 
 15   diag(i) = sqrt (a(i))
c
c ... scale rhs, u, and ubar. 
c
      do 20 i = 1,n 
 20   u(i) = diag(i)*u(i)
      if (iflag .eq. 0) go to 30
      do 25 i = 1,n 
 25   ubar(i) = diag(i)*ubar(i)
 30   do 35 i = 1,n 
 35   diag(i) = 1.0d0/diag(i)
      do 40 i = 1,n 
 40   rhs(i) = diag(i)*rhs(i) 
      if (keygs .eq. 2) go to 60
c
c ... using gathers.
c
      ist = 1
 45   ied = min (ist-1+n,nz) 
      if (ied .lt. ist) return
         len = ied - ist + 1
         call vgathr (len,diag,ia(ist),work)
         do 50 i = ist,ied
 50      a(i) = a(i)*work(i-ist+1)
         call vgathr (len,diag,ja(ist),work)
         do 55 i = ist,ied
 55      a(i) = a(i)*work(i-ist+1)
      ist = ied + 1 
      go to 45
c
c ... not using gathers.
c
 60   do 65 i = 1,nz
 65   a(i) = a(i)*diag(ia(i))*diag(ja(i))
      return
      end 
      subroutine sorstp (n,u,ubar,dnrm,ccon)
      implicit double precision (a-h, o-z)
c
c ... sorstp performs a test to see if the sor
c     method has converged to a solution inside the error
c     tolerance, zeta.
c
c ... parameters -- 
c
c          n      order of system
c          u      present solution estimate
c          ubar   exact solution
c          dnrm   inner product of pseudo-residuals at preceding
c                    iteration
c          con    stopping test parameter (= ccon)
c
c ... specifications for parameters
c
      dimension u(1), ubar(1) 
      logical q1
      save    q1
c
c *** begin -- itpack 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
c
c *** end   -- itpack common
c
      con = ccon
      halt = .false.
      if (ntest .eq. 6) go to 25
c
c ... special procedure for zeroth iteration.
c
      if (in .ge. 1) go to 5
      q1 = .false.
      udnm = 1.0d0
      stptst = 1000.0d0
      return
c
c ... test if udnm needs to be recomputed
c
 5    if (q1) go to 15
      if ((in .gt. 5)  .and.  (mod(in,5) .ne. 0)) go to 15
      uold = udnm
      udnm = 0.0d0
      do 10 i = 1,n 
 10   udnm = udnm + u(i)*u(i) 
      if (udnm .eq. 0.0d0) udnm = 1.0d0
      if ((in .gt. 5) .and.
     a       (abs (udnm-uold) .le. udnm*zeta)) q1 = .true.
c
c ... compute stopping test
c
 15   tr = sqrt (udnm)
      tl = 1.0d0
      if (con .eq. 1.0d0)  go to 20
      tl = sqrt (dnrm)
      tr = tr*(1.0d0 - con)
 20   stptst = tl/tr
      if (tl .ge. tr*zeta) return
      halt = .true. 
      return
c
c ... second test.
c
 25   if (in .eq. 0) ubarnm = sqrt (vdot(n,ubar,ubar))
      sum = 0.0d0
      do 30 i = 1,n 
 30   sum = sum + (u(i) - ubar(i))**2
      tl = sqrt (sum)
      tr = ubarnm
      stptst = tl/tr
      if (tl .lt. tr*zeta) halt = .true.
      return
      end 
      subroutine sords (ndim,nn,maxtt,jt,d,t,omegaa,irwise, 
     a                  u,rhs,unew,iwksp)
      implicit double precision (a-h, o-z)
c
c ... sords does an sor solve (natural ordering,
c     symmetric diagonal storage).
c
c        unew = inv(d + w*l)*((1-w)*d*un + w*(rhs - u*un))
c
c ... parameters -- 
c
c        ndim   row dimension of t array
c        n      order of system (= nn)
c        maxt   number of columns in t array
c        jt     integer vector of length maxt giving the diagonal
c                indices of the corresponding columns in t
c        d      vector of length n giving the diagonal elements
c                of the matrix
c        t      array of active size n by maxt giving the super-
c                diagonals of the matrix
c        omega  over-relaxation factor
c        irwise rowwise algorithm switch
c                = 0  use diagonal algorithm
c                = 1  use row-wise algorithm
c        u      current solution vector 
c        rhs    right hand side
c        unew   updated solution vector 
c        iwksp  integer workspace of length maxt
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), u(1), unew(1), rhs(1)
      integer   jt(1), iwksp(1)
c
c
      n = nn
      maxt = maxtt
      omega = omegaa
c
c ... rhs = (1-w)*d*un + w*(rhs - u*un) 
c
      call vsubd (ndim,1,n,n,maxt,t,jt,rhs,u,0)
      con = 1.0d0 - omega
      do 10 i = 1,n 
 10   rhs(i) = con*d(i)*u(i) + omega*rhs(i)
c
c ... rhs = inv(i+w*l*inv(d))*rhs
c
c ... select rowwise or diagonal-wise algorithm.
c
      if (irwise .eq. 1) go to 50
c
c ... diagonal-wise algorithm.
c
      do 15 i = 1,maxt
 15   iwksp(i) = jt(i) + 1
c
c ... determine nc, imin.
c
 20   nc = n
      do 25 i = 1,maxt
         nterm = iwksp(i) - 1 
         if (nterm .ge. nc) go to 25
         nc = nterm 
         imin = i
 25   continue
      if (nc .ge. n) go to 70 
      ndel = jt(imin)
      ibeg = nc + 1 
      if (ndel .gt. 1) go to 40
c
c ... special case for first minor subdiagonal.
c
      nc1 = n
      do 30 i = 1,maxt
         if (i .eq. imin) go to 30
         if (iwksp(i) .lt. nc1) nc1 = iwksp(i)
 30   continue
      iwksp(imin) = nc1 + 1
      do 35 j = ibeg,nc1
 35   rhs(j) = rhs(j) - omega*t(j-1,imin)*rhs(j-1)/d(j-1)
      go to 20
c
c ... far diagonals  (do vector computations).
c
 40   iwksp(imin) = iwksp(imin) + ndel
      iend = min (ibeg+ndel-1,n)
      do 45 i = ibeg,iend
 45   rhs(i) = rhs(i) - omega*t(i-ndel,imin)*rhs(i-ndel)/d(i-ndel)
      go to 20
c
c ... rowwise algorithm.
c
 50   do 65 i = 1,n 
         do 55 j = 1,maxt
 55      iwksp(j) = min (n,i+jt(j))
         term = omega*rhs(i)/d(i)
         do 60 j = 1,maxt
 60      rhs(iwksp(j)) = rhs(iwksp(j)) - t(i,j)*term
 65   continue
c
c ... unew = inv(d)*rhs
c
 70   do 75 i = 1,n 
 75   unew(i) = rhs(i)/d(i)
      return
      end 
      subroutine sordn (ndim,nn,maxtt,maxbb,jt,jb,d,t,b,omegaa,
     a                  irwise,u,rhs,unew,iwksp)
      implicit double precision (a-h, o-z)
c
c ... sordn does an sor solve (natural ordering,
c     nonsymmetric diagonal storage).
c
c        unew = inv(d + w*l)*((1-w)*d*un + w*(rhs - u*un))
c
c ... parameters -- 
c
c        ndim   row dimension of t array
c        n      order of system (= nn)
c        maxt   number of columns in t array
c        maxb   number of columns in b array
c        jt     integer vector of length maxt giving the diagonal
c                indices of the corresponding columns in t
c        jb     integer vector of length maxb giving the diagonal
c                indices of the corresponding columns in b
c        d      vector of length n giving the diagonal elements
c                of the matrix
c        t      array of active size n by maxt giving the super-
c                diagonals of the matrix
c        b      array of active size n by maxb giving the sub-
c                diagonals of the matrix
c        omega  over-relaxation factor
c        irwise rowwise algorithm switch
c                = 0  use diagonal algorithm
c                = 1  use row-wise algorithm
c        u      current solution vector 
c        rhs    right hand side
c        unew   updated solution vector 
c        iwksp  integer workspace of length maxt
c
c ... specifications for parameters
c
      dimension d(1), t(ndim,1), b(ndim,1), u(1), unew(1), rhs(1)
      integer   jt(1), jb(1), iwksp(1)
c
c
      n = nn
      maxt = maxtt
      maxb = maxbb
      omega = omegaa
c
c ... rhs = (1-w)*d*un + w*(rhs - u*un) 
c
      call vsubd (ndim,1,n,n,maxt,t,jt,rhs,u,0)
      con = 1.0d0 - omega
      do 10 i = 1,n 
 10   rhs(i) = con*d(i)*u(i) + omega*rhs(i)
c
c