subroutine jcg (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp,
     *                iparm,rparm,ier)
c
c     itpackv 2d main routine  jcg  (jacobi conjugate gradient)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          jcg drives the jacobi conjugate gradient algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  jacobi conjugate
c                 gradient needs this to be in length at least
c                 4*n + 4*itmax.  here itmax = iparm(1) is the
c                 maximum allowable number of iterations.
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... jcg module references --
c
c         from itpackv    chgcon, determ, dfault, echall,
c                         eigvns, eigvss, eqrt1s, iterm ,
c                         itjcg , parcon, permat, perror,
c                         pervec, pjac  , pmult , prbndx, pstop ,
c                         sbelm , scal  , sum3  , unscal,
c                         vout  , zbrent
c          system         abs, alog10, amax0, amax1, mod, sqrt
c
c ... local itpackv references --
c
c          echall, itjcg , permat,
c          perror, pervec, pjac  , prbndx, sbelm , scal  ,
c          unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      j c g      ')
      if (iparm(1).le.0) go to 370
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *      '    in itpackv routine jcg'/1x,
     *      '    rparm(1) =',e10.3,' (zeta)'/1x,
     *      '    a value this small may hinder convergence '/1x,
     *      '    since machine precision srelpr =',e10.3/1x,
     *      '    zeta reset to ',e10.3) 
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 11
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *        '    called from itpackv routine jcg '/1x,
     *        '    invalid matrix dimension, n =',i8)
      go to 370
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jcg '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 370
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
      ib4 = ib3+n
      ib5 = ib4+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jcg  '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 350
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 14
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jcg      '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 350
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... check for sufficient workspace.
c
  170 iparm(8) = 4*n+4*itmax
      if (nw.ge.iparm(8)) go to 190
      ier = 12
      if (level.ge.0) write (nout,180) nw,iparm(8)
  180 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jcg '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 330
c
  190 continue
      if (level.le.2) go to 220
      write (nout,200)
  200 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
      if (adapt) write (nout,210)
  210 format (1x,'cme is the estimate of the largest eigenvalue of',
     *   ' the jacobi matrix')
  220 if (iparm(11).eq.0) timi1 = timer(0.0)
c
c ... compute initial pseudo-residual
c
      do 230 i = 1,nw
         wksp(i) = 0.0
  230 continue
      call scopy (n,rhs,1,wksp(ib2),1)
      call pjac (n,ndim,maxnz,jcoef,coef,u,wksp(ib2))
      do 240 i = 1,n
         wksp(n+i) = wksp(n+i)-u(i)
  240 continue
c
c ... iteration sequence
c
      itmax1 = itmax+1
      do 260 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 250
c
c ... code for the even iterations.
c
c     u           = u(in)             wksp(ib2) = del(in)
c     wksp(ib1)   = u(in-1)           wksp(ib3) = del(in-1) 
c
         call itjcg (n,ndim,maxnz,jcoef,coef,u,wksp(ib1),wksp(ib2),
     *      wksp(ib3),wksp(ib4),wksp(ib5))
c
         if (halt) go to 290
         go to 260
c
c ... code for the odd iterations.
c
c     u           = u(in-1)           wksp(ib2) = del(in-1) 
c     wksp(ib1)   = u(in)             wksp(ib3) = del(in)
c
  250    call itjcg (n,ndim,maxnz,jcoef,coef,wksp(ib1),u,wksp(ib3),
     *      wksp(ib2),wksp(ib4),wksp(ib5))
c
         if (halt) go to 290
  260 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 270
      timi2 = timer(0.0)
      time1 = timi2-timi1
  270 ier = 13
      if (level.ge.1) write (nout,280) itmax
  280 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine jcg'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 320
c
c ... method has converged
c
  290 if (iparm(11).ne.0) go to 300
      timi2 = timer(0.0)
      time1 = timi2-timi1
  300 if (level.ge.1) write (nout,310) in
  310 format (/1x,'jcg  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  320 continue
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
c
c ... un-permute matrix,rhs, and solution
c
  330 if (iparm(9).ne.-1) go to 340
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp(ib4),
     *   iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib4))
      call pervec (n,iwksp(ib2),u,wksp(ib4))
      if (ier.eq.12) go to 350
c
c ... optional error analysis 
c
  340 idgts = iparm(12)
      if (idgts.lt.0) go to 350
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  350 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      iparm(8) = iparm(8)-4*(itmax-in)
      if (iparm(11).ne.0) go to 360
      timj2 = timer(0.0)
      time2 = timj2-timj1
  360 if (iparm(3).ne.0) go to 370
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  370 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,380)
  380 format (/1x,'execution successful')
c
      return
      end 
      subroutine jsi (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp,
     *                iparm,rparm,ier)
c
c     itpackv 2d main routine  jsi  (jacobi semi-iterative)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          jsi drives the jacobi semi-iteration algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  jacobi si 
c                 needs this to be in length at least 2*n.
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... jsi module references --
c
c         from itpackv    cheby , chgsi , chgsme, dfault, echall,
c                         iterm , itjsi ,
c                         parsi , permat, perror, pervec, pjac  ,
c                         pmult , prbndx, pstop , pvtbv ,
c                         sbelm , scal  , sum3  , 
c                         tstchg, unscal,
c                         vout
c          system         abs, alog10, amax0, amax1, float, mod, sqrt 
c
c ... local itpackv references --
c
c          echall, itjsi , permat,
c          perror, pervec, prbndx, sbelm , scal  ,
c          unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      j s i      ')
      if (iparm(1).le.0) go to 350
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine jsi'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 21
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jsi '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 350
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jsi '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 350
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jsi  '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 330
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 24
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jsi      '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 330
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... check for sufficient workspace.
c
  170 iparm(8) = 2*n
      if (nw.ge.iparm(8)) go to 190
      ier = 22
      if (level.ge.0) write (nout,180) nw,iparm(8)
  180 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine jsi '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 310
c
  190 continue
      if (level.le.2) go to 210
      write (nout,200)
  200 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
  210 if (iparm(11).eq.0) timi1 = timer(0.0)
      do 220 i = 1,nw
         wksp(i) = 0.0
  220 continue
c
c ... iteration sequence
c
      itmax1 = itmax+1
      do 240 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 230
c
c ... code for the even iterations.
c
c     u           = u(in)
c     wksp(ib1)   = u(in-1)
c
         call itjsi (n,ndim,maxnz,jcoef,coef,rhs,u,wksp(ib1),wksp(ib2),
     *      icnt)
c
         if (halt) go to 270
         go to 240
c
c ... code for the odd iterations.
c
c     u           = u(in-1)
c     wksp(ib1)   = u(in)
c
  230    call itjsi (n,ndim,maxnz,jcoef,coef,rhs,wksp(ib1),u,wksp(ib2),
     *      icnt)
c
         if (halt) go to 270
  240 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 250
      timi2 = timer(0.0)
      time1 = timi2-timi1
  250 ier = 23
      if (level.ge.1) write (nout,260) itmax
  260 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine jsi'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 300
c
c ... method has converged
c
  270 if (iparm(11).ne.0) go to 280
      timi2 = timer(0.0)
      time1 = timi2-timi1
  280 if (level.ge.1) write (nout,290) in
  290 format (/1x,'jsi  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  300 continue
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
c
c ... un-permute matrix,rhs, and solution
c
  310 if (iparm(9).ne.-1) go to 320
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp(ib2),
     *   iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.22) go to 330
c
c ... optional error analysis 
c
  320 idgts = iparm(12)
      if (idgts.lt.0) go to 330
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  330 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      if (iparm(11).ne.0) go to 340
      timj2 = timer(0.0)
      time2 = timj2-timj1
  340 if (iparm(3).ne.0) go to 350
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  350 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,360)
  360 format (/1x,'execution successful')
c
      return
      end 
      subroutine sor (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp,
     *                iparm,rparm,ier)
c
c     itpackv 2d main routine  sor  (successive overrelaxation)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          sor drives the successive overrelaxation algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  sor
c                 needs this to be in length at least 2*n.
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... sor module references --
c
c         from itpackv    dfault, echall, ipstr , iterm ,
c                         itsor , move  , permat, 
c                         perror, pervec, pfsor1, pmult , prbndx,
c                         pstop , sbelm , scal  , 
c                         tau   , unscal,
c                         vout
c          system         abs, alog10, amax0, amax1, float, iabs, mod,
c                         sqrt
c
c ... local itpackv references --
c
c          echall, itsor , move  ,
c          permat, perror, pervec, prbndx, sbelm , scal  ,
c          unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c ... specifications for local variables
c
      logical moved 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      s o r      ')
      if (iparm(1).le.0) go to 360
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine sor'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 31
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine sor '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 360
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine sor '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 360
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine sor  '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 340
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 34
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine sor      '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 340
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... re-arrange data structure if natural ordering is used.
c
  170 if (nb.ge.0) go to 180
      call move (n,ndim,maxnz,nlow,nup,jcoef,coef,level,nout,moved,
     *   iwksp)
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,1)
c
c ... check for sufficient workspace.
c
  180 iparm(8) = 2*n
      if (nw.ge.iparm(8)) go to 200
      ier = 32
      if (level.ge.0) write (nout,190) nw,iparm(8)
  190 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine sor '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 310
