subroutine nspcg (precon,accel,ndimm,mdimm,nn,maxnzz,coef,
     a                  jcoef,p,ip,u,ubar,rhs,wksp,iwksp,nw,inw,
     a                  iparm,rparm,ier)
c
c ... nspcg is the driver for the nspcg package.
c
c ... parameters -- 
c
c       precon    preconditioning module
c       accel     acceleration module
c       coef      real matrix data array
c       jcoef     integer matrix data array
c       n         input integer.  order of the system (= nn)
c       u         input/output vector.  on input, u contains the
c                 initial guess to the solution.  on output, it
c                 contains the latest estimate to the solution.
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       wksp      real workspace vector of length nw
c       iwksp     integer workspace vector of length inw
c       nw        length of wksp upon input, amount used upon
c                  output
c       inw       length of iwksp upon input, amount used upon
c                  output
c       iparm     integer vector of length 30.  allows user to
c                 specify some integer parameters which affect
c                 the method. 
c       rparm     real vector of length 30.  allows user to 
c                 specify some real parameters which affect 
c                 the method. 
c       ier       output integer.  error flag.
c
c ... specifications for parameters
c
      external  accel, precon 
      integer   iparm(30), jcoef(2), p(1), ip(1), iwksp(1)
      dimension coef(1), rhs(1), u(1), ubar(1), rparm(30), wksp(1)
c
c *** begin -- package common 
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      ier = 0
      ndim = ndimm
      mdim = mdimm
      n = nn
      maxnz = maxnzz
      lenr = nw
      leni = inw
      irmax = 0
      iimax = 0
      t1 = timer (dummy)
      call echall (n,iparm,rparm,1,1,ier)
      if (ier .lt. 0) return
      timfac = 0.0
      call pointr (1,wksp,iwksp,ier)
c
c ... call preparatory routines.
c
c ... remove zeros from jcoef for purdue data structure.
c
      if (nstore .eq. 1) call adjust (n,ndim,maxnz,jcoef,1) 
      call prep (coef,jcoef,wksp(irpnt),iwksp(iipnt),n,nstore,ier)
      if (ier .lt. 0) then
         call ershow (ier,'nspcg')
         go to 20
      endif
c
c ... eliminate penalty-method dirichlet points, if requested.
c
      ielim = iparm(24)
      tol = rparm(15)
      if (ielim .eq. 1) call elim (n,jcoef,coef,rhs,wksp,iwksp,
     a                             tol) 
c
c ... determine symmetry of matrix.
c
      if (nstore .eq. 1 .and. isymm .eq. 2) call detsym
     a                (ndim,maxnz,coef,jcoef,n,isymm)
c
c ... scale matrix. 
c
      call scale (coef,jcoef,wksp,1,n,u,ubar,rhs,ier)
      if (ier .lt. 0) go to 20
c
c ... permute matrix.
c
      call permut (coef,jcoef,p,ip,wksp,iwksp,1,n,u,ubar,rhs,ier)
      if (ier .lt. 0) go to 15
c
c ... call iterative routine. 
c
      call precon (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a             iparm,rparm,ier)
c
c ... unpermute matrix.
c
      call permut (coef,jcoef,p,ip,wksp,iwksp,2,n,u,ubar,rhs,ier)
c
c ... unscale matrix.
c
 15   call scale (coef,jcoef,wksp,2,n,u,ubar,rhs,ier)
c
c ... restore zeros to jcoef for purdue data structure.
c
 20   if (nstore .eq. 1) call adjust (n,ndim,maxnz,jcoef,2) 
      t2 = timer (dummy)
      timtot = t2 - t1
      iparm(18) = ipropa
      iparm(23) = isymm
      rparm(13) = timfac
      rparm(14) = timtot
      call echall (n,iparm,rparm,2,1,ier)
c
      call pointr (2,wksp,iwksp,ier)
      nw = irmax
      inw = iimax
      maxnzz = maxnz
      return
      end 
      subroutine rsnsp (precon,accel,ndimm,mdimm,nn,maxnzz,coef,
     a                  jcoef,p,ip,u,ubar,rhs,wksp,iwksp,nw,inw,
     a                  iparm,rparm,ier)
c
c ... rsnsp is the driver for the nspcg package for methods 
c     applied to the explicitly computed reduced system.
c
c ... parameters -- 
c
c       precon    preconditioning module
c       accel     acceleration module
c       coef      real matrix data array
c       jcoef     integer matrix data array
c       n         input integer.  order of the system (= nn)
c       u         input/output vector.  on input, u contains the
c                 initial guess to the solution.  on output, it
c                 contains the latest estimate to the solution.
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       wksp      real workspace of length nw
c       iwksp     integer workspace of length inw 
c       nw        length of wksp upon input, amount used upon
c                  output
c       inw       length of iwksp upon input, amount used upon
c                  output
c       iparm     integer vector of length 30.  allows user to
c                 specify some integer parameters which affect
c                 the method. 
c       rparm     real vector of length 30.  allows user to 
c                 specify some real parameters which affect 
c                 the method. 
c       ier       output integer.  error flag.
c
c ... specifications for parameters
c
      external  accel, precon 
      integer   iparm(30), jcoef(2), p(1), ip(1), iwksp(1)
      dimension coef(1), rhs(1), u(1), ubar(1), rparm(30), wksp(1)
c
c *** begin -- package common 
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      ier = 0
      ndim = ndimm
      mdim = mdimm
      n = nn
      maxnz = maxnzz
      lenr = nw
      leni = inw
      irmax = 0
      iimax = 0
      t1 = timer (dummy)
      call echall (n,iparm,rparm,1,1,ier)
      timfac = 0.0
      call pointr (1,wksp,iwksp,ier)
c
c ... call preparatory routines.
c
c ... remove zeros from jcoef for purdue data structure.
c
      if (nstore .eq. 1) call adjust (n,ndim,maxnz,jcoef,1) 
      call prep (coef,jcoef,wksp(irpnt),iwksp(iipnt),n,nstore,ier)
      if (ier .lt. 0) then
         call ershow (ier,'rsnsp')
         go to 20
      endif
c
c ... eliminate penalty-method dirichlet points, if requested.
c
      ielim = iparm(24)
      tol = rparm(15)
      if (ielim .eq. 1) call elim (n,jcoef,coef,rhs,wksp,iwksp,
     a                             tol) 
c
c ... determine symmetry of matrix.
c
      if (nstore .eq. 1 .and. isymm .eq. 2) call detsym
     a                (ndim,maxnz,coef,jcoef,n,isymm)
c
c ... form reduced system matrix.
c
      call rsprep (coef,jcoef,wksp,iwksp,n,rhs,u,ubar,
     a             p,ip,nr,irs,ijcrs,irsrhs,ier)
c
c ... scale matrix. 
c
      call scale (wksp(irs),iwksp(ijcrs),wksp,1,nr,u,ubar,
     a            wksp(irsrhs),ier)
      if (ier .lt. 0) go to 20
c
c ... call iterative routine. 
c
      call precon (accel,wksp(irs),iwksp(ijcrs),nr,u,ubar,
     a             wksp(irsrhs),wksp,iwksp,iparm,rparm,ier) 
c
c ... unscale matrix.
c
      call scale (wksp(irs),iwksp(ijcrs),wksp,2,nr,u,ubar,
     a            wksp(irsrhs),ier)
c
c ... restore to original system.
c
      call rspost (coef,jcoef,wksp,iwksp,n,rhs,u,ubar,
     a             p,ip,nr,irs,ijcrs,ier)
c
c ... restore zeros to jcoef for purdue data structure.
c
 20   if (nstore .eq. 1) call adjust (n,ndim,maxnz,jcoef,2) 
      t2 = timer (dummy)
      timtot = t2 - t1
      iparm(18) = ipropa
      iparm(23) = isymm
      rparm(13) = timfac
      rparm(14) = timtot
      call echall (n,iparm,rparm,2,1,ier)
c
      call pointr (2,wksp,iwksp,ier)
      nw = irmax
      inw = iimax
      maxnzz = maxnz
      return
      end 
      subroutine prep (coef,jcoef,wksp,iwksp,nn,nstore,ier) 
c
c ... prep puts the diagonal entries of the matrix into column
c     one of coef.
c
c ... parameters -- 
c
c         n       dimension of matrix
c         jcoef   integer matrix representation array
c         coef    matrix representation array
c         wksp    workspace array of size n
c         iwksp   integer workspace
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), iwksp(1)
      dimension coef(1), wksp(1)
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cmpart / mpstrt, mpart
c
      n = nn
      go to (5,10,10,15,15), nstore
 5    call prep1 (n,ndim,maxnz,jcoef,coef,ier)
      return
 10   call prep2 (n,ndim,maxnz,jcoef,coef,wksp,ier)
      return
 15   call needw ('prep',1,iipnt,2*n+1,ier)
      if (ier .lt. 0) return
      call prep3 (n,maxnz,jcoef,jcoef(ndim+1),coef,mpart,
     a            iwksp,iwksp(n+2))
      mpstrt = iipnt
      iipnt = iipnt + mpart + 1
      return
      end 
      subroutine pointr (icall,wksp,iwksp,ier)