c
  200 continue
      if (level.le.2) go to 230
      if (adapt) write (nout,210)
  210 format (///1x,'cme is the estimate of the largest eigenvalue of',
     *   ' the jacobi matrix')
      write (nout,220)
  220 format (1x,'omega is the relaxation factor')
  230 if (iparm(11).eq.0) timi1 = timer(0.0)
      do 240 i = 1,nw
         wksp(i) = 0.0
  240 continue
c
c ... iteration sequence
c
      itmax1 = itmax+1
      do 250 loop = 1,itmax1
         in = loop-1
c
c ... code for one iteration. 
c
c         u          = u(in)
c
         call itsor (n,nb,ndim,maxnz,jcoef,coef,rhs,u,nup,moved,
     *      wksp,iwksp(ib3))
c
         if (halt) go to 280
  250 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 260
      timi2 = timer(0.0)
      time1 = timi2-timi1
  260 ier = 33
      if (level.ge.1) write (nout,270) itmax
  270 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine sor'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 310
c
c ... method has converged
c
  280 if (iparm(11).ne.0) go to 290
      timi2 = timer(0.0)
      time1 = timi2-timi1
  290 if (level.ge.1) write (nout,300) in
  300 format (/1x,'sor  has converged in ',i5,' iterations')
c
c ... un-permute matrix,rhs, and solution
c
  310 if (iparm(9).ne.-1) go to 320
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp,iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.32) go to 340
c
c ... undo wave front reordering.
c
  320 if (nb.ge.0) go to 330
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,2)
c
c ... optional error analysis 
c
  330 idgts = iparm(12)
      if (idgts.lt.0) go to 340
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  340 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      if (iparm(11).ne.0) go to 350
      timj2 = timer(0.0)
      time2 = timj2-timj1
  350 if (iparm(3).ne.0) go to 360
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(5) = omega
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  360 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,370)
  370 format (/1x,'execution successful')
c
      return
      end 
      subroutine ssorcg (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp, 
     *                   iparm,rparm,ier)
c
c     itpackv 2d main routine  ssorcg  (symmetric successive over-
c                                     relaxation conjugate gradient)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          ssorcg drives the symmetric sor-cg algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  ssorcg needs
c                 this to be in length at least 6*n + 4*itmax.
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... ssorcg module references --
c
c         from itpackv    chgcon, determ, dfault, echall,
c                         eigvns, eigvss, eqrt1s, iterm ,
c                         itsrcg, move  , omeg  , omgchg,
c                         omgstr, parcon, pbeta , pbsor , permat,
c                         perror, pervec, pfsor , pjac  , pmult ,
c                         prbndx, pstop , pvtbv , sbelm , scal  ,
c                         sum3  , unscal,
c                           vout  ,
c                         zbrent
c          system         abs, alog, alog10, amax0, amax1, amin1,
c                         mod, sqrt
c
c ... local itpackv references --
c
c          echall, itsrcg, move  ,
c          omeg  , pbeta , permat, perror, pervec, pfsor ,
c          prbndx, sbelm , scal  , unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c ... specifications for local variables
c
      logical moved 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (iparm(9).ge.-1) iparm(6) = 2
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      s s o r c g')
      if (iparm(1).le.0) go to 420
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine ssorcg'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 41
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorcg '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 420
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorcg '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 420
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
      ib4 = ib3+n
      ib5 = ib4+n
      ib6 = ib5+n
      ib7 = ib6+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorcg  '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 400
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 44
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorcg   '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 400
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... re-arrange data structure if natural ordering is used.
c
  170 if (nb.ge.0) go to 180
      call move (n,ndim,maxnz,nlow,nup,jcoef,coef,level,nout,moved,
     *   iwksp)
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,1)
c
c ... check for sufficient workspace.
c
  180 iparm(8) = 6*n+4*itmax
      if (nw.ge.iparm(8)) go to 200
      ier = 42
      if (level.ge.0) write (nout,190) nw,iparm(8)
  190 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorcg '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 370
c
  200 continue
      if (level.le.2) go to 230
      write (nout,210)
  210 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
      write (nout,220)
  220 format (1x,'s-prime is an initial estimate for new cme')
  230 if (iparm(11).eq.0) timi1 = timer(0.0)
c
c ... special procedure for fully adaptive case.
c
      do 240 i = 1,nw
         wksp(i) = 0.0
  240 continue
      if (.not.adapt) go to 270
      if (.not.betadt) go to 260
      do 250 i = 1,n
         wksp(i) = 1.0
  250 continue
      betnew = pbeta(n,nb,ndim,maxnz,jcoef,coef,wksp(ib1),nup,nlow,moved
     *   ,wksp(ib2),wksp(ib3))/float(n) 
      betab = amax1(betab,0.25,betnew)
  260 call omeg (0.0,1)
      is = 0
c
c ... initialize forward pseudo-residual
c
  270 call sorscl (n,ndim,maxnz,coef,rhs,omega)
      call scopy (n,rhs,1,wksp,1)
      call scopy (n,u,1,wksp(ib2),1)
      call pfsor (n,nb,ndim,maxnz,jcoef,coef,wksp(ib2),wksp(ib1),nup, 
     *   moved,iwksp(ib3))
      do 280 i = 1,n
         wksp(n+i) = wksp(n+i)-u(i)
  280 continue
c
c ... iteration sequence
c
      itmax1 = itmax+1
      do 300 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 290
c
c ... code for the even iterations.
c
c     u           = u(in)       wksp(ib2) = c(in) 
c     wksp(ib1)   = u(in-1)     wksp(ib3) = c(in-1)
c
         call itsrcg (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,u,
     *      wksp(ib1),wksp(ib2),wksp(ib3),wksp(ib4),wksp(ib5),moved,
     *      wksp(ib6),wksp(ib7),iwksp(ib3))
c
         if (halt) go to 330
         go to 300
c
c ... code for the odd iterations.
c
c     u           = u(in-1)     wksp(ib2) = c(in-1)
c     wksp(ib1)   = u(in)       wksp(ib3) = c(in) 
c
  290    call itsrcg (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,wksp(ib1),
     *      u,wksp(ib3),wksp(ib2),wksp(ib4),wksp(ib5),moved,wksp(ib6),
     *      wksp(ib7),iwksp(ib3))
c
         if (halt) go to 330
  300 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 310
      timi2 = timer(0.0)
      time1 = timi2-timi1
  310 ier = 43
      if (level.ge.1) write (nout,320) itmax
  320 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine ssorcg'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 360
c
c ... method has converged
c
  330 if (iparm(11).ne.0) go to 340
      timi2 = timer(0.0)
      time1 = timi2-timi1
  340 if (level.ge.1) write (nout,350) in
  350 format (/1x,'ssorcg  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  360 factor = 1.0/omega
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
c
c ... un-permute matrix,rhs, and solution
c
  370 if (iparm(9).ne.-1) go to 380
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp,iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.42) go to 400
c
c ... undo wave front reordering.
c
  380 if (nb.ge.0) go to 390
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,2)
c
c ... optional error analysis 
c
  390 idgts = iparm(12)
      if (idgts.lt.0) go to 400
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  400 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      if (iparm(11).ne.0) go to 410
      timj2 = timer(0.0)
      time2 = timj2-timj1
  410 iparm(8) = iparm(8)-4*(itmax-in)
      if (iparm(3).ne.0) go to 420
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(5) = omega
      rparm(6) = specr
      rparm(7) = betab
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  420 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,430)
  430 format (/1x,'execution successful')
c
      return
      end 
      subroutine ssorsi (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp, 
     *                   iparm,rparm,ier)
c
c     itpackv 2d main routine  ssorsi  (symmetric successive over-
c                                     relaxation semi-iteration)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          ssorsi drives the symmetric sor-si algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  ssorsi
c                 needs this to be in length at least 5*n
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... ssorsi module references --
c
c         from itpackv    cheby , chgsi , dfault, echall,
c                         iterm , itsrsi, move  , 
c                         omeg  , omgstr, parsi , pbeta , pbsor ,
c                         permat, perror, pervec, pfsor , pjac  ,
c                         pmult , prbndx, pstop , pvtbv , sbelm ,
c                         scal  , sum3  , tstchg, 
c                         unscal,
c                         vout
c          system         abs, alog, alog10, amax0, amax1, float,
c                         mod, sqrt
c
c ... local itpackv references --
c
c          echall, itsrsi, move  ,
c          omeg  , pbeta , permat, perror, pervec, prbndx,
c          sbelm , scal  , unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c ... specifications for local variables
c
      logical moved 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      if (iparm(9).ge.-1) iparm(6) = 2
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      s s o r s i')
      if (iparm(1).le.0) go to 400
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine ssorsi'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 51
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorsi '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 400
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorsi '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 400
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
      ib4 = ib3+n
      ib5 = ib4+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorsi  '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 380
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 54
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorsi   '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 380
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... re-arrange data structure if natural ordering is used.
c
  170 if (nb.ge.0) go to 180
      call move (n,ndim,maxnz,nlow,nup,jcoef,coef,level,nout,moved,
     *   iwksp)
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,1)
c
c ... check for sufficient workspace.
c
  180 iparm(8) = 5*n
      if (nw.ge.iparm(8)) go to 200
      ier = 52
      if (level.ge.0) write (nout,190) nw,iparm(8)
  190 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine ssorsi '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 350
c
  200 continue
      if (level.le.2) go to 220
      write (nout,210)
  210 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
  220 if (iparm(11).eq.0) timi1 = timer(0.0)
c
c ... special procedure for fully adaptive case.
c
      do 230 i = 1,nw
         wksp(i) = 0.0
  230 continue
      if (.not.adapt) go to 260
      if (.not.betadt) go to 250
      do 240 i = 1,n
         wksp(i) = 1.0
  240 continue
      betnew = pbeta(n,nb,ndim,maxnz,jcoef,coef,wksp(ib1),nup,nlow,moved
     *   ,wksp(ib2),wksp(ib3))/float(n) 
      betab = amax1(betab,0.25,betnew)
  250 call omeg (0.0,1)
      is = 0
c
c ... iteration sequence
c
  260 call sorscl (n,ndim,maxnz,coef,rhs,omega)
      itmax1 = itmax+1
      do 280 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 270
c
c ... code for the even iterations.
c
c     u           = u(in)
c     wksp(ib1)   = u(in-1)
c
         call itsrsi (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,u,
     *      wksp(ib1),wksp(ib2),wksp(ib3),wksp(ib4),moved,wksp(ib5),
     *      iwksp(ib3))
c
         if (halt) go to 310
         go to 280
c
c ... code for the odd iterations.
c
c     u           = u(in-1)
c     wksp(ib1)   = u(in)
c
  270    call itsrsi (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,wksp(ib1),
     *      u,wksp(ib2),wksp(ib3),wksp(ib4),moved,wksp(ib5),iwksp(ib3))
c
         if (halt) go to 310
  280 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 290
      timi2 = timer(0.0)
      time1 = timi2-timi1
  290 ier = 53
      if (level.ge.1) write (nout,300) itmax
  300 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine ssorsi'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 340
c
c ... method has converged
c
  310 if (iparm(11).ne.0) go to 320
      timi2 = timer(0.0)
      time1 = timi2-timi1
  320 if (level.ge.1) write (nout,330) in
  330 format (/1x,'ssorsi  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  340 factor = 1.0/omega
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
c
c ... un-permute matrix,rhs, and solution
c
  350 if (iparm(9).ne.-1) go to 360
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp,iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.52) go to 380
c
c ... undo wave front reordering.
c
  360 if (nb.ge.0) go to 370
      call sorwav (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,iwksp(ib2),
     *   iwksp(ib3),wksp,numwav,2)
c
c ... optional error analysis 
c
  370 idgts = iparm(12)
      if (idgts.lt.0) go to 380
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  380 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      if (iparm(11).ne.0) go to 390
      timj2 = timer(0.0)
      time2 = timj2-timj1
  390 if (iparm(3).ne.0) go to 400
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(5) = omega
      rparm(6) = specr
      rparm(7) = betab
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  400 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,410)
  410 format (/1x,'execution successful')
c
      return
      end 
      subroutine rscg (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp,
     *                 iparm,rparm,ier) 
c
c     itpackv 2d main routine  rscg  (reduced system conjugate
c                                    gradient)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          rscg drives the reduced system cg algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  rscg needs this
c                 to be in length at least n+3*nb+4*itmax.  here
c                 itmax = iparm(1) and nb is the order of the black
c                 subsystem.  itmax is the maximum allowable number of
c                 iterations. 
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... rscg module references --
c
c         from itpackv    chgcon, determ, dfault, echall,
c                         eigvns, eigvss, eqrt1s, iterm ,
c                         itrscg, parcon, permat, perror,
c                         pervec, pmult , prbndx, prsblk, prsred,
c                         pstop , sbelm , scal  , 
c                         sum3  , unscal,
c                         vout  , zbrent
c          system         abs, alog10, amax0, amax1, mod, sqrt
c
c ... local itpackv references --
c
c          echall, itrscg, permat,
c          perror, pervec, prbndx, prsblk, prsred, sbelm ,
c          scal  , unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      r s c g    ')
      if (iparm(1).le.0) go to 380
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine rscg'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 61
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rscg '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 380
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rscg '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 380
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rscg '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 360
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 64
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rscg     '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 360
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... finish workspace base addresses
c
  170 jb2 = ib1+nb
      jb3 = jb2+nb
      jb4 = jb3+nb
      ib5 = n+3*nb+1
      nr = n-nb
      nrp1 = nr+1
c
c ... check for sufficient workspace.
c
      iparm(8) = n+3*nb+4*itmax
      if (nw.ge.iparm(8)) go to 190
      ier = 62
      if (level.ge.0) write (nout,180) nw,iparm(8)
  180 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rscg '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 340
c
  190 continue
      if (level.le.2) go to 220
      write (nout,200)
  200 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
      if (adapt) write (nout,210)
  210 format (1x,'cme is the estimate of the largest eigenvalue of',
     *   ' the jacobi matrix')
  220 if (iparm(11).eq.0) timi1 = timer(0.0)
c
c ... initialize forward pseudo-residual
c
      if (n.gt.1) go to 230
      u(1) = rhs(1) 
      go to 300
  230 do 240 i = 1,nw
         wksp(i) = 0.0
  240 continue
      call scopy (nr,rhs,1,wksp,1)
      call prsred (nr,ndim,maxnz,jcoef,coef,u,wksp)
      call scopy (nb,rhs(nrp1),1,wksp(ib2),1)
      call prsblk (nb,nr,ndim,maxnz,jcoef,coef,wksp,wksp(jb2))
      do 250 i = 1,nb
         wksp(n+i) = wksp(n+i)-u(nr+i)
  250 continue
c
c ... iteration sequence
c
      itmax1 = itmax+1
      do 270 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 260
c
c ... code for the even iterations.
c
c     u           = u(in)             wksp(jb2+nr) = d(in)
c     wksp(ib1)   = u(in-1)           wksp(jb3+nr) = d(in-1)
c
         call itrscg (n,nb,ndim,maxnz,jcoef,coef,u,wksp,wksp(jb2),
     *      wksp(jb3),wksp(jb4),wksp(ib5))
c
         if (halt) go to 300
         go to 270
c
c ... code for the odd iterations.
c
c     u           = u(in-1)           wksp(jb2+nr) = d(in-1)
c     wksp(ib1)   = u(in)             wksp(jb3+nr) = d(in)
c
  260    call itrscg (n,nb,ndim,maxnz,jcoef,coef,wksp,u,wksp(jb3),
     *      wksp(jb2),wksp(jb4),wksp(ib5))
c
         if (halt) go to 300
  270 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 280
      timi2 = timer(0.0)
      time1 = timi2-timi1
  280 ier = 63
      if (level.ge.1) write (nout,290) itmax
  290 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine rscg'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 330
c
c ... method has converged
c
  300 if (iparm(11).ne.0) go to 310
      timi2 = timer(0.0)
      time1 = timi2-timi1
  310 if (level.ge.1) write (nout,320) in
  320 format (/1x,'rscg  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  330 continue
      if (n.eq.1) go to 340
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
      call scopy (nr,rhs,1,u,1)
      call prsred (nr,ndim,maxnz,jcoef,coef,u,u)
c
c ... un-permute matrix,rhs, and solution
c
  340 if (iparm(9).ne.-1) go to 350
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp,iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.62) go to 360
c
c ... optional error analysis 
c
  350 idgts = iparm(12)
      if (idgts.lt.0) go to 360
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  360 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      iparm(8) = iparm(8)-4*(itmax-in)
      if (iparm(11).ne.0) go to 370
      timj2 = timer(0.0)
      time2 = timj2-timj1
  370 if (iparm(3).ne.0) go to 380
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  380 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,390)
  390 format (/1x,'execution successful')
c
      return
      end 
      subroutine rssi (n,ndim,maxnz,jcoef,coef,rhs,u,iwksp,nw,wksp,
     *                 iparm,rparm,ier) 