c
c ... pointr adjusts pointers according to ifact. 
c
c ... parameters -- 
c
c       icall     indicates beginning or ending call
c                  = 1  for beginning
c                  = 2  for ending
c       wksp      real workspace vector 
c       iwksp     integer workspace vector
c
c ... specifications for parameters
c
      integer   iwksp(1)
      dimension wksp(1)
c
c *** begin -- package common 
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      if (icall .eq. 2) go to 15
c
c ... initialize pointers.
c
      if (ifact .eq. 0) return
      iipnt = 1
      irpnt = 1
      nfactr = 0
      nfacti = 0
      ifactr = 1
      ifacti = 1
      return
c
c ... reset pointers for return
c
 15   if (ier .lt. 0) return
      if (nfacti .eq. 0) go to 20
      call vicopy (nfacti,iwksp(ifacti),iwksp)
      iipnt = nfacti + 1
      ifacti = 1
 20   if (nfactr .eq. 0) return
      call vcopy (nfactr,wksp(ifactr),wksp)
      iwkpt2 = iwkpt2 - ifactr + 1
      irpnt = nfactr + 1
      ifactr = 1
      return
      end 
      subroutine needw (subnam,isw,istart,length,ier)
c
c ... needw determines if enough integer or real workspace
c     is available. 
c
c ... parameters -- 
c
c        subnam  name of calling routine
c        isw     switch for real or integer workspace check 
c                 = 0     real
c                 = 1     integer
c        istart  starting address
c        length  length desired
c        ier     error indicator (output)
c                 = -2    insufficient real workspace
c                 = -3    insufficient integer workspace
c
c ... specifications for parameters
c
      character*10 subnam
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
c
      newlen = istart + length - 1
      if (isw .eq. 1) go to 10
      if (lenr .ge. newlen) go to 5
      ier = -2
      call ershow (ier,subnam)
 5    irmax = max0 (irmax,newlen)
      return
 10   if (leni .ge. newlen) go to 15
      ier = -3
      call ershow (ier,subnam)
 15   iimax = max0 (iimax,newlen)
      return
      end 
      subroutine scale (coef,jcoef,wksp,icall,n,u,ubar,rhs,ier)