c
c     itpackv 2d main routine  rssi  (reduced system semi-iterative)
c     each of the main routines --
c           jcg, jsi, sor, ssorcg, ssorsi, rscg, rssi
c     can be used independently of the others
c
c ... function --
c
c          rssi drives the reduced system si algorithm.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for sparse matrix representation.
c          coef   array for sparse matrix representation.
c                 jcoef and coef use the ellpack data structure.
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 initial guess to the solution. on output, it contains
c                 the latest estimate to the solution.
c          iwksp  integer vector workspace of length 3*n
c          nw     input integer.  length of available wksp.  on output,
c                 iparm(8) is amount used.
c          wksp   vector used for working space.  rssi
c                 needs this to be in length at least  n + nb
c                 nb is the order of the black subsystem
c          iparm  integer vector of length 12.  allows user to specify
c                 some integer parameters which affect the method.
c          rparm  vector of length 12. allows user to specify some
c                 parameters which affect the method.
c          ier    output integer.  error flag.
c
c ... rssi module references --
c
c         from itpackv    cheby , chgsi , dfault, echall,
c                         iterm , itrssi, parsi , 
c                         permat, perror, pervec, pmult , prbndx,
c                         prsblk, prsred, pstop , sbelm ,
c                         scal  , sum3  , tstchg, 
c                         unscal,
c                         vout
c          system         abs, alog10, amax0, amax1, float, mod,
c                         sqrt
c
c ... local itpackv references --
c
c          echall, itrssi, permat,
c          perror, pervec, prbndx, prsred, sbelm , scal  ,
c          unscal
c
c     version -  itpackv 2d (january 1990) 
c
c     code written by - david kincaid, roger grimes, john respess
c                       center for numerical analysis
c                       university of texas
c                       austin, tx  78712
c                       (512) 471-1242
c
c     for additional details on the
c          (a) routine    see toms article 1982
c          (b) algorithm  see cna report 150
c
c     based on theory by - david young, david kincaid, lou hageman
c
c     reference the book - applied iterative methods
c                          l. hageman, d. young
c                          academic press, 1981
c
c     **************************************************
c     *               important note                   *
c     *                                                *
c     *      when installing itpackv routines on a     *
c     *  different computer, reset some of the values  *
c     *  in  subroutne dfault.   most important are    *
c     *                                                *
c     *   srelpr      machine relative precision       *
c     *   rparm(1)    stopping criterion               *
c     *                                                *
c     *   also change system-dependent routine         *
c     *   second used in timer                         *
c     *                                                *
c     **************************************************
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwksp(1),iparm(12)
      dimension coef(ndim,1),rhs(n),u(n),wksp(nw),rparm(12) 
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c ... variables in common block - itcom1
c
c     in     - iteration number
c     is     - iteration number when parameters last changed
c     isym   - symmetric/nonsymmetric case switch 
c     itmax  - maximum number of iterations allowed
c     level  - level of output control switch
c     nout   - output unit number
c
c ... variables in common block - itcom2
c
c     adapt  - fully adaptive procedure switch
c     betadt - switch for adaptive determination of beta
c     caseii - adaptive procedure case switch
c     halt   - stopping test switch
c     partad - partially adaptive procedure switch
c
c ... variables in common block - itcom3
c
c     bdelnm - two norm of b times delta-super-n
c     betab  - estimate for the spectral radius of lu matrix
c     cme    - estimate of largest eigenvalue
c     delnnm - inner product of pseudo-residual at iteration n
c     delsnm - inner product of pseudo-residual at iteration s
c     ff     - adaptive procedure damping factor
c     gamma  - acceleration parameter
c     omega  - overrelaxation parameter for sor and ssor
c     qa     - pseudo-residual ratio
c     qt     - virtual spectral radius
c     rho    - acceleration parameter
c     rrr    - adaptive parameter
c     sige   - parameter sigma-sub-e
c     sme    - estimate of smallest eigenvalue
c     specr  - spectral radius estimate for ssor
c     srelpr - machine relative precision
c     stptst - stopping parameter
c     udnm   - two norm of u
c     zeta   - stopping criterion
c
c ... initialize common blocks
c
      ier = 0
      level = iparm(2)
      nout = iparm(4)
      if (level.ge.1) write (nout,10)
   10 format (///1x,'i t p a c k      r s s i    ')
      if (iparm(1).le.0) go to 360
      if (iparm(11).eq.0) timj1 = timer(0.0)
      call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm,rparm,1)
      temp = 500.0*srelpr
      if (zeta.ge.temp) go to 30
      if (level.ge.1) write (nout,20) zeta,srelpr,temp
   20 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine rssi'/1x,
     *   '    rparm(1) =',e10.3,' (zeta)'/1x,
     *   '    a value this small may hinder convergence '/1x,
     *   '    since machine precision srelpr =',e10.3/1x,
     *   '    zeta reset to ',e10.3)
      zeta = temp
   30 continue
      time1 = rparm(9)
      time2 = rparm(10)
      digit1 = rparm(11)
      digit2 = rparm(12)
c
c ... verify n
c
      if (n.gt.0) go to 50
      ier = 71
      if (level.ge.0) write (nout,40) n 
   40 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rssi '/1x,
     *   '    invalid matrix dimension, n =',i8)
      go to 360
c
c ... scale linear system, u, and rhs by the square root of the
c ... diagonal elements.
c
   50 continue
      call scal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,ier)
      if (ier.eq.0) go to 70
      if (level.ge.0) write (nout,60) ier
   60 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rssi '/1x,
     *   '    error detected in routine  scal  '/1x,
     *   '    which scales the system   '/1x,
     *   '    ier = ',i5)
      go to 360
c
c ... remove rows and columns if requested
c
   70 continue
      if (iparm(10).eq.0) go to 80
      tol = rparm(8)
      call sbelm (n,ndim,maxnz,jcoef,coef,rhs,wksp,tol)
c
c ... initialize wksp base addresses.
c
   80 ib1 = 1
      ib2 = ib1+n
      ib3 = ib2+n
c
c ... permute to  red-black system if requested
c
      nb = iparm(9) 
      if (nb.ge.0) go to 110
      if (nb.le.-2) go to 170 
      n3 = n*3
      do 90 i = 1,n3
         iwksp(i) = 0
   90 continue
      call prbndx (n,ndim,maxnz,jcoef,iwksp,iwksp(ib2),nb,level,nout, 
     *   ier)
      if (ier.eq.0) go to 110 
      if (level.ge.0) write (nout,100) ier,nb
  100 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rssi '/1x,
     *   '    error detected in routine  prbndx'/1x,
     *   '    which computes the red-black indexing'/1x,
     *   '    ier = ',i5,' iparm(9) = ',i5,' (nb)')
      go to 340
  110 if (nb.ge.0.and.nb.le.n) go to 130
      ier = 74
      if (level.ge.0) write (nout,120) ier,nb
  120 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rssi     '/1x,
     *   '    error detected in red-black subsystem index'/1x,
     *   '    ier = ',i5,' iparm(9) =',i5,' (nb)')
      go to 340
  130 if (nb.ne.0.and.nb.ne.n) go to 150
      nbo = nb
      nb = n/2
      if (level.ge.2) write (nout,140) nbo,nb
  140 format (/10x,' nb = ',i5,' implies matrix is diagonal'/10x,
     *   ' nb reset to ',i5)
c
c ... permute matrix and rhs
c
  150 if (level.ge.2) write (nout,160) nb
  160 format (/10x,'order of black subsystem = ',i5,' (nb)')
      if (iparm(9).ge.0) go to 170
      call permat (n,ndim,maxnz,jcoef,coef,iwksp,wksp,iwksp(ib3))
      call pervec (n,iwksp,rhs,wksp)
      call pervec (n,iwksp,u,wksp)
c
c ... finish workspace base addresses
c
  170 jb2 = ib1+nb
      nr = n-nb
c
c ... check for sufficient workspace.
c
      iparm(8) = n+nb
      if (nw.ge.iparm(8)) go to 190
      ier = 72
      if (level.ge.0) write (nout,180) nw,iparm(8)
  180 format (/1x,'*** f a t a l     e r r o r ************'//1x,
     *   '    called from itpackv routine rssi '/1x,
     *   '    not enough workspace at ',i10/1x,
     *   '    set iparm(8) =',i10,' (nw)')
      go to 320
c
  190 continue
      if (level.le.2) go to 210
      write (nout,200)
  200 format (///1x,'in the following, rho and gamma are',
     *   ' acceleration parameters')
  210 if (iparm(11).eq.0) timi1 = timer(0.0)
      do 220 i = 1,nw
         wksp(i) = 0.0
  220 continue
c
c ... iteration sequence
c
      if (n.gt.1) go to 230
      u(1) = rhs(1) 
      go to 280
c
  230 itmax1 = itmax+1
      do 250 loop = 1,itmax1
         in = loop-1
         if (mod(in,2).eq.1) go to 240
c
c ... code for the even iterations.
c
c     u           = u(in)
c     wksp(ib1)   = u(in-1)
c
         call itrssi (n,nb,ndim,maxnz,jcoef,coef,rhs,u,wksp,wksp(jb2))
c
         if (halt) go to 280
         go to 250
c
c ... code for the odd iterations.
c
c     u           = u(in-1)
c     wksp(ib1)   = u(in)
c
  240    call itrssi (n,nb,ndim,maxnz,jcoef,coef,rhs,wksp,u,wksp(jb2))
c
         if (halt) go to 280
  250 continue
c
c ... itmax has been reached
c
      if (iparm(11).ne.0) go to 260
      timi2 = timer(0.0)
      time1 = timi2-timi1
  260 ier = 73
      if (level.ge.1) write (nout,270) itmax
  270 format (/1x,'*** w a r n i n g ************'//1x,
     *   '    in itpackv routine rssi'/1x,
     *   '    failure to converge in',i5,' iterations')
      if (iparm(3).eq.0) rparm(1) = stptst
      go to 310
c
c ... method has converged
c
  280 if (iparm(11).ne.0) go to 290
      timi2 = timer(0.0)
      time1 = timi2-timi1
  290 if (level.ge.1) write (nout,300) in
  300 format (/1x,'rssi  has converged in ',i5,' iterations')
c
c ... put solution into u if not already there.
c
  310 continue
      if (n.eq.1) go to 320
      if (mod(in,2).eq.1) call scopy (n,wksp,1,u,1)
      call scopy (nr,rhs,1,u,1)
      call prsred (nr,ndim,maxnz,jcoef,coef,u,u)
c
c ... un-permute matrix,rhs, and solution
c
  320 if (iparm(9).ne.-1) go to 330
      call permat (n,ndim,maxnz,jcoef,coef,iwksp(ib2),wksp,iwksp(ib3))
      call pervec (n,iwksp(ib2),rhs,wksp(ib2))
      call pervec (n,iwksp(ib2),u,wksp(ib2))
      if (ier.eq.72) go to 340
c
c ... optional error analysis 
c
  330 idgts = iparm(12)
      if (idgts.lt.0) go to 340
      if (iparm(2).le.0) idgts = 0
      call perror (n,ndim,maxnz,jcoef,coef,rhs,u,wksp,digit1,digit2,
     *   idgts)
c
c ... unscale the matrix, solution, and rhs vectors.
c
  340 continue
      call unscal (n,ndim,maxnz,jcoef,coef,rhs,u,wksp)
c
c ... set return parameters in iparm and rparm
c
      if (iparm(11).ne.0) go to 350
      timj2 = timer(0.0)
      time2 = timj2-timj1
  350 if (iparm(3).ne.0) go to 360
      iparm(1) = in 
      iparm(9) = nb 
      rparm(2) = cme
      rparm(3) = sme
      rparm(9) = time1
      rparm(10) = time2
      rparm(11) = digit1
      rparm(12) = digit2
c
  360 continue
      if (level.ge.3) call echall (n,ndim,maxnz,jcoef,coef,rhs,iparm, 
     *   rparm,2)
      if (ier.eq.0.and.level.ge.1) write (nout,370)
  370 format (/1x,'execution successful')
c
      return
      end 
      subroutine itjcg (n,ndim,maxnz,jcoef,coef,u,u1,d,d1,dtwd,tri)
c
c ... itjcg performs one iteration of the jacobi conjugate gradient
c     algorithm.  it is called by jcg.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer sparse matrix representation
c          coef   sparse matrix representation
c          u      input vector.  contains the value of the
c                 solution vector at the end of in iterations.
c          u1     input/output vector.  on input, it contains
c                 the value of the solution at the end of the in-1
c                 iteration.  on output, it will contain the newest
c                 estimate for the solution vector.
c          d      input vector.  contains the pseudo-residual
c                 vector after in iterations.
c          d1     input/output vector.  on input, d1 contains
c                 the pseudo-residual vector after in-1 iterations.  on
c                 output, it will contain the newest pseudo-residual
c                 vector.
c          dtwd   array.  used in the computations of the
c                 acceleration parameter gamma and the new pseudo-
c                 residual.
c          tri    array.  stores the tridiagonal matrix associated
c                 with the eigenvalues of the conjugate gradient
c                 polynomial. 
c
c ... local itpackv references --
c
c          chgcon, iterm , parcon, pjac  , pstop ,
c          sum3
c
c ... specifications for arguments
c
      integer jcoef(ndim,1)
      dimension u(n),u1(n),d(n),d1(n),dtwd(n),tri(2,1),coef(ndim,1)
c
c ... specifications for local variables
c
      logical q1
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine jcg
c
c ... compute new estimate for cme if adapt = .true.
c
      save
      if (adapt) call chgcon (itmax,tri,gamold,rhoold,1)
c
c ... test for stopping
c
      delnnm = sdot(n,d,1,d,1)
      dnrm = delnnm 
      con = cme
      call pstop (n,u,dnrm,con,1,q1)
      if (halt) go to 50
c
c ... compute rho and gamma - acceleration parameters
c
      do 10 i = 1,n 
         dtwd(i) = 0.0
   10 continue
      call pjac (n,ndim,maxnz,jcoef,coef,d,dtwd)
      dtnrm = sdot(n,d,1,dtwd,1)
      if (isym.eq.0) go to 20 
      rhotmp = sdot(n,dtwd,1,d1,1)
      call parcon (dtnrm,c1,c2,c3,c4,gamold,rhotmp,1)
      rhoold = rhotmp
      go to 30
   20 call parcon (dtnrm,c1,c2,c3,c4,gamold,rhoold,1)
c
c ... compute u(in+1) and d(in+1)
c
   30 do 40 i = 1,n 
         u1(i) = c1*d(i)+c2*u(i)+c3*u1(i)
         d1(i) = c1*dtwd(i)+c4*d(i)+c3*d1(i)
   40 continue
c
c ... output intermediate information
c
   50 call iterm (n,coef,u,dtwd,1)
c
      return
      end 
      subroutine itjsi (n,ndim,maxnz,jcoef,coef,rhs,u,u1,d,icnt)
c
c ... function --
c
c          itjsi performs one iteration of the
c          jacobi semi-iterative algorithm.  it is called by jsi.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                  routine
c          maxnz  maximum number of nonzero entries per row 
c          jcoef  integer array for sparse matrix representation
c          coef   array for sparse matrix representation
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input vector.  contains the estimate for the
c                 solution vector after in iterations.
c          u1     input/output vector.  on input, u1 contains the
c                 solution vector after in-1 iterations.  on output,
c                 it will contain the newest estimate for the solution
c                 vector.
c          d      array.  d is used for the computation of the
c                 pseudo-residual array for the current iteration.
c          icnt   number of iterations since last change of sme
c
c ... local itpackv references --
c
c          chgsi , chgsme, iterm , parsi , pjac  , pstop , pvtbv ,
c          sum3  , tstchg
c
c ... specifications for arguments
c
      integer jcoef(ndim,1)
      dimension coef(ndim,1),rhs(1),u(1),u1(1),d(1)
c
c ... specifications for local variables
c
      logical q1
c
c ... specifications for function subprograms
c
      logical tstchg,chgsme
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine jsi
c
      if (in.eq.0) icnt = 0
c
c ... compute pseudo-residuals
c
      call scopy (n,rhs,1,d,1)
      call pjac (n,ndim,maxnz,jcoef,coef,u,d)
      do 10 i = 1,n 
         d(i) = d(i)-u(i)
   10 continue
c
c ... stopping and adaptive change tests
c
      oldnrm = delnnm
      delnnm = sdot(n,d,1,d,1)
      dnrm = delnnm 
      con = cme
      call pstop (n,u,dnrm,con,1,q1)
      if (halt) go to 70
      if (.not.adapt) go to 50
      if (.not.tstchg(1)) go to 20
c
c ... change iterative parameters (cme) 
c
      dtnrm = pvtbv(n,ndim,maxnz,jcoef,coef,d)
      call chgsi (dtnrm,1)
      if (.not.adapt) go to 50
      go to 30
c
c ... test if sme needs to be changed and change if necessary.
c
   20 continue
      if (caseii) go to 50
      if (.not.chgsme(oldnrm,icnt)) go to 50
      icnt = 0
c
c ... compute u(in+1) after change of parameters
c
   30 do 40 i = 1,n 
         u1(i) = u(i)+gamma*d(i)
   40 continue
      go to 70
c
c ... compute u(in+1) without change of parameters
c
   50 call parsi (c1,c2,c3,1) 
      do 60 i = 1,n 
         u1(i) = c1*d(i)+c2*u(i)+c3*u1(i)
   60 continue
c
c ... output intermediate information
c
   70 call iterm (n,coef,u,d,2)
c
      return
      end 
      subroutine itsor (n,nb,ndim,maxnz,jcoef,coef,rhs,u,nup,
     *                  moved,wk,iwk)
c
c          itsor performs one iteration of the successive
c          overrelaxation algorithm.  it is called by sor.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          nb     dimension of the black subsystem
c          ndim   row dimension of coef and jcoef arrays in calling
c                  routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array of coefficient columns
c          coef   array of coefficients 
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input/output vector.  on input, u contains the
c                 solution vector after in iterations.  on output,
c                 it will contain the newest estimate for the solution
c                 vector.
c          nup    maximum number of nonzeros per row in upper
c                  triangle u (for natural ordering only)
c          moved  logical variable indicating whether or not matrix
c                  was reshuffled
c          wk     array.  work vector of length 2*n.
c          iwk    integer array giving wave populations.
c
c ... local itpackv references --
c
c          ipstr , iterm , pfsor1, pstop , tau
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwk(1)
      dimension coef(ndim,1),rhs(1),u(1),wk(1)
      logical moved 
c
c ... specifications for local variables
c
      logical change,q1
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine sor
c
c ... set initial parameters not already set
c
      save
      ib2 = n+1
      if (in.ne.0) go to 20
      call pstop (n,u,0.0,0.0,0,q1)
      if (adapt) go to 10
      change = .false.
      ip = 0
      iphat = 2
      iss = 0
      call sorscl (n,ndim,maxnz,coef,rhs,omega)
      go to 30
c
   10 change = .true.
      ip = 0
      omegap = omega
      omega = 1.0
      iss = 0
      iphat = 2
      ipstar = 4
      if (omegap.le.1.0) change = .false.
c
c ... reset omega, iphat, and ipstar (circle a in flowchart)
c
   20 if (.not.change) go to 30
      change = .false.
      is = is+1
      ip = 0
      iss = 0
      omgsav = omega
      omega = amin1(omegap,tau(is))
      iphat = max0(3,int((omega-1.0)/(2.0-omega)))
      ipstar = ipstr(omega)
      factor = omega/omgsav
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
c
c ... compute u (in + 1) and norm of del(s,p) - circle b in flow chart
c
   30 continue
      delsnm = delnnm
      spcrm1 = specr
      call scopy (n,rhs,1,wk,1)
      call pfsor1 (n,nb,ndim,maxnz,jcoef,coef,u,wk,nup,moved,
     *   wk(ib2),iwk)
      if (delnnm.eq.0.0) go to 40
      if (in.ne.0) specr = delnnm/delsnm
      if (ip.lt.iphat) go to 80
c
c ... stopping test, set h
c
      if (specr.ge.1.0) go to 80
      if (.not.(specr.gt.(omega-1.0))) go to 40
      h = specr
      go to 50
   40 iss = iss+1
      h = omega-1.0 
c
c ... perform stopping test.
c
   50 continue
      dnrm = delnnm**2
      call pstop (n,u,dnrm,h,1,q1)
      if (.not.halt) go to 60 
      factor = 1.0/omega
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
      go to 80