c
c ... scale scales the matrix, u, ubar, and rhs.
c
c ... parameters -- 
c
c       icall     key to indicate whether scaling (icall=1) or
c                  unscaling (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -4   nonpositive diagonal element
c
c ... specifications for parameters
c
      integer jcoef(2)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(1)
c
c ... data common blocks
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      if (iscale .ne. 1) return
      go to (5,10,10,15,15), nstore
 5    call scalep (coef,jcoef,wksp,icall,n,u,ubar,rhs,ier)
      return
 10   call scaled (coef,jcoef,wksp,icall,n,u,ubar,rhs,ier)
      return
 15   call scales (coef,jcoef,wksp,icall,n,u,ubar,rhs,ier)
      return
      end 
      subroutine permut (coef,jcoef,p,ip,wksp,iwksp,icall,n,u,
     a                   ubar,rhs,ier)
c
c ... permut permutes the matrix, u, ubar, and rhs.
c
c ... parameters -- 
c
c       icall     key to indicate whether permuting (icall=1) or
c                  unpermuting (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -2   insufficient real space to permute system
c                  =  -3   insufficient integer space to permute
c                          system
c
c ... specifications for parameters
c
      integer jcoef(2), p(1), ip(1), iwksp(1)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(1)
c
c ... data common blocks
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      if (iperm .ne. 1) return
      go to (5,10,10,15,15), nstore
 5    call permp (coef,jcoef,p,ip,wksp,iwksp,
     a            icall,n,u,ubar,rhs,ier)
      return
 10   call permd (coef,jcoef,p,ip,wksp,iwksp,
     a            icall,n,u,ubar,rhs,ier)
      return
 15   call perms (coef,jcoef,p,ip,wksp,iwksp,
     a            icall,n,u,ubar,rhs,ier)
      return
      end 
      subroutine elim (n,jcoef,coef,rhs,wksp,iwksp,toll)
c
c ... elim removes rows of the matrix for which the ratio of the
c     sum of off-diagonal elements to the diagonal element is
c     small (less than tol) in absolute value.
c     this is to take care of matrices arising from finite
c     element discretizations of partial differential equations
c     with dirichlet boundary conditions implemented by penalty
c     methods.  any such rows and corresponding columns are then
c     eliminated (set to the identity after correcting the rhs).
c
c ... parameter list --
c
c         n       dimension of matrix
c         jcoef   integer array of matrix representation
c         coef    array of sparse matrix representation
c         rhs     right hand side of matrix problem
c         wksp    wksp array of length n
c         tol     tolerance factor  (= toll)
c
c ... specifications for arguments
c
      common / cmpart / mpstrt, mpart
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      integer   jcoef(2), iwksp(1)
      dimension coef(1), rhs(1), wksp(1)
c
      tol = toll
      go to (5,10,15,20,25), nstore
 5    call elim1 (n,ndim,maxnz,jcoef,coef,rhs,wksp(irpnt),tol)
      return
 10   call elim2 (n,ndim,maxnz,jcoef,coef,rhs,wksp(irpnt),tol)
      return
 15   call elim3 (n,ndim,maxnz,jcoef,coef,rhs,wksp(irpnt),tol)
      return
 20   call elim4 (mpart,iwksp(mpstrt),jcoef,jcoef(ndim+1),
     a            coef,rhs,wksp(irpnt),tol)
      return
 25   call elim5 (mpart,iwksp(mpstrt),jcoef,jcoef(ndim+1),
     a            coef,rhs,wksp(irpnt),tol)
      return
      end 
      subroutine scaled (coef,jcoef,wksp,icall,nn,u,ubar,rhs,ier)
c
c ... scaled scales the matrix, u, ubar, and rhs. 
c     (symmetric or nonsymmetric diagonal format) 
c
c ... parameters -- 
c
c       icall     key to indicate whether scaling (icall=1) or
c                  unscaling (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -4   nonpositive diagonal element
c
c ... specifications for parameters
c
      integer jcoef(2)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      iflag = 0
      if (ntest .eq. 6) iflag = 1
      if (icall .eq. 2) go to 20
c
c ... scale system. 
c
c ... check for sufficient room.
c
      call needw ('scaled',0,irpnt,n,ier)
      if (ier .lt. 0) return
      iptscl = irpnt
      irpnt = irpnt + n
      call scal2 (n,ndim,maxnz,jcoef,coef,rhs,u,ubar,
     a            wksp(iptscl),iflag,ier)
      if (ier .lt. 0) call ershow (ier,'scaled')
      return
c
c ... unscale system.
c
 20   call uscal2 (n,ndim,maxnz,jcoef,coef,rhs,u,ubar,
     a             wksp(iptscl),iflag)
      return
      end 
      subroutine scalep (coef,jcoef,wksp,icall,nn,u,ubar,rhs,ier)
c
c ... scalep scales the matrix, u, ubar, and rhs. 
c     (purdue format)
c
c ... parameters -- 
c
c       icall     key to indicate whether scaling (icall=1) or
c                  unscaling (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -4   nonpositive diagonal element
c
c ... specifications for parameters
c
      integer jcoef(2)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      iflag = 0
      if (ntest .eq. 6) iflag = 1
      if (icall .eq. 2) go to 20
c
c ... scale system. 
c
c ... check for sufficient room.
c
      call needw ('scalep',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iptscl = irpnt
      irpnt = irpnt + n
      call scal1 (n,ndim,maxnz,jcoef,coef,rhs,u,ubar,
     a            wksp(iptscl),wksp(irpnt),iflag,ier)
      if (ier .lt. 0) call ershow (ier,'scalep')
      return
c
c ... unscale system.
c
 20   call uscal1 (n,ndim,maxnz,jcoef,coef,rhs,u,ubar,
     a             wksp(iptscl),wksp(irpnt),iflag)
      return
      end 
      subroutine scales (coef,jcoef,wksp,icall,nn,u,ubar,rhs,ier)
c
c ... scales scales the matrix, u, ubar, and rhs. 
c     (sparse format)
c
c ... parameters -- 
c
c       icall     key to indicate whether scaling (icall=1) or
c                  unscaling (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -4   nonpositive diagonal element
c
c ... specifications for parameters
c
      integer jcoef(2)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      iflag = 0
      if (ntest .eq. 6) iflag = 1
      if (icall .eq. 2) go to 10
c
c ... scale system. 
c
c ... check for sufficient room.
c
      call needw ('scales',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iptscl = irpnt
      irpnt = irpnt + n
      call scal3 (n,maxnz,jcoef,jcoef(ndim+1),coef,rhs,u,ubar,
     a            wksp(iptscl),wksp(irpnt),iflag,ier)
      if (ier .lt. 0) call ershow (ier,'scales')
      return
c
c ... unscale system.
c
 10   call uscal3 (n,maxnz,jcoef,jcoef(ndim+1),coef,rhs,u,ubar,
     a             wksp(iptscl),wksp(irpnt),iflag)
      return
      end 
      subroutine permd (coef,jcoef,p,ip,wksp,iwksp,icall,nn,u,
     a                  ubar,rhs,ier)
c
c ... permd permutes the matrix, u, ubar, and rhs.
c     (diagonal format)
c
c ... parameters -- 
c
c       icall     key to indicate whether permuting (icall=1) or
c                  unpermuting (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -2   insufficient real space to permute system
c                  =  -3   insufficient integer space to permute
c                          system
c
c ... specifications for parameters
c
      integer jcoef(2), p(1), ip(1), iwksp(1)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
c
      n = nn
      if (icall .eq. 2) go to 65
c
c ... permute system.
c
c ... check for sufficient storage to permute matrix
c
      call needw ('permd',0,irpnt,n,ier)
      if (ier .lt. 0) return
      nc = iipnt
      call needw ('permd',1,nc,n,ier)
      if (ier .lt. 0) return
      call pgen (n,p,ip,iwksp(nc),ncolor)
      ipt = nc + ncolor
      ncmax = 0
      do 16 i = nc,ipt-1
         if (ncmax .lt. iwksp(i)) ncmax = iwksp(i)
 16   continue
      call needw ('permd',1,ipt,ncolor+1,ier)
      if (ier .lt. 0) return
      call iptgen (ncolor,iwksp(ipt),iwksp(nc))
      maxnew = ipt + ncolor + 1
      jcnew = maxnew + ncolor 
      lbhb = jcnew + ncolor*mdim
      call needw ('permd',1,maxnew,ncolor+ncolor*mdim+n,ier)
      if (ier .lt. 0) return
      isym = nstore - 2
      call pmdg (ndim,mdim,n,maxnz,jcoef,coef,ncolor,iwksp(nc),
     a           p,ip,maxd,iwksp(maxnew),
     a           iwksp(jcnew),wksp(irpnt),iwksp(lbhb),isym,ier)
      if (ier .lt. 0) then
         call ershow (ier,'permd')
         return
      endif
      lbhb = jcnew + ncolor*maxd
      iblock = lbhb + ncolor
      call move4 (ndim,n,iwksp(maxnew),iwksp(jcnew),coef,ncolor,
     a            iwksp(nc),wksp(irpnt),iwksp(lbhb))
      call needw ('permd',1,lbhb,ncolor+3*ncolor*(maxd+1),ier)
      if (ier .lt. 0) return
      call define (ndim,iwksp(maxnew),iwksp(jcnew),coef,ncolor,
     a             iwksp(nc),iwksp(iblock),iwksp(lbhb))
      lbhbm = iwksp(lbhb)
      do 45 j = 1,ncolor-1
         if (lbhbm .lt. iwksp(lbhb+j)) lbhbm = iwksp(lbhb+j)
 45   continue
      is1 = iblock + 3*ncolor*lbhbm
      is2 = is1 + ncolor
      call needw ('permd',1,is1,2*ncolor,ier)
      if (ier .lt. 0) return
      call prbblk (ncolor,ncolor,iwksp(iblock),iwksp(lbhb), 
     a             iwksp(is1),iwksp(is2),propa)
      if (propa) ipropa = 1
      if (.not. propa) ipropa = 0
      iipnt = is1
      call pervec (n,p,rhs,wksp(irpnt)) 
      call pervec (n,p,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,p,ubar,wksp(irpnt))
      return
c
c ... unpermute system.
c
 65   call needw ('permd',1,iipnt,2*n,ier)
      if (ier .lt. 0) return
      isym = nstore - 2
      call unpmdg (ndim,n,maxnz,jcoef,coef,ncolor,iwksp(nc),
     a             p,ip,maxd,iwksp(maxnew),
     a             iwksp(jcnew),wksp(irpnt),iwksp(iipnt),isym)
      call pervec (n,ip,rhs,wksp(irpnt))
      call pervec (n,ip,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,ip,ubar,wksp(irpnt)) 
      return
      end 
      subroutine permp (coef,jcoef,p,ip,wksp,iwksp,icall,nn,u,
     a                  ubar,rhs,ier)
c
c ... permp permutes the matrix, u, ubar, and rhs.
c     (purdue format)
c
c ... parameters -- 
c
c       icall     key to indicate whether permuting (icall=1) or
c                  unpermuting (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -2   insufficient real space to permute system
c                  =  -3   insufficient integer space to permute
c                          system
c
c ... specifications for parameters
c
      integer jcoef(2), p(1), ip(1), iwksp(1)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
c
      n = nn
      if (icall .eq. 2) go to 40
c
c ... permute system.
c
c ... check for sufficient storage to permute matrix
c
      call needw ('permp',0,irpnt,n,ier)
      if (ier .lt. 0) return
      nc = iipnt
      call needw ('permp',1,nc,n,ier)
      if (ier .lt. 0) return
      call pgen (n,p,ip,iwksp(nc),ncolor)
      ipt = nc + ncolor
      ncmax = 0
      do 20 i = nc,ipt-1
         if (ncmax .lt. iwksp(i)) ncmax = iwksp(i)
 20   continue
      call needw ('permp',1,ipt,ncolor+1,ier)
      if (ier .lt. 0) return
      call iptgen (ncolor,iwksp(ipt),iwksp(nc))
      iipnt = iipnt + 2*ncolor + 1
      call needw ('permp',1,iipnt,n,ier)
      if (ier .lt. 0) return
      call permat (ndim,maxnz,coef,jcoef,
     a                    wksp(irpnt),iwksp(iipnt),n,p)
      call needw ('permp',1,iipnt,2*ncolor,ier)
      if (ier .lt. 0) return
      ndt = iipnt
      ndb = iipnt + ncolor
      call move3 (ndim,mdim,n,maxnz,jcoef,coef,iwksp(ndt),
     a            iwksp(ndb),ncolor,iwksp(nc),ier)
      iipnt = iipnt + 2*ncolor
      if (ier .lt. 0) then
         call ershow (ier,'permp')
         return
      endif
      call pervec (n,p,rhs,wksp(irpnt)) 
      call pervec (n,p,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,p,ubar,wksp(irpnt))
      return
c
c ... unpermute system.
c
 40   call needw ('permp',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call needw ('permp',1,iipnt,n,ier)
      if (ier .lt. 0) return
      call permat (ndim,maxnz,coef,jcoef,
     a                    wksp(irpnt),iwksp(iipnt),n,ip)
      call pervec (n,ip,rhs,wksp(irpnt))
      call pervec (n,ip,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,ip,ubar,wksp(irpnt)) 
      return
      end 
      subroutine perms (coef,jcoef,p,ip,wksp,iwksp,icall,nn,u,
     a                  ubar,rhs,ier)
c
c ... perms permutes the matrix, u, ubar, and rhs.
c     (sparse format)
c
c ... parameters -- 
c
c       icall     key to indicate whether permuting (icall=1) or
c                  unpermuting (icall=2) is to be done
c       n         order of system
c       u         current solution estimate
c       ubar      input vector containing the true solution 
c                  (optional) 
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       ier       error flag
c                  =   0   no errors detected
c                  =  -2   insufficient real space to permute system
c                  =  -3   insufficient integer space to permute
c                          system
c
c ... specifications for parameters
c
      integer jcoef(2), p(1), ip(1), iwksp(1)
      dimension rhs(1), u(1), ubar(1), coef(1), wksp(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
c
c *** end   -- package common 
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      isym = 0
      if (nstore .eq. 5) isym = 1
      if (icall .eq. 2) go to 10
c
c ... permute system.
c
c ... check for sufficient storage to permute matrix
c
      call needw ('perms',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call needw ('perms',1,iipnt,n,ier)
      if (ier .lt. 0) return
      call pgen (n,p,ip,iwksp(iipnt),ncolor)
      call permas (isym,n,maxnz,jcoef,jcoef(ndim+1),
     a             coef,wksp(irpnt),p)
      call pervec (n,p,rhs,wksp(irpnt)) 
      call pervec (n,p,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,p,ubar,wksp(irpnt))
      return
c
c ... unpermute system.
c
 10   call needw ('perms',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call needw ('perms',1,iipnt,n,ier)
      if (ier .lt. 0) return
      call permas (isym,n,maxnz,jcoef,jcoef(ndim+1),
     a             coef,wksp(irpnt),ip) 
      call pervec (n,ip,rhs,wksp(irpnt))
      call pervec (n,ip,u,wksp(irpnt))
      if (ntest .eq. 6) call pervec (n,ip,ubar,wksp(irpnt)) 
      return
      end 
      subroutine rsprep (coef,jcoef,wksp,iwksp,nn,rhs,u,ubar,
     a                   p,ip,nrr,irs,ijcrs,irsrhs,ier)
c
c ... rsprep is the preprocessor for methods using the
c     explicitly-computed reduced system.
c
c ... parameters -- 
c
c       coef      real matrix data array
c       jcoef     integer matrix data array
c       n         input integer.  order of the system (= nn)
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       u         current solution estimate
c       ubar      exact solution vector (if known)
c       nr        order of the reduced system upon output
c       irs       pointer into wksp for reduced system matrix
c       ijcrs     pointer into wksp for reduced system integer
c                  array
c       irsrhs    pointer into wksp for reduced system rhs
c       ier       output integer.  error flag.
c
c ... specifications for parameters
c
      integer jcoef(2), iwksp(1), p(1), ip(1)
      dimension coef(1), rhs(1), u(1), ubar(1), wksp(1)
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / rscons / ndimrs, mdimrs, maxzrs
c
      n = nn
c
c ... permute matrix.
c
      call permut (coef,jcoef,p,ip,wksp,iwksp,1,
     a             n,u,ubar,rhs,ier)
      if (ier .lt. 0) return
c
c ... form reduced system matrix.
c
      nr = iwksp(nc)
      nb = iwksp(nc+1)
      irs = irpnt
      ijcrs = iipnt 
      length = lenr - irpnt + 1
      call vfill (length,wksp(irpnt),0.0)
      if (nstore .ge. 2) go to 30
c
c ... purdue storage.
c
      call needw ('rsprep',0,irpnt,3*nr,ier)
      if (ier .lt. 0) return
      call needw ('rsprep',1,iipnt,2*nr,ier)
      if (ier .lt. 0) return
      lim1 = (lenr - 2*nr - irpnt + 1)/nr
      lim2 = (leni - nr - iipnt + 1)/nr 
      maxlim = min0(lim1,lim2)
      ip1 = irpnt + nr*maxlim 
      ip2 = iipnt + nr*maxlim 
      call rsmatp (ndim,nr,maxnz,jcoef,coef,maxrs,iwksp(ijcrs),
     a             wksp(irs),maxlim,wksp(ip1),iwksp(ip2),ier)
      if (ier .lt. 0) then
         call ershow (ier,'rsprep')
         return
      endif
      irpnt = irpnt + nr*maxrs
      iipnt = iipnt + nr*maxrs
      go to 45
c
c ... diagonal storage.
c
 30   call needw ('rsprep',0,irpnt,nr,ier)
      if (ier .lt. 0) return
      call needw ('rsprep',1,iipnt,nr,ier)
      if (ier .lt. 0) return
      maxlim = length/nr
      isym = 0
      if (nstore .eq. 3) isym = 1
      call rsmatd (ndim,nr,nb,iwksp(maxnew),iwksp(jcnew),coef,
     a             coef(ndim+1),coef(nr+ndim+1),coef(nr+1),maxrs,
     a             iwksp(ijcrs),wksp(irs),maxlim,isym,ier)
      if (ier .lt. 0) then
         call ershow (ier,'rsprep')
         return
      endif
      irpnt = irpnt + nr*maxrs
      iipnt = iipnt + maxrs
c
c ... form reduced system rhs.
c
 45   irsrhs = irpnt
      ip1 = irpnt + nr
      call needw ('rsprep',0,irpnt,n+nr,ier)
      if (ier .lt. 0) return
      if (nstore .eq. 1) call rsbegp (n,nr,ndim,maxnz,jcoef,
     a                     coef,wksp(irsrhs),rhs,wksp(ip1)) 
      if (nstore .ge. 2) call rsrhsd (n,nr,ndim,iwksp(maxnew),
     a                    iwksp(jcnew),coef,wksp(irsrhs),rhs,
     a                    wksp(ip1))
      irpnt = irpnt + nr
c
c ... update constants.
c
      ndimrs = ndim 
      mdimrs = mdim 
      maxzrs = maxnz
      ndim = nr
      mdim = maxrs
      maxnz = maxrs 
      nrr = nr
      return
      end 
      subroutine rspost (coef,jcoef,wksp,iwksp,nn,rhs,u,ubar,
     a                   p,ip,nrr,irs,ijcrs,ier)
c
c ... rspost is the postprocessor for methods using the
c     explicitly-computed reduced system.
c
c ... parameters -- 
c
c       coef      real matrix data array
c       jcoef     integer matrix data array
c       n         input integer.  order of the system (= nn)
c       rhs       input vector.  contains the right hand side
c                 of the matrix problem.
c       u         current solution estimate
c       ubar      exact solution vector (if known)
c       nr        order of the reduced system upon input
c       irs       pointer into wksp for reduced system matrix
c       ijcrs     pointer into wksp for reduced system integer
c                  array
c       ier       output integer.  error flag.
c
c ... specifications for parameters
c
      integer jcoef(2), iwksp(1), p(1), ip(1)
      dimension coef(1), rhs(1), u(1), ubar(1), wksp(1)
c
c ... data common blocks
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / rscons / ndimrs, mdimrs, maxzrs
c
      n = nn
      nr = nrr
      nb = n - nr
c
c ... update constants.
c
      ndim = ndimrs 
      mdim = mdimrs 
      maxnz = maxzrs
      irpnt = irs
      iipnt = ijcrs 
c
c ... compute xb.
c
      call needw ('rspost',0,irpnt,nb,ier)
      if (ier .lt. 0) return
      if (nstore .eq. 1) call rsendp (n,nr,ndim,maxnz,jcoef,
     a                     coef,u,rhs,wksp(irpnt))
      if (nstore .ge. 2) call rsxbd (n,nr,ndim,iwksp(maxnew),
     a                    iwksp(jcnew),coef,u,rhs)
c
c ... unpermute matrix.
c
      call permut (coef,jcoef,p,ip,wksp,iwksp,2,
     a             n,u,ubar,rhs,ier)
      if (ier .lt. 0) return
      return
      end 
      subroutine redblk (ndim,n,maxnz,coef,jcoef,p,ip,nstore,
     a                   iwksp,ier)
c
c ... redblk determines if the matrix has property a.
c
c ... parameters -- 
c
c        n        problem size
c        nstore   storage mode
c                  = 1  purdue format
c                  = 2  symmetric diagonal format 
c                  = 3  nonsymmetric diagonal format
c                  = 4  symmetric sparse format
c                  = 5  nonsymmetric sparse format
c        iwksp    integer workspace vector of length n
c        ier      error code
c                  =  0   no errors detected
c                  = -8   matrix does not have property a
c
c ... common blocks 
c
      integer jcoef(2), p(1), ip(1), iwksp(1)
      dimension coef(1)
      logical           propal
c
      go to (5,5,5,10,10), nstore
 5    call prbndx (n,ndim,maxnz,jcoef,coef,p,ip,propal,nstore)
      go to 15
 10   call bicol (n,maxnz,jcoef,jcoef(ndim+1),p,ip,iwksp,propal)
 15   if (propal) ier = 0
      if (.not. propal) ier = -8
      if (propal) return
      call ershow (ier,'redblk')
      return
      end 
      subroutine noadp (coef,jcoef,wksp,iwksp,n,p,r,pdp,pldup)
c
c ... noadp is a dummy routine to do no adaption. 
c
c ... specifications for parameters
c
c
      integer jcoef(2), iwksp(1)
      dimension p(1), r(1), coef(1), wksp(1)
c
      return
      end 
      subroutine copy (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... copy does a vector copy (null preconditioner)
c
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      do 10 i = 1,n 
 10   z(i) = r(i)
      return
      end 
      subroutine split (accel,suba,subat,subq,subqt,subql,subqlt,
     a                  subqr,subqrt,subadp,
     a                  coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... split determines how to apply the splitting based on
c     iqlr.
c
      external accel, suba, subat, subq, subqt, subql, subqlt,
     a         subqr, subqrt, subadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(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
c
c *** end   -- package common 
c
      if (iqlr .eq. 0) then
         call accel (suba,subat,copy,copy,copy,copy,subadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a            iparm,rparm,jer)
      endif
      if (iqlr .eq. 1) then
         call accel (suba,subat,subq,subqt,copy,copy,subadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a            iparm,rparm,jer)
      endif
      if (iqlr .eq. 2) then
         call accel (suba,subat,copy,copy,subq,subqt,subadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a            iparm,rparm,jer)
      endif
      if (iqlr .eq. 3) then
         call accel (suba,subat,subql,subqlt,subqr,subqrt,
     a            subadp,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a            iparm,rparm,jer)
      endif
      if (jer .ne. 0) ier = jer
      return
      end 
      subroutine rich1 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... rich1 drives the richardson preconditioner. 
c
      external accel, suba8, suba9, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + n
      call split (accel,suba8,suba9,copy,copy,copy,copy,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,
     a            ier)
      if (keygs .eq. 1) irpnt = irpnt - n
      return
      end 
      subroutine jac1 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... jac1 drives the jacobi preconditioner.
c
      external accel, suba8, suba9, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + n
      call split (accel,suba8,suba9,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (keygs .eq. 1) irpnt = irpnt - n
      return
      end 
      subroutine sor1 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... sor1 drives the point sor method. 
c
      external accel, suba8, suba9, subq78, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call move1 (ndim,mdim,n,maxnz,jcoef,coef,maxt,maxb,ier)
      if (ier .lt. 0) then
         call ershow (ier,'sor1')
         return
      endif
      call split (accel,suba8,suba9,subq78,subq78,subq78,subq78,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine ssor1 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ssor1 drives the point ssor method.
c
      external accel, suba8, suba9, subq79, subq80, subq81, subq82,
     a         subq83, subq84, subq85
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      iwkpt1 = irpnt
      irpnt = irpnt + n
      if (isymm .ne. 0) irpnt = irpnt + n
      call move1 (ndim,mdim,n,maxnz,jcoef,coef,maxt,maxb,ier)
      if (ier .lt. 0) then
         call ershow (ier,'ssor1')
         return
      endif
      call split (accel,suba8,suba9,subq79,subq80,subq81,subq82,
     a            subq83,subq84,subq85, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (isymm .ne. 0) irpnt = irpnt - n
      irpnt = irpnt - n
      return
      end 
      subroutine ic1 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... ic1 drives the ic preconditioner. 
c
      external accel, suba8, suba9, subq86, subq87, subq88, 
     a         subq89, subq90, subq91, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c
      n = nn
      if (ifact .eq. 0 .and. lvfill .gt. 0) go to 20
      call move1 (ndim,mdim,n,maxnz,jcoef,coef,maxt,maxb,ier)
      if (ier .lt. 0) then
         call ershow (ier,'ic1') 
         return
      endif
 20   t1 = timer (dummy)
      if (ifact .eq. 1) call pfact1 (coef,jcoef,wksp,iwksp,n,1,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba8,suba9,subq86,subq87,subq88,subq89,
     a            subq90,subq91,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine mic1 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... mic1 drives the mic preconditioner.
c
      external accel, suba8, suba9, subq86, subq87, subq88, 
     a         subq89, subq90, subq91, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c
      n = nn
      if (ifact .eq. 0 .and. lvfill .gt. 0) go to 20
      call move1 (ndim,mdim,n,maxnz,jcoef,coef,maxt,maxb,ier)
      if (ier .lt. 0) then
         call ershow (ier,'mic1')
         return
      endif
 20   t1 = timer (dummy)
      if (ifact .eq. 1) call pfact1 (coef,jcoef,wksp,iwksp,n,2,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba8,suba9,subq86,subq87,subq88,subq89,
     a            subq90,subq91,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine lsp1 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... lsp1 drives the least squares polynomial preconditioner.
c
      external accel, suba8, suba9, subq92, subq93, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('lsp1',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      call ainfn (n,ndim,maxnz,jcoef,coef,1,ainf,wksp(irpnt))
      iwkpt2 = irpnt
      irpnt = irpnt + 2*n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + n
      call split (accel,suba8,suba9,subq92,subq93,subq92,subq93,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      if (keygs .eq. 1) irpnt = irpnt - n
      return
      end 
      subroutine neu1 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... neu1 drives the neumann polynomial preconditioner.
c
      external accel, suba8, suba9, subq94, subq95, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('neu1',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      iwkpt2 = irpnt
      irpnt = irpnt + n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + n
      call split (accel,suba8,suba9,subq94,subq95,subq94,subq95,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      if (keygs .eq. 1) irpnt = irpnt - n
      return
      end 
      subroutine rich2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... rich2 drives the richardson preconditioner. 
c
      external accel, suba1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      call split (accel,suba1,suba1,copy,copy,copy,copy,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine jac2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... jac2 drives the jacobi preconditioner.
c
      external accel, suba1, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      call split (accel,suba1,suba1,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine ljac2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ljac2 drives the line jacobi preconditioner.
c
      external accel, suba1, subq2, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq2,subq2,subq2,subq2,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine ljacx2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... ljacx2 drives the line jacobi preconditioner.
c
      external accel, suba1, subq4, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call linv (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq4,subq4,subq4,subq4,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine sor2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... sor2 drives the point sor method. 
c
      external accel, suba1, subq6, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call rowise (maxnz,jcoef,irwise)
      call needw ('sor2',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      iwkpt1 = iipnt
      iipnt = iipnt + maxnz
      call split (accel,suba1,suba1,subq6,subq6,subq6,subq6,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - maxnz
      return
      end 
      subroutine ssor2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ssor2 drives the point ssor method.
c
      external accel, suba1, subq7, subq8, subq9, subq10,
     a         subq11, subq12 
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call rowise (maxnz,jcoef,irwise)
      call needw ('ssor2',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      iwkpt1 = iipnt
      iipnt = iipnt + maxnz
      call split (accel,suba1,suba1,subq7,subq7,subq8,subq9,
     a            subq10,subq11,subq12, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - maxnz
      return
      end 
      subroutine ic2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... ic2 drives the symmetric ic preconditioner. 
c
      external accel, suba1, subq13, subq14, subq15, subq16,
     a         subq17, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call pfact2 (coef,jcoef,wksp,iwksp,n,1,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      leniw = max0 (maxnz,nfacti)
      iwkpt1 = iipnt
      iipnt = iipnt + leniw
      call split (accel,suba1,suba1,subq13,subq13,subq14,subq15,
     a            subq16,subq17,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - leniw
      return
      end 
      subroutine mic2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... mic2 drives the symmetric mic preconditioner.
c
      external accel, suba1, subq13, subq14, subq15, subq16,
     a         subq17, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / dscons / ndim, mdim, maxnz
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call pfact2 (coef,jcoef,wksp,iwksp,n,2,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      leniw = max0 (maxnz,nfacti)
      iwkpt1 = iipnt
      iipnt = iipnt + leniw
      call split (accel,suba1,suba1,subq13,subq13,subq14,subq15,
     a            subq16,subq17,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - leniw
      return
      end 
      subroutine lsp2 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... lsp2 drives the least squares polynomial preconditioner.
c
      external accel, suba1, subq18, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      call needw ('lsp2',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      call ainfn (n,ndim,maxnz,jcoef,coef,2,ainf,wksp(irpnt))
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba1,suba1,subq18,subq18,subq18,subq18,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine neu2 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... neu2 drives the neumann polynomial preconditioner.
c
      external accel, suba1, subq19, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      call needw ('neu2',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba1,suba1,subq19,subq19,subq19,subq19,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine lsor2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... lsor2 drives the line sor method. 
c
      external accel, suba1, subq20, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq20,subq20,subq20,subq20,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine lssor2 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... lssor2 drives the line ssor method.
c
      external accel, suba1, subq21, subq22, copy 
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      n = nn
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + n
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq21,subq21,subq21,subq21,
     a            copy,copy,subq22,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine bic2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bic2 drives the block factorization (version 1) method.
c
      external accel, suba1, subq25, copy, noadp
      external ibfcs1
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacs (1,ibfcs1,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + kblsz
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq25,subq25,subq25,subq25,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - kblsz
      return
      end 
      subroutine mbic2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbic2 drives the block factorization (version 1, modified)
c     method.
c
      external accel, suba1, subq25, copy, noadp
      external ibfcs3
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacs (2,ibfcs3,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + kblsz
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq25,subq25,subq25,subq25,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - kblsz
      return
      end 
      subroutine bicx2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bicx2 drives the block factorization (version 2) method.
c
      external accel, suba1, subq25, copy, noadp
      external ibfcs2
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacs (3,ibfcs2,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + kblsz
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq25,subq25,subq25,subq25,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - kblsz
      return
      end 
      subroutine mbicx2 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbicx2 drives the block factorization (version 2, modified)
c     method.
c
      external accel, suba1, subq25, copy, noadp
      external ibfcs4
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacs (4,ibfcs4,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + kblsz
      if (ier .lt. 0) return
      call split (accel,suba1,suba1,subq25,subq25,subq25,subq25,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - kblsz
      return
      end 
      subroutine llsp2 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... llsp2 drives the line least squares polynomial preconditioner.
c
      external accel, suba1, subq23, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      call needw ('llsp2',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call adinfn (n,ndim,maxnz,jcoef,coef,2,ainf,wksp(irpnt))
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call needw ('llsp2',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba1,suba1,subq23,subq23,subq23,subq23,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine lneu2 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... lneu2 drives the line neumann polynomial preconditioner.
c
      external accel, suba1, subq24, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call needw ('lneu2',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba1,suba1,subq24,subq24,subq24,subq24,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine rich3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... rich3 drives the richardson preconditioner. 
c
      external accel, suba4, suba5, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      call split (accel,suba4,suba5,copy,copy,copy,copy,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine jac3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... jac3 drives the jacobi preconditioner.
c
      external accel, suba4, suba5, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      call split (accel,suba4,suba5,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine ljac3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ljac3 drives the line jacobi preconditioner.
c
      external accel, suba4, suba5, subq2, subq3, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba4,suba5,subq2,subq3,subq2,subq3,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine ljacx3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... ljacx3 drives the line jacobi preconditioner.
c
      external accel, suba4, suba5, subq4, subq5, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call linv (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba4,suba5,subq4,subq5,subq4,subq5,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine sor3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... sor3 drives the point sor method. 
c
      external accel, suba4, suba5, subq40, noadp, copy
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call rowise (maxnz,jcoef,irwise)
      call needw ('sor3',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      call needw ('sor3',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      call move2 (ndim,n,maxnz,jcoef,coef,wksp(irpnt),
     a            iwksp(iipnt),maxt,maxb)
      iwkpt1 = iipnt
      iipnt = iipnt + maxnz
      call split (accel,suba4,suba5,subq40,subq40,subq40,subq40,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - maxnz
      return
      end 
      subroutine ssor3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ssor3 drives the point ssor method.
c
      external accel, suba4, suba5, subq41, subq42, subq43, subq44,
     a         subq45, subq46, subq47
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      call rowise (maxnz,jcoef,irwise)
      call needw ('ssor3',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      call needw ('ssor3',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call move2 (ndim,n,maxnz,jcoef,coef,wksp(irpnt),
     a            iwksp(iipnt),maxt,maxb)
      iwkpt1 = irpnt
      irpnt = irpnt + n
      iwkpt2 = iipnt
      iipnt = iipnt + maxnz
      call split (accel,suba4,suba5,subq41,subq42,subq43,subq44,
     a            subq45,subq46,subq47, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      iipnt = iipnt - maxnz
      return
      end 
      subroutine ic3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... ic3 drives the nonsymmetric ic preconditioner.
c
      external accel, suba4, suba5, subq48, subq49, subq50, 
     a         subq51, subq52, subq53, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / dscons / ndim, mdim, maxnz
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
c
      n = nn
      call needw ('ic3',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      call needw ('ic3',0,irpnt,n,ier)
      if (ier .lt. 0) return
      if (ifact .eq. 0 .and. lvfill .gt. 0) go to 20
      call move2 (ndim,n,maxnz,jcoef,coef,wksp(irpnt),
     a            iwksp(iipnt),maxt,maxb)
 20   t1 = timer (dummy)
      if (ifact .eq. 1) call pfact3 (coef,jcoef,wksp,iwksp,n,1,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      leniw = max0 (maxnz,nfacti)
      iwkpt1 = iipnt
      iipnt = iipnt + leniw
      call split (accel,suba4,suba5,subq48,subq49,subq50,subq51,
     a            subq52,subq53,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - leniw
      return
      end 
      subroutine mic3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... mic3 drives the nonsymmetric mic preconditioner.
c
      external accel, suba4, suba5, subq48, subq49, subq50, 
     a         subq51, subq52, subq53, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / dscons / ndim, mdim, maxnz
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
c
      n = nn
      call needw ('mic3',1,iipnt,maxnz,ier)
      if (ier .lt. 0) return
      call needw ('mic3',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      if (ifact .eq. 0 .and. lvfill .gt. 0) go to 20
      call move2 (ndim,n,maxnz,jcoef,coef,wksp(irpnt),
     a            iwksp(iipnt),maxt,maxb)
 20   t1 = timer (dummy)
      if (ifact .eq. 1) call pfact3 (coef,jcoef,wksp,iwksp,n,2,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      leniw = max0 (maxnz,nfacti)
      iwkpt1 = iipnt
      iipnt = iipnt + leniw
      call split (accel,suba4,suba5,subq48,subq49,subq50,subq51,
     a            subq52,subq53,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      iipnt = iipnt - leniw
      return
      end 
      subroutine lsp3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... lsp3 drives the least squares polynomial preconditioner.
c
      external accel, suba4, suba5, subq54, subq55, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      call needw ('lsp3',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      call ainfn (n,ndim,maxnz,jcoef,coef,3,ainf,wksp(irpnt))
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba4,suba5,subq54,subq55,subq54,subq55,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine neu3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... neu3 drives the neumann polynomial preconditioner.
c
      external accel, suba4, suba5, subq56, subq57, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      call needw ('neu3',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba4,suba5,subq56,subq57,subq56,subq57,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine lsor3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... lsor3 drives the line sor method. 
c
      external accel, suba4, suba5, subq58, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba4,suba5,subq58,subq58,subq58,subq58,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine lssor3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... lssor3 drives the line ssor method.
c
      external accel, suba4, suba5, subq59, subq60, subq61, 
     a         subq62, subq63, subq64, subq65
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      n = nn
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      iwkpt1 = irpnt
      irpnt = irpnt + n
      if (ier .lt. 0) return
      call split (accel,suba4,suba5,subq59,subq60,subq61,subq62,
     a            subq63,subq64,subq65, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine bic3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bic3 drives the block factorization (version 1) method.
c
      external accel, suba4, suba5, subq70, subq71, subq72, 
     a         subq73, subq74, subq75, noadp
      external ibfcn1
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmz (1,ibfcn1,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*kblsz 
      call split (accel,suba4,suba5,subq70,subq71,subq72,subq73,
     a            subq74,subq75,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*kblsz 
      return
      end 
      subroutine mbic3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbic3 drives the block factorization (version 1, modified)
c     method.
c
      external accel, suba4, suba5, subq70, subq71, subq72, 
     a         subq73, subq74, subq75, noadp
      external ibfcn3
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmz (2,ibfcn3,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*kblsz 
      call split (accel,suba4,suba5,subq70,subq71,subq72,subq73,
     a            subq74,subq75,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*kblsz 
      return
      end 
      subroutine bicx3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bicx3 drives the block factorization (version 2)
c     method.
c
      external accel, suba4, suba5, subq70, subq71, subq72, 
     a         subq73, subq74, subq75, noadp
      external ibfcn2
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmz (3,ibfcn2,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*kblsz 
      call split (accel,suba4,suba5,subq70,subq71,subq72,subq73,
     a            subq74,subq75,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*kblsz 
      return
      end 
      subroutine mbicx3 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbicx3 drives the block factorization (version 2, modified)
c     method.
c
      external accel, suba4, suba5, subq70, subq71, subq72, 
     a         subq73, subq74, subq75, noadp
      external ibfcn4
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      call blkdef (coef,jcoef,wksp,iwksp,n,ier)
      if (ier .lt. 0) return
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmz (4,ibfcn4,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*kblsz 
      call split (accel,suba4,suba5,subq70,subq71,subq72,subq73,
     a            subq74,subq75,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*kblsz 
      return
      end 
      subroutine llsp3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... llsp3 drives the line least squares polynomial preconditioner.
c
      external accel, suba4, suba5, subq66, subq67, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      call needw ('llsp3',0,irpnt,n,ier)
      if (ier .lt. 0) return
      call adinfn (n,ndim,maxnz,jcoef,coef,3,ainf,wksp(irpnt))
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call needw ('llsp3',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba4,suba5,subq66,subq67,subq66,subq67,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine lneu3 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... lneu3 drives the line neumann polynomial preconditioner.
c
      external accel, suba4, suba5, subq68, subq69, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call lfact (coef,jcoef,wksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call needw ('lneu3',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*n
      call split (accel,suba4,suba5,subq68,subq69,subq68,subq69,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      return
      end 
      subroutine rich4 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... rich4 drives the richardson preconditioner. 
c
      external accel, suba12, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba12,suba12,copy,copy,copy,copy,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine jac4 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... jac4 drives the jacobi preconditioner.
c
      external accel, suba12, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba12,suba12,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine lsp4 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... lsp4 drives the least squares polynomial preconditioner.
c
      external accel, suba12, sub110, copy, noadp 
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('lsp4',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      call ainfn (n,ndim,maxnz,jcoef,coef,4,ainf,wksp(irpnt))
      iwkpt2 = irpnt
      irpnt = irpnt + 2*n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba12,suba12,sub110,sub110,sub110,sub110,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine neu4 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... neu4 drives the neumann polynomial preconditioner.
c
      external accel, suba12, sub111, copy, noadp 
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('neu4',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      iwkpt2 = irpnt
      irpnt = irpnt + n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba12,suba12,sub111,sub111,sub111,sub111,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine rich5 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... rich5 drives the richardson preconditioner. 
c
      external accel, suba13, suba14, copy, noadp 
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba13,suba14,copy,copy,copy,copy,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine jac5 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... jac5 drives the jacobi preconditioner.
c
      external accel, suba13, suba14, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom4 / keygs, srelpr, keyzer
c
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba13,suba14,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine lsp5 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... lsp5 drives the least squares polynomial preconditioner.
c
      external accel, suba13, suba14, sub112, sub113, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / itcom8 / ainf
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('lsp5',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      call ainfn (n,ndim,maxnz,jcoef,coef,5,ainf,wksp(irpnt))
      iwkpt2 = irpnt
      irpnt = irpnt + 2*n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba13,suba14,sub112,sub113,sub112,sub113,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*n
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine neu5 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... neu5 drives the neumann polynomial preconditioner.
c
      external accel, suba13, suba14, sub114, sub115, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom4 / keygs, srelpr, keyzer
c
      n = nn
      call needw ('neu5',0,irpnt,n,ier) 
      if (ier .lt. 0) return
      iwkpt2 = irpnt
      irpnt = irpnt + n
      iwkpt1 = irpnt
      if (keygs .eq. 1) irpnt = irpnt + 2*n
      call split (accel,suba13,suba14,sub114,sub115,sub114,sub115,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      if (keygs .eq. 1) irpnt = irpnt - 2*n
      return
      end 
      subroutine sor6 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... sor6 drives the multi-color sor method.
c
      external accel, suba8, subq96, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba8,suba8,subq96,subq96,subq96,subq96,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine ssor6 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ssor6 drives the multi-color ssor method.
c
      external accel, suba8, suba9, subq97, subq98, subq99, sub100,
     a         sub101, sub102, sub103
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
c
      iwkpt1 = irpnt
      irpnt = irpnt + n + ncmax
      call split (accel,suba8,suba9,subq97,subq98,subq99,sub100,
     a            sub101,sub102,sub103, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n - ncmax
      return
      end 
      subroutine ic6 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... ic6 drives the ic preconditioner. 
c     (multi-color ordering)
c
      external accel, suba8, suba9, sub104, sub105, sub106, 
     a         sub107, sub108, sub109, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call pfactc (coef,jcoef,wksp,iwksp,n,1,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba8,suba9,sub104,sub105,sub106,sub107,
     a            sub108,sub109,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine mic6 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... mic6 drives the mic preconditioner.
c     (multi-color ordering)
c
      external accel, suba8, suba9, sub104, sub105, sub106, 
     a         sub107, sub108, sub109, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
c
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call pfactc (coef,jcoef,wksp,iwksp,n,2,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba8,suba9,sub104,sub105,sub106,sub107,
     a            sub108,sub109,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine rs6 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... rs6 drives the reduced system method (purdue storage
c     with red-black coloring).
c
      external accel, suba10, suba11, subq1, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
c
c ... compute red-black rhs.
c
      nr = iwksp(nc)
      nb = n - nr
      call needw ('rs6',0,irpnt,2*n,ier)
      if (ier .lt. 0) return
      irhs = irpnt
      irpnt = irpnt + nr
      call vfill (2*n,wksp(irhs),0.0)
      call rsbegp (n,nr,ndim,maxnz,jcoef,coef,wksp(irhs),rhs,
     a             wksp(irpnt))
      iwkpt1 = irpnt
      irpnt = irpnt + n + nb
      call split (accel,suba10,suba11,subq1,subq1,subq1,subq1,
     a            copy,copy,noadp,
     a            coef,jcoef,nr,u,ubar,wksp(irhs),wksp,iwksp,
     a            iparm,rparm,ier)
      call rsendp (n,nr,ndim,maxnz,jcoef,coef,u,rhs,wksp(iwkpt1))
      irpnt = irpnt - 2*n
      return
      end 
      subroutine sor7 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                 iparm,rparm,ier) 
c
c ... sor7 drives the multi-color sor method.
c
      external accel, suba2, subq26, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call mfact (coef,jcoef,wksp,iwksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      call split (accel,suba2,suba2,subq26,subq26,subq26,subq26,
     a            copy,copy,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      return
      end 
      subroutine ssor7 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... ssor7 drives the multi-color ssor method.
c
      external accel, suba2, suba3, subq27, subq28, subq29, 
     a         subq30, subq31, subq32, subq33
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call mfact (coef,jcoef,wksp,iwksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + n
      call split (accel,suba2,suba3,subq27,subq28,subq29,subq30,
     a            subq31,subq32,subq33, 
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - n
      return
      end 
      subroutine bic7 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bic7 drives the block factorization (version 1) method.
c     (multi-color ordering)
c
      external accel, suba2, suba3, subq34, subq35, subq36, 
     a         subq37, subq38, subq39, noadp
      external ibfcn1
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmy (1,ibfcn1,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*ncmax 
      call split (accel,suba2,suba3,subq34,subq35,subq36,subq37,
     a            subq38,subq39,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*ncmax 
      return
      end 
      subroutine mbic7 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbic7 drives the block factorization (version 1, modified)
c     method (multi-color ordering)
c
      external accel, suba2, suba3, subq34, subq35, subq36, 
     a         subq37, subq38, subq39, noadp
      external ibfcn3
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmy (2,ibfcn3,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*ncmax 
      call split (accel,suba2,suba3,subq34,subq35,subq36,subq37,
     a            subq38,subq39,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*ncmax 
      return
      end 
      subroutine bicx7 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                  iparm,rparm,ier)
c
c ... bicx7 drives the block factorization (version 2)
c     method (multi-color ordering)
c
      external accel, suba2, suba3, subq34, subq35, subq36, 
     a         subq37, subq38, subq39, noadp
      external ibfcn2
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmy (3,ibfcn2,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*ncmax 
      call split (accel,suba2,suba3,subq34,subq35,subq36,subq37,
     a            subq38,subq39,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*ncmax 
      return
      end 
      subroutine mbicx7 (accel,coef,jcoef,n,u,ubar,rhs,wksp,iwksp,
     a                   iparm,rparm,ier)
c
c ... mbicx7 drives the block factorization (version 2, modified)
c     method (multi-color ordering)
c
      external accel, suba2, suba3, subq34, subq35, subq36, 
     a         subq37, subq38, subq39, noadp
      external ibfcn4
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
c
      t1 = timer (dummy)
      if (ifact .eq. 1) call bfacmy (4,ibfcn4,coef,jcoef,wksp,iwksp,
     a                              n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
      iwkpt1 = irpnt
      irpnt = irpnt + 2*ncmax 
      call split (accel,suba2,suba3,subq34,subq35,subq36,subq37,
     a            subq38,subq39,noadp,
     a            coef,jcoef,n,u,ubar,rhs,wksp,iwksp,iparm,rparm,ier) 
      irpnt = irpnt - 2*ncmax 
      return
      end 
      subroutine rs7 (accel,coef,jcoef,nn,u,ubar,rhs,wksp,iwksp,
     a                iparm,rparm,ier)
c
c ... rs7 drives the reduced system method (diagonal storage
c     with red-black coloring).
c
      external accel, suba6, suba7, subq76, subq77, copy, noadp
      integer   iparm(30), jcoef(2), iwksp(1)
      dimension rhs(1), u(1), ubar(1), rparm(30), coef(1), wksp(1)
c
      common / dscons / ndim, mdim, maxnz
      common / cwkcon / lenr, irpnt, irmax, leni, iipnt, iimax
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
c
      n = nn
      t1 = timer (dummy)
      if (ifact .eq. 1) call mfact (coef,jcoef,wksp,iwksp,n,ier)
      t2 = timer (dummy)
      timfac = t2 - t1
      if (ier .lt. 0) return
c
c ... compute red-black rhs.
c
      nr = iwksp(nc)
      nb = n - nr
      call needw ('rs7',0,irpnt,n,ier)
      if (ier .lt. 0) return
      irhs = irpnt
      irpnt = irpnt + nr
      call rsbegd (n,n,nr,ndim,iwksp(maxnew),ndt,ndb,iwksp(jcnew),
     a             coef,wksp(irhs),rhs,wksp(ifactr),wksp(irpnt))
      iwkpt1 = irpnt
      irpnt = irpnt + nb
      call split (accel,suba6,suba7,subq76,subq77,subq76,subq77,
     a            copy,copy,noadp,
     a            coef,jcoef,nr,u,ubar,wksp(irhs),wksp,iwksp,
     a            iparm,rparm,ier)
      call rsendd (n,n,nr,ndim,iwksp(maxnew),ndt,ndb,iwksp(jcnew),
     a             coef,u,rhs,wksp(ifactr))
      irpnt = irpnt - n
      return
      end 
      subroutine suba1 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba1 calls mult2s.
c
      common / dscons / ndim, mdim, maxnz
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mult2s (ndim,maxnz,coef,jcoef,n,x,y)
      return
      end 
      subroutine suba2 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba2 calls muldc.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call muldc (ndim,n,coef,ncolor,iwksp(nc),iwksp(maxnew),
     a            iwksp(jcnew),x,y)
      return
      end 
      subroutine suba3 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba3 calls muldct.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call muldct (ndim,n,coef,ncolor,iwksp(nc),iwksp(maxnew),
     a             iwksp(jcnew),x,y)
      return
      end 
      subroutine suba4 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba4 calls mult2n.
c
      common / dscons / ndim, mdim, maxnz
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mult2n (ndim,maxnz,coef,jcoef,n,x,y)
      return
      end 
      subroutine suba5 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba5 calls mul2nt.
c
      common / dscons / ndim, mdim, maxnz
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mul2nt (ndim,maxnz,coef,jcoef,n,x,y)
      return
      end 
      subroutine suba6 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba6 calls rsad.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      nr = iwksp(nc)
      nb = iwksp(nc+1)
      nbig = nr + nb
      call rsad (nbig,n,n,ndim,iwksp(maxnew),ndt,ndb,
     a           iwksp(jcnew),coef,y,x,wksp(ifactr),wksp(iwkpt1))
      return
      end 
      subroutine suba7 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba7 calls rsatd.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      nr = iwksp(nc)
      nb = iwksp(nc+1)
      nbig = nr + nb
      call rsatd (nbig,n,n,ndim,iwksp(maxnew),ndt,ndb,
     a           iwksp(jcnew),coef,y,x,wksp(ifactr),wksp(iwkpt1))
      return
      end 
      subroutine suba8 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba8 calls mult1.
c
      common / dscons / ndim, mdim, maxnz
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mult1 (ndim,maxnz,coef,jcoef,wksp(iwkpt1),n,x,y) 
      return
      end 
      subroutine suba9 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba9 calls mul1t.
c
      common / dscons / ndim, mdim, maxnz
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mul1t (ndim,maxnz,coef,jcoef,wksp(iwkpt1),n,x,y) 
      return
      end 
      subroutine suba10 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba10 calls rsap.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      nr = iwksp(nc)
      nb = iwksp(nc+1)
      nbig = nr + nb
      call rsap (ndim,nbig,n,maxnz,jcoef,coef,x,y,wksp(iwkpt1))
      return
      end 
      subroutine suba11 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba11 calls rsatp.
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      nr = iwksp(nc)
      nb = iwksp(nc+1)
      nbig = nr + nb
      if (isymm .eq. 0) call rsap (ndim,nbig,n,maxnz,jcoef,coef,
     a                                   x,y,wksp(iwkpt1))
      if (isymm .eq. 1) call rsatp (ndim,nbig,n,maxnz,jcoef,coef,
     a                                   x,y,wksp(iwkpt1))
      return
      end 
      subroutine suba12 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba12 calls mult3.
c
      common / dscons / ndim, mdim, maxnz
      common / cmpart / mpstrt, mpart
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mult3 (mpart,iwksp(mpstrt),coef,jcoef,jcoef(ndim+1),
     a            wksp(iwkpt1),x,y)
      return
      end 
      subroutine suba13 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba13 calls mult3n.
c
      common / dscons / ndim, mdim, maxnz
      common / cmpart / mpstrt, mpart
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mult3n (mpart,iwksp(mpstrt),coef,jcoef,jcoef(ndim+1),
     a             wksp(iwkpt1),x,y)
      return
      end 
      subroutine suba14 (coef,jcoef,wksp,iwksp,n,x,y)
c
c ... suba14 calls mul3nt.
c
      common / dscons / ndim, mdim, maxnz
      common / cmpart / mpstrt, mpart
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension x(1), y(1), coef(1), wksp(1)
c
      call mul3nt (mpart,iwksp(mpstrt),coef,jcoef,jcoef(ndim+1),
     a             wksp(iwkpt1),x,y)
      return
      end 
      subroutine subq1 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq1 calls pjac, for jacobi preconditioning.
c
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      call pjac (coef,n,r,z)
      return
      end 
      subroutine subq2 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq2 calls bdsol, for line jacobi preconditioning.
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (nstore .eq. 2) isym = 0
      if (nstore .eq. 3) isym = 1
      call bdsol (n,n,kblsz,ndt,ndb,wksp(ifactr),r,z,isym)
      return
      end 
      subroutine subq3 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq3 calls bdsolt, for line jacobi preconditioning.
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      call bdsolt (n,n,kblsz,ndt,ndb,wksp(ifactr),r,z)
      return
      end 
      subroutine subq4 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq4 call bmul or bmuln, for line jacobi preconditioning
c     (approximate inverse)
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (nstore .eq. 2) isym = 0
      if (nstore .eq. 3) isym = 1
      ift = ifactr + n
      ifb = ifactr + n*(ndt + 1)
      if (isym .eq. 0) call bmul (n,n,ndt,wksp(ifactr),wksp(ift),r,z) 
      if (isym .eq. 1) call bmuln (n,n,ndt,ndb,wksp(ifactr),
     a                             wksp(ift),wksp(ifb),r,z) 
      return
      end 
      subroutine subq5 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq5 call bmul or bmulnt, for line jacobi preconditioning
c     (approximate inverse)
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (nstore .eq. 2) isym = 0
      if (nstore .eq. 3) isym = 1
      ift = ifactr + n
      ifb = ifactr + n*(ndt + 1)
      if (isym .eq. 0) call bmul (n,n,ndt,wksp(ifactr),wksp(ift),r,z) 
      if (isym .eq. 1) call bmulnt (n,n,ndt,ndb,wksp(ifactr),
     a                             wksp(ift),wksp(ifb),r,z) 
      return
      end 
      subroutine subq6 (coef,jcoef,wksp,iwksp,n,u,rhs,unew) 
c
c ... subq6 calls the basic sor iterative step
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension u(1), rhs(1), unew(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call sords (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a            irwise,u,rhs,unew,iwksp(iwkpt1))
      return
      end 
      subroutine subq7 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq7 calls the ssor preconditioner.
c
c
c *** begin -- itpack common
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c *** end   -- itpack common
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call srs (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a          irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq8 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq8 calls the ssor preconditioner.
c
c
c *** begin -- itpack common
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c *** end   -- itpack common
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call srs1 (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a          irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq9 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq9 calls the ssor preconditioner.
c
c
c *** begin -- itpack common
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c *** end   -- itpack common
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call srs3 (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a          irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq10 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq10 calls the ssor preconditioner.
c
c
c *** begin -- itpack common
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c *** end   -- itpack common
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call srs2 (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a          irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq11 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq11 calls the ssor preconditioner.
c
c
c *** begin -- itpack common
c
      logical           omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c *** end   -- itpack common
c
      common / dscons / ndim, mdim, maxnz
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call srs4 (ndim,n,maxt,jcoef(2),coef,coef(ndim+1),omega,
     a          irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq12 (coef,jcoef,wksp,iwksp,n,p,r,pdp,pldup)
c
c ... subq12 calls the ssor adaption routine.
c
c
      common / dscons / ndim, mdim, maxnz
      integer jcoef(2), iwksp(1)
      dimension p(1), r(1), coef(1), wksp(1)
c
      maxt = maxnz - 1
      call ssord (ndim,maxt,jcoef(2),coef,coef(ndim+1),n,p,r,
     a            pdp,pldup)
      return
      end 
      subroutine subq13 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq13 calls ics, for ic(s) preconditioning.
c
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (propa) call ics (ndim,n,maxt,jcoef(2),wksp(ifactr),
     a                     coef(ndim+1),1,irwise,iwksp(iwkpt1),r,z)
      if (.not. propa) call ics (n,n,maxt,iwksp(ifacti+1),
     a                           wksp(ifactr),wksp(ifactr+n),
     a                           0,irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq14 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq14 calls ics1, for ic(s) preconditioning.
c
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (propa) call ics1 (ndim,n,maxt,jcoef(2),wksp(ifactr),
     a                     coef(ndim+1),1,irwise,iwksp(iwkpt1),r,z)
      if (.not. propa) call ics1 (n,n,maxt,iwksp(ifacti+1), 
     a                           wksp(ifactr),wksp(ifactr+n),
     a                           0,irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq15 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq15 calls ics3, for ic(s) preconditioning.
c
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (propa) call ics3 (ndim,n,maxt,jcoef(2),wksp(ifactr),
     a                     coef(ndim+1),1,irwise,iwksp(iwkpt1),r,z)
      if (.not. propa) call ics3 (n,n,maxt,iwksp(ifacti+1), 
     a                           wksp(ifactr),wksp(ifactr+n),
     a                           0,irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq16 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq16 calls ics2, for ic(s) preconditioning.
c
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (propa) call ics2 (ndim,n,maxt,jcoef(2),wksp(ifactr),
     a                     coef(ndim+1),1,irwise,iwksp(iwkpt1),r,z)
      if (.not. propa) call ics2 (n,n,maxt,iwksp(ifacti+1), 
     a                           wksp(ifactr),wksp(ifactr+n),
     a                           0,irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq17 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq17 calls ics4, for ic(s) preconditioning.
c
c
      common / dscons / ndim, mdim, maxnz
      logical           propa 
      common / cblock / propa, ncolor, maxd, nc, ipt, maxnew,
     a                  jcnew, lbhb, iblock, ncmax
      common / intern / ndt, ndb, maxt, maxb, ivers, irwise 
      common / cfactr / nfactr, nfacti, ifactr, ifacti, timfac
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1), wksp(1)
c
      if (propa) call ics4 (ndim,n,maxt,jcoef(2),wksp(ifactr),
     a                     coef(ndim+1),1,irwise,iwksp(iwkpt1),r,z)
      if (.not. propa) call ics4 (n,n,maxt,iwksp(ifacti+1), 
     a                           wksp(ifactr),wksp(ifactr+n),
     a                           0,irwise,iwksp(iwkpt1),r,z)
      return
      end 
      subroutine subq18 (coef,jcoef,wksp,iwksp,n,r,z)
c
c ... subq18 calls ppii, for lspoly preconditioning.
c
c
      common / itcom6 / method, iscale, iperm, nstore,
     a                  ifact, kblsz, lvfill, ltrunc, ndeg, 
     a                  ipropa, isymm, ifctv
      common / point  / iptscl, iwkpt1, iwkpt2, iwkpt3
      common / itcom8 / ainf
      integer jcoef(2), iwksp(1)
      dimension r(1), z(1), coef(1