c
c ... method has not converged yet, test for changing omega 
c
   60 if (.not.adapt) go to 80
      if (ip.lt.ipstar) go to 80
      if (omega.gt.1.0) go to 70
      cme = sqrt(abs(specr))
      omegap = 2.0/(1.0+sqrt(abs(1.0-specr)))
      change = .true.
      go to 80
   70 if (iss.ne.0) go to 80
      if (specr.le.(omega-1.0)**ff) go to 80
      if ((specr+0.00005).le.spcrm1) go to 80
c
c ... change parameters
c
      cme = (specr+omega-1.0)/(sqrt(abs(specr))*omega)
      omegap = 2.0/(1.0+sqrt(abs(1.0-cme*cme)))
      change = .true.
c
c ... output intermediate information
c
   80 call iterm (n,coef,u,wk,3)
      ip = ip+1
c
      return
      end 
      subroutine itsrcg (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,u,
     *                   u1,c,c1,d,dl,moved,wk,tri,iwk)
c
c ... function --
c
c          itsrcg performs one iteration of the symmetric sor
c          conjugate gradient algorithm.  it is called by ssorcg.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          nb     order of the black subsystem
c          nup    maximum number of nonzeros per row in the upper
c                  triangle u  (input for natural ordering only)
c          nlow   maximum number of nonzeros per row in the lower
c                  triangle l  (input for natural ordering only)
c          ndim   row dimension of jcoef and coef arrays in calling
c                  routine
c          maxnz  maximum number of nonzeros per row in the entire
c                  array
c          jcoef  integer array for coefficient columns
c          coef   array for coefficients
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input vector.  contains the estimate of the
c                 solution vector after in iterations.
c          u1     input/output vector.  on input, u1 contains the
c                 the estimate for the solution after in-1 iterations.
c                 on output, u1 contains the updated estimate.
c          c      input vector.  contains the forward residual
c                 after in iterations.
c          c1     input/output vector.  on input, c1 contains
c                 the forward residual after in-1 iterations.  on
c                 output, c1 contains the updated forward residual.
c          d      vector.  is used to compute the backward pseudo-
c                 residual vector for the current iteration.
c          dl     vector.  is used in the computations of the
c                 acceleration parameters.
c          moved  logical constant indicating whether or not matrix
c                  has been reshuffled
c          wk     workspace vector of length n.
c          tri    vector. stores the tridiagonal matrix associated
c                 with the conjugate gradient acceleration. 
c          iwk    integer vector giving wave population counts
c
c ... local itpackv references --
c
c        chgcon, iterm , omeg  , omgchg, omgstr, parcon, pbeta ,
c        pbsor , pfsor , pjac  , pstop , pvtbv ,
c        sum3
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwk(1)
      dimension coef(ndim,1),rhs(1),u(1),u1(1),c(1),c1(1),d(1),dl(1), 
     *   wk(1),tri(2,1)
c
c ... specifications for local variables
c
      logical q1,moved
c
c ... specifications for function subprograms
c
      logical omgchg,omgstr
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine ssorcg
c
c ... calculate s-prime for adaptive procedure.
c
      save
      if (adapt.or.partad) call chgcon (itmax,tri,gamold,rhoold,3)
c
c ... compute backward residual
c
      call scopy (n,rhs,1,wk,1)
      do 10 i = 1,n 
         d(i) = c(i)+u(i)
   10 continue
      call pbsor (n,nb,ndim,maxnz,jcoef,coef,d,wk,nup,nlow,moved,iwk) 
      do 20 i = 1,n 
         d(i) = d(i)-u(i)
   20 continue
c
c ... compute acceleration parameters and then u(in+1) (in u1)
c
      call scopy (n,d,1,dl,1) 
      do 30 i = 1,n 
         wk(i) = 0.0
   30 continue
      call pfsor (n,nb,ndim,maxnz,jcoef,coef,dl,wk,nup,moved,iwk)
      do 40 i = 1,n 
         dl(i) = d(i)-dl(i)
   40 continue
      delnnm = sdot(n,c,1,c,1)
      if (delnnm.eq.0.0) go to 80
      dnrm = sdot(n,c,1,dl,1) 
      if (dnrm.eq.0.0) go to 80
      if (isym.eq.0) go to 50 
      rhotmp = sdot(n,c,1,c1,1)-sdot(n,dl,1,c1,1) 
      call parcon (dnrm,t1,t2,t3,t4,gamold,rhotmp,3)
      rhoold = rhotmp
      go to 60
   50 call parcon (dnrm,t1,t2,t3,t4,gamold,rhoold,3)
   60 do 70 i = 1,n 
         u1(i) = t1*d(i)+t2*u(i)+t3*u1(i)
   70 continue
c
c ... test for stopping
c
   80 bdelnm = sdot(n,d,1,d,1)
      dnrm = bdelnm 
      con = specr
      call pstop (n,u,dnrm,con,1,q1)
      if (halt) go to 190
c
c ... if non- or partially-adaptive, compute c(in+1) and exit.
c
      if (adapt) go to 100
      do 90 i = 1,n 
         c1(i) = -t1*dl(i)+t2*c(i)+t3*c1(i)
   90 continue
      go to 190
c
c ... fully adaptive procedure
c
  100 continue
      omgsav = omega
      if (omgstr(1)) go to 170
      if (omgchg(1)) go to 120
c
c ... parameters have been unchanged.  compute c(in+1) and exit.
c
      do 110 i = 1,n
         c1(i) = -t1*dl(i)+t2*c(i)+t3*c1(i)
  110 continue
      go to 190
c
c ... it has been decided to change parameters
c        (1) compute new betab if betadt = .true. 
c
  120 continue
      if (.not.betadt) go to 130
      betnew = pbeta(n,nb,ndim,maxnz,jcoef,coef,d,nup,nlow,moved,wk,c1)/
     *   (bdelnm*omega**2)
      betab = amax1(betab,0.25,betnew)
c
c ...    (2) compute new cme, omega, and specr
c
  130 continue
      if (caseii) go to 140
      dnrm = pvtbv(n,ndim,maxnz,jcoef,coef,d)/omega
      go to 160
  140 do 150 i = 1,n
         wk(i) = 0.0
  150 continue
      call pjac (n,ndim,maxnz,jcoef,coef,d,wk)
      dnrm = sdot(n,wk,1,wk,1)/(omega**2)
  160 call omeg (dnrm,3)
c
c ...    (3) compute new forward residual since omega has been changed.
c
  170 continue
      factor = omega/omgsav
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
      call scopy (n,rhs,1,wk,1)
      call scopy (n,u1,1,c1,1)
      call pfsor (n,nb,ndim,maxnz,jcoef,coef,c1,wk,nup,moved,iwk)
      do 180 i = 1,n
         c1(i) = c1(i)-u1(i)
  180 continue
c
c ... output intermediate results.
c
  190 call iterm (n,coef,u,wk,4)
c
      return
      end 
      subroutine itsrsi (n,nb,nup,nlow,ndim,maxnz,jcoef,coef,rhs,u,
     *                   u1,c,d,ctwd,moved,wk,iwk)
c
c ... function --
c
c          itsrsi performs one iteration of the symmetric sor
c          semi-iteration algorithm.  it is called by ssorsi.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          nb     order of the black subsystem (for red-black order-
c                   ing only) 
c          nup    maximum number of nonzeros per row in the upper
c                   triangle u (for natural ordering only)
c          nlow   maximum number of nonzeros per row in the lower
c                   triangle l (for natural ordering only)
c          ndim   row dimension of coef and jcoef arrays in calling
c                   routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer array for coefficient columns
c          coef   array for coefficients
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          u      input vector.  contains the estimate of the
c                 solution vector after in iterations.
c          u1     input/output vector.  on input, u1 contains the
c                 the estimate for the solution after in-1 iterations.
c                 on output, u1 contains the updated estimate.
c          c      vector.  is used to compute the forward pseudo-
c                 residual vector for the current iteration.
c          d      vector.  is used to compute the backward pseudo-
c                 residual vector for the current iteration.
c          ctwd   vector.  is used in the computations of the
c                 acceleration parameters.
c          moved  logical constant indicating whether or not matrix
c                  has been reshuffled
c          wk     workspace vector of length n.
c          iwk    integer vector giving wave population counts
c
c ... local itpackv references --
c
c         chgsi , iterm , omeg  , omgstr, parsi , pbeta , pbsor ,
c         pfsor , pjac  , pstop , pvtbv , sum3  , 
c         tstchg
c
c ... specifications for arguments
c
      integer jcoef(ndim,1),iwk(1)
      dimension coef(ndim,1),rhs(1),u(1),u1(1),c(1),d(1),ctwd(1),wk(1)
c
c ... specifications for local variables
c
      logical q1,moved
c
c ... specifications for function subprograms
c
      logical omgstr,tstchg
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine ssorsi
c
c ... compute pseudo-residuals (forward and backward)
c
      call scopy (n,rhs,1,wk,1)
      call scopy (n,u,1,ctwd,1)
      call pfsor (n,nb,ndim,maxnz,jcoef,coef,ctwd,wk,nup,moved,iwk)
      do 10 i = 1,n 
         c(i) = ctwd(i)-u(i)
   10 continue
      call scopy (n,rhs,1,wk,1)
      call pbsor (n,nb,ndim,maxnz,jcoef,coef,ctwd,wk,nup,nlow,moved,iwk)
      do 20 i = 1,n 
         d(i) = ctwd(i)-u(i)
   20 continue
c
c ... compute u(in+1) -- contained in the vector u1.
c
      call parsi (c1,c2,c3,3) 
      do 30 i = 1,n 
         u1(i) = c1*d(i)+c2*u(i)+c3*u1(i)
   30 continue
c
c ... test for stopping
c
      bdelnm = sdot(n,d,1,d,1)
      dnrm = bdelnm 
      con = specr
      call pstop (n,u,dnrm,con,1,q1)
      if (halt.or..not.(adapt.or.partad)) go to 100
c
c ... adaptive procedure
c
      omgsav = omega
      if (omgstr(1)) go to 100
      delnnm = sdot(n,c,1,c,1)
      if (in.eq.is) delsnm = delnnm
      if (in.eq.0.or..not.tstchg(1)) go to 100
c
c ... it has been decided to change parameters.
c ...    (1) compute ctwd
c
      call scopy (n,d,1,ctwd,1)
      do 40 i = 1,n 
         wk(i) = 0.0
   40 continue
      call pfsor (n,nb,ndim,maxnz,jcoef,coef,ctwd,wk,nup,moved,iwk)
      do 50 i = 1,n 
         ctwd(i) = ctwd(i)+c(i)-d(i)
   50 continue
c
c ...    (2) compute new spectral radius for current omega. 
c
      dnrm = sdot(n,c,1,ctwd,1)
      call chgsi (dnrm,3)
      if (.not.adapt) go to 100
c
c ...    (3) compute new betab if betadt = .true. 
c
      if (.not.betadt) go to 60
      betnew = pbeta(n,nb,ndim,maxnz,jcoef,coef,d,nup,nlow,moved,wk,
     *   ctwd)/(bdelnm*omega**2)
      betab = amax1(betab,0.25,betnew)
c
c ...    (4) compute new cme, omega, and specr.
c
   60 continue
      if (caseii) go to 70
      dnrm = pvtbv(n,ndim,maxnz,jcoef,coef,d)/omega
      go to 90
   70 do 80 i = 1,n 
         wk(i) = 0.0
   80 continue
      call pjac (n,ndim,maxnz,jcoef,coef,d,wk)
      dnrm = sdot(n,wk,1,wk,1)/(omega**2)
   90 call omeg (dnrm,3)
      factor = omega/omgsav
      call sorscl (n,ndim,maxnz,coef,rhs,factor)
c
c ... output intermediate information
c
  100 call iterm (n,coef,u,wk,5)
c
      return
      end 
      subroutine itrscg (n,nb,ndim,maxnz,jcoef,coef,ub,ub1, 
     *                   db,db1,wb,tri) 
c
c ... itrscg performs one iteration of the reduced system conjugate
c ... gradient algorithm.  it is called by rscg.
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          nb     input integer.  contains the number of black points 
c                 in the red-black matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                 routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer data structure for coefficient columns
c          coef   data structure for matrix coefficients
c          ub     input vector.  contains the estimate for the
c                 solution on the black points after in iterations.
c          ub1    input/output vector.  on input, ub1 contains
c                 the solution vector after in-1 iterations.  on
c                 output, it will contain the newest estimate for
c                 the solution vector.  this is only for the black
c                 points.
c          db     input array.  db(nrp1) contains the value of
c                 the current pseudo-residual on the black points.
c          db1    input/output array.  db1(nrp1) contains the
c                 pseudo-residual on the black points for the in-1
c                 iteration.
c                 on input.  on output, it is for the in+1 iteration. 
c          wb     array.  wb is used for computations involving
c                 black vectors.
c          tri    array.  stores the tridiagonal matrix associated
c                 with conjugate gradient acceleration.
c
c ... local itpackv references --
c
c          chgcon, iterm , parcon, prsblk, prsred, pstop ,
c          sum3
c
c ... specifications for arguments
c
      integer jcoef(ndim,1)
      dimension coef(ndim,1),ub(1),ub1(1),db(1),db1(1),wb(1),tri(2,1) 
c
c ... specifications for local variables
c
      logical q1
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine rscg
c
c ... compute new estimate for cme if adapt = .true.
c
      save
      nr = n-nb
      nrp1 = nr+1
      if (adapt) call chgcon (itmax,tri,gamold,rhoold,2)
c
c ... test for stopping
c
      delnnm = sdot(nb,db(nrp1),1,db(nrp1),1)
      dnrm = delnnm 
      con = cme
      call pstop (nb,ub(nrp1),dnrm,con,2,q1)
      if (halt) go to 70
c
c ... compute acceleration parameters
c
      do 10 i = 1,nr
         ub1(i) = 0.0
   10 continue
      call prsred (nr,ndim,maxnz,jcoef,coef,db,ub1)
      do 20 i = nrp1,n
         wb(i) = 0.0
   20 continue
      call prsblk (nb,nr,ndim,maxnz,jcoef,coef,ub1,wb)
      dnrm = sdot(nb,db(nrp1),1,wb(nrp1),1)
      if (isym.eq.0) go to 30 
      rhotmp = sdot(nb,wb(nrp1),1,db1(nrp1),1)
      call parcon (dnrm,c1,c2,c3,c4,gamold,rhotmp,2)
      rhoold = rhotmp
      go to 40
   30 call parcon (dnrm,c1,c2,c3,c4,gamold,rhoold,2)
c
c ... compute ub(in+1) and db(in+1)
c
   40 do 50 i = nrp1,n
         ub1(i) = c1*db(i)+c2*ub(i)+c3*ub1(i)
   50 continue
      do 60 i = nrp1,n
         db1(i) = c1*wb(i)+c4*db(i)+c3*db1(i)
   60 continue
c
c ... output intermediate information
c
   70 call iterm (nb,coef(nrp1,1),ub(nrp1),wb(nrp1),6)
      return
c
      end 
      subroutine itrssi (n,nb,ndim,maxnz,jcoef,coef,rhs,ub,ub1,db)
c
c ... itrssi performs one iteration of the reduced system semi-
c ... iteration algorithm.  it is called by rssi. 
c
c ... parameter list --
c
c          n      input integer.  dimension of the matrix.
c          nb     input integer.  contains the number of black points 
c                 in the red-black matrix.
c          ndim   row dimension of jcoef and coef arrays in calling
c                 routine
c          maxnz  maximum number of nonzeros per row
c          jcoef  integer data structure for coefficient columns
c          coef   data structure for matrix coefficients
c          rhs    input vector.  contains the right hand side
c                 of the matrix problem.
c          ub     input vector.  contains the estimate for the
c                 solution on the black points after in iterations.
c          ub1    input/output vector.  on input, ub1 contains
c                 the solution vector after in-1 iterations.  on
c                 output, it will contain the newest estimate for
c                 the solution vector.  this is only for the black
c                 points.
c          db     input array.  db(nrp1) contains the value of
c                 the current pseudo-residual on the black points.
c
c ... local itpackv references --
c
c          chgsi , iterm , parsi , prsblk, prsred, pstop ,
c          sum3  , tstchg
c
c ... specifications for arguments
c
      integer jcoef(ndim,1)
      dimension coef(ndim,1),rhs(1),ub(1),ub1(1),db(1)
c
c ... specifications for local variables
c
      logical q1
c
c ... specifications for function subprograms
c
      logical tstchg
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in routine rssi
c
c ... compute ur(in) into ub
c
      nr = n-nb
      nrp1 = nr+1
      call scopy (nr,rhs,1,ub,1)
      call prsred (nr,ndim,maxnz,jcoef,coef,ub,ub)
c
c ... compute pseudo-residual, db(in)
c
      call scopy (nb,rhs(nrp1),1,db(nrp1),1)
      call prsblk (nb,nr,ndim,maxnz,jcoef,coef,ub,db)
      do 10 i = nrp1,n
         db(i) = db(i)-ub(i)
   10 continue
c
c ... test for stopping
c
      delnnm = sdot(nb,db(nrp1),1,db(nrp1),1)
      dnrm = delnnm 
      const = cme
      call pstop (nb,ub(nrp1),dnrm,const,2,q1)
      if (halt) go to 60
      if (.not.adapt) go to 40
c
c ... test to change parameters
c
      if (.not.tstchg(2)) go to 40
c
c ... change parameters
c
      do 20 i = 1,nr
         ub1(i) = 0.0
   20 continue
      call prsred (nr,ndim,maxnz,jcoef,coef,db,ub1)
      dnrm = sdot(nr,ub1,1,ub1,1)
      call chgsi (dnrm,2)
      if (.not.adapt) go to 40
c
c ... compute ub(n+1) after changing parameters
c
      do 30 i = nrp1,n
         ub1(i) = ub(i)+gamma*db(i)
   30 continue
      go to 60
c
c ... compute ub(n+1) without change of parameters
c
   40 call parsi (c1,c2,c3,2) 
      do 50 i = nrp1,n
         ub1(i) = c1*db(i)+c2*ub(i)+c3*ub1(i)
   50 continue
c
c ... output intermediate information
c
   60 call iterm (nb,coef(nrp1,1),ub(nrp1),db(nrp1),7)
      return
c
      end 
      function cheby (qa,qt,rrr,ip,cme,sme)
c
c ... cheby computes the solution to the chebyshev equation 
c
c ... parameter list --
c
c          qa     ratio of pseudo-residuals
c          qt     virtual spectral radius
c          rrr    adaptive parameter
c          ip     number of iterations since last change of 
c                     parameters
c          cme,   estimates for the largest and smallest eigen-
c          sme      values of the iteration matrix
c
      z = 0.5*(qa+sqrt(abs(qa**2-qt**2)))*(1.0+rrr**ip)
      x = z**(1.0/float(ip))
      y = (x+rrr/x)/(1.0+rrr) 
c
      cheby = 0.5*(cme+sme+y*(2.0-cme-sme))
c
      return
      end 
      subroutine chgcon (ldt,tri,gamold,rhoold,ibmth)
c
c ... chgcon computes the new estimate for the largest eigenvalue
c     for conjugate gradient acceleration.
c
c ... parameter list --
c
c          ldt    leading dimension of tri
c          tri    tridiagonal matrix associated with the eigenvalues
c                    of the conjugate gradient polynomial
c          gamold
c            and
c          rhoold previous values of acceleration parameters
c          ibmth  indicator of basic method being accelerated by cg
c                      ibmth = 1,  jacobi
c                            = 2,  reduced system 
c                            = 3,  ssor 
c
c ... local itpackv references --
c
c          eigvns, eigvss
c
c ... specifications for arguments
c
      dimension tri(ldt,4)
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in main routine
c
      save
      go to (10,20,30), ibmth 
c
c ... jacobi conjugate gradient
c
   10 start = cme
      ip = in
      go to 40
c
c ... reduced system cg
c
   20 start = cme**2
      ip = in
      go to 40
c
c ... ssor cg
c
   30 if (adapt) start = spr
      if (.not.adapt) start = specr
      ip = in-is
c
c ... define the matrix
c
   40 if (ip.ge.2) go to 60
      if (ip.eq.1) go to 50
c
c ... ip = 0
c
      end = 0.0
      cmold = 0.0
      go to 110
c
c ... ip = 1
c
   50 end = 1.0-1.0/gamma
      tri(1,1) = end
      tri(1,2) = 0.0
      go to 110
c
c ... ip > 1
c
   60 if (abs(start-cmold).le.zeta*start) go to 120
      cmold = start 
c
c ... compute the largest eigenvalue
c
      tri(ip,1) = 1.0-1.0/gamma
      tri(ip,2) = (1.0-rho)/(rho*rhoold*gamma*gamold)
      if (isym.ne.0) go to 80 
      end = eigvss(ip,tri,start,zeta,itmax,ier)
      if (ier.eq.0) go to 100 
      if (level.ge.2) write (nout,70) ier
   70 format (/10x,'difficulty in computation of maximum eigenvalue'/ 
     *         15x,'of iteration matrix'/
     *         10x,'routine zbrent returned ier =',i5)
      go to 100
   80 continue
      end = eigvns(ldt,ip,tri,tri(1,3),tri(1,4),ier)
      if (ier.eq.0) go to 100 
      if (level.ge.2) write (nout,90) ier
   90 format (/10x,'difficulty in computation of maximum eigenvalue'/ 
     *         15x,'of iteration matrix'/
     *         10x,'routine eqrt1s returned ier =',i5)
  100 continue
      if (ier.ne.0) go to 130 
c
c ... set spectral radius for the various methods 
c
  110 if (ibmth.eq.1) cme = end
      if (ibmth.eq.2) cme = sqrt(abs(end))
      if (ibmth.eq.3.and.adapt) spr = end
      if (ibmth.eq.3.and..not.adapt) specr = end
      return
c
c ... relative change in cme is less than zeta.  therefore stop
c     changing.
c
  120 adapt = .false.
      partad = .false.
      return
c
c ... estimate for cme > one.  therefore need to stop adaptive
c     procedure and keep old value of cme.
c
  130 adapt = .false.
      partad = .false.
      if (level.ge.2) write (nout,140) in,start
  140 format (/10x,'estimate of maximum eigenvalue of jacobi   '/15x, 
     *   'matrix (cme) not accurate'/10x,
     *   'adaptive procedure turned off at iteration ',i5/10x,
     *   'final estimate of maximum eigenvalue =',e15.7/)
c
      return
      end 
      subroutine chgsi (dtnrm,ibmth)
c
c ... chgsi computes new chebyshev acceleration parameters adaptively.
c
c ... parameter list --
c
c          dtnrm  numerator of rayleigh quotient
c          ibmth  indicator of basic method being accelerated by si
c                      ibmth = 1,   jacobi
c                            = 2,   reduced system
c                            = 3,   symmetric sor 
c
c ... local itpackv references --
c
c          cheby
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in main routine
c
      go to (10,30,50), ibmth 
c
c     ---------------------
c ... jacobi semi-iterative
c     ---------------------
c
c ... chebyshev equation
c
   10 continue
      if (in.eq.0) zm1 = cme
      if (in.ne.0) zm1 = cheby(qa,qt,rrr,in-is,cme,sme)
c
c ... rayleigh quotient
c
      zm2 = dtnrm/delnnm
c
c ... computation of iterative parameters
c
      cmold = cme
      cme = amax1(zm1,zm2,cmold)
      if (cme.ge.1.0) go to 20
      if (caseii) sme = -cme
      sige = (cme-sme)/(2.0-cme-sme)
      gamma = 2.0/(2.0-cme-sme)
      rrr = (1.0-sqrt(abs(1.0-sige*sige)))/(1.0+sqrt(abs(1.0-sige*sige))
     *   )
      is = in
      delsnm = delnnm
      rho = 1.0
      if (level.ge.2) write (nout,90) in,zm1,zm2,cme,gamma,cme
      return
c
c ... adaptive procedure failed for jacobi si
c
   20 cme = cmold
      adapt = .false.
      if (level.ge.2) write (nout,110) in,cme
      return
c
c     -----------------------------
c ... reduced system semi-iterative
c     -----------------------------
c
c ... chebyshev equation
c
   30 continue
      if (in.eq.0) zm1 = cme
      if (in.ne.0) zm1 = cheby(qa,qt,rrr,2*(in-is),0.0,0.0) 
c
c ... rayleigh quotient
c
      zm2 = sqrt(abs(dtnrm/delnnm))
c
c ... computation of new iterative parameters
c
      cmold = cme
      cme = amax1(zm1,zm2,cmold)
      if (cme.ge.1.0) go to 40
      sige = cme*cme/(2.0-cme*cme)
      gamma = 2.0/(2.0-cme*cme)
      rrr = (1.0-sqrt(abs(1.0-cme*cme)))/(1.0+sqrt(abs(1.0-cme*cme))) 
      is = in
      delsnm = delnnm
      rho = 1.0
      if (level.ge.2) write (nout,90) in,zm1,zm2,cme,gamma,cme
      return
c
c ... adaptive procedure failed for reduced system si
c
   40 cme = cmold
      adapt = .false.
      if (level.ge.2) write (nout,110) in,cme
      return
c
c     -----------------------------
c ... symmetric sor semi-iterative
c     ----------------------------
c
   50 continue
      if (specr.eq.0.0) specr = .171572875
      if (in.eq.0) go to 60
      zm1 = cheby(qa,qt,rrr,in-is,specr,0.0)
      go to 70
   60 zm1 = specr
      spr = specr
c
c ... rayleigh quotient
c
   70 zm2 = dtnrm/delnnm
c
c ... computation of new estimate for spectral radius
c
      if (adapt) go to 80
c
c ... partially adaptive ssor si
c
      specr = amax1(zm1,zm2,specr)
      is = in+1
      delsnm = delnnm
      if (level.ge.2) write (nout,100) in,zm1,zm2,cme,specr 
      return
c
c ... fully adaptive ssor si
c
   80 spr = amax1(zm1,zm2,spr)
      return
c
c ... format statements
c
   90 format (/5x,'parameters were changed at iteration no.',i5/10x,
     *   'solution to chebyshev eqn.       =',e15.7/10x,
     *   'solution to rayleigh quotient    =',e15.7/10x,
     *   'new estimate for cme             =',e15.7/10x,
     *   'new estimate for gamma           =',e15.7/10x,
     *   'new estimate for spectral radius =',e15.7/)
c
  100 format (/5x,'parameters were changed at iteration no.',i5/10x,
     *   'solution to chebyshev eqn.       =',e15.7/10x,
     *   'solution to rayleigh quotient    =',e15.7/10x,
     *   'new estimate for cme             =',e15.7/10x,
     *   'new estimate for spectral radius =',e15.7/)
c
  110 format (/10x,'estimate of maximum eigenvalue of jacobi   '/10x, 
     *   'matrix (cme) too large'/10x,
     *   'adaptive procedure turned off at iteration ',i5/10x,
     *   'final estimate of maximum eigenvalue =',e15.7/)
c
      end 
      logical function chgsme (oldnrm,icnt)
c
c ... chgsme tests for jacobi si whether sme should be changed
c ... when caseii = .false..  if the test is positive the new value
c ... of sme is computed.
c
c ... parameter list --
c
c          oldnrm square of the norm of the pseudo-residual 
c                    at the last iteration
c          icnt   number of iterations since last change of 
c                    parameters
c
c *** begin -- itpackv common
c
      common /itcom1/ in,is,isym,itmax,level,nout,numwav
c
      logical adapt,betadt,caseii,halt,partad
      common /itcom2/ adapt,betadt,caseii,halt,partad
c
      common /itcom3/ bdelnm,betab,cme,delnnm,delsnm,ff,gamma,omega,qa,
     *   qt,rho,rrr,sige,sme,specr,spr,srelpr,stptst,udnm,zeta
c
c *** end   -- itpackv common
c
c     description of variables in common blocks in main r