c-------------------------------------------------------c
c     List of Routines:
c     
c     gmres
c     
c-------------------------------------------------------c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Non-symmetric solvers                               c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      
      subroutine gmres(x,d,m,nit,tol,ibnd,prod,zsol,nvec,tol2,iout,it,
     *     sum,sparse,irwste, indexe)
      implicit double precision (a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     preconditioned gmres.                                 c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++c
c     Ref. Saad, Y. & Schultz M.H., 1986.   GMRES :         c
c     A Generalised                                         c
c     Minimal Residual Algorithm for Solving Non-symmetric  c
c     Linear Systems. SIAM J. Sci. Stat. Comput. V7 #3 July c
c     1986 pp 856-869.                                      c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++c
c     This version uses Fast Householder instead of Givens  c
c     rotations to perform the QR factorisation in the      c
c     residual minimisation.                                c
c-------------------------------------------------------c
c     solution of the m-vector linear system Ax = d or      c
c     x is feasible solution on input/        c
c     solution vector on output             c
c     d is the right hand vector for the      c
c     system Ax = d                         c
c     A is a non-singular non-symmetric matrixc
c     is implicitly defined in the segments prod(h,z,m)     c
c     which compute the product of A with an m-vector       c
c     z to produce an m-vector h.                           c
c     this allows the user                                  c
c     to define their own (sparse) data structure.          c
c     the pre-conditioning matrix is defined by :           c
c     zsol(z,r,m) must be                                   c
c     written to solve the system  LRz = r.,                c
c     The preconditioned system is then                     c
c     (LR)**-1 A x = (LR)**-1 d                            c
c     where LR is an approximation to A                     c
c     m is the number of independent          c
c     variables.                            c
c     nit is the maximum number of            c
c     iterations allowed.                 c
c     tol is the stopping tolerance           c
c     ibnd is an integer array. normally set ibnd(i)=1      c
c     for i=1,2,3,..,m. this is useful for finite           c
c     element problems where A is singular, but some        c
c     of the variable x(i) are known. usually in such       c
c     problems one would "overwrite the boundary            c
c     conditions" and overwrite A. This overwriting can     c
c     be completely avoided if one follows the following    c
c     rules :                                               c
c     (1) enter the routine with each of the known         c
c     values x(i) in their correct positions in the         c
c     solution vector x.                                    c
c     (2) the integer array ibnd should have the values    c
c     ibnd(i)=1 if x(i) is one of the unknown values  c
c     ibnd(i)=0 if x(i) is a known value.             c
c     the routine will then correctly find the solution     c
c     using the matrix A without overwriting.               c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Author : N.R.C Birkett                                c
c     Oxford University Computing Laboratory       c
c     11 Keble Road                                c
c     Oxford OX1 3QD                               c
c     email  : nrcb@uk.ac.oxford.na                         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      include 'params.for'
      parameter  (mdim=411)
      dimension x(*),d(*),ibnd(*)
      dimension sparse(nmax*nzeromax),  indexe(nmax*nzeromax), 
     *     irwste(nmax+1)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Version December 1992 : Economised storage            c
c     Continuous Householder updates => residual calculated c
c     from H matrix without need to calculate  x at each    c
c     sweep.                                                c
c     Max vectors allowed for is GMRES(mdim-1)              c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   Big vectors
ccc
      dimension ik(nmax)
      common/cgr/qd(nmax,1:mdim),vk(nmax)
ccc
ccc   Small vectors
ccc
      dimension amat(mdim,mdim),bvec(mdim),vv(mdim)
      dimension hmat(mdim,mdim),hvec(mdim)
      external prod,zsol
ccc
ccc nvec < mdim
ccc
cccc
cccc  ioption = 1 gives convergence if r'r < tol**2 *d'd
cccc  AND r'r < tol2**2
cccc  ioption = 2 gives convergence if r'r < tol**2 *r0'r0
cccc                              AND r'r < tol2**2
cccc
cccc  iout    = 1 writes residual at each sweep to standard output
cccc  iout    = 0 suppresses printing.
cccc  
ccccc idump=1 to dump Upper Hessenberg matrix
ccccc in MATLAB format
      
ccc   for preconditioning, two simple versions
c      call MATRIX_DEFINE(m)
c      call MATRIX_DEFINE_DIAG(m)
      

c      print*,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(10)
c      print*,d(1),d(2),d(3),d(4),d(5),d(6),d(7),d(10)

      if(m.gt.nmax)stop'Array error in gmres. Too many variables'
      idump=0
      ioption=2
c      iout=1

      ihist=0
      if(ihist.eq.1) then
         open(22,file='history.cg',status='UNKNOWN')
      endif



cccc  h(1,1) contains the norm of the residual at the
cccc  start of nvec sweeps.
cccc  abs(h(i+1,i+1)) contains the norm of the residual
cccc  after i sweeps.
cccc  

      if(iout.eq.1) then
         write(*,400)nvec
      endif
 400  format('GMRES3(',I4,')')
      if(nvec.gt.mdim-1) then
         print *,'Array bound error in gmres3'
         print *,'re-set mdim >',nvec
         print *,'Solving with ',mdim-1,' vectors'
         nvec=mdim-1
      endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This algorithm is GMRES(nvec). Increasing nvec(< mdim) c
c will in general reduce the number of iterations to     c
c     convergence, but increase the execution time per step  c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

ccc
ccc   Initialisation
ccc   
      
      iflag=0
      it=0
      nk=0
ceee ...we do not use it
c      do 77 j=1,m
c         if(ibnd(j).eq.0) then
c            nk=nk+1
c            ik(nk)=j
c         endif
c 77   continue
      
ccc
ccc   Loop back to  78  instead of 500
ccc   if non-recursive residual required
ccc
 78   continue
      
      dnorm=0.0
ccc   
ccc   Get new rhs (LU)**-1 d
ccc   

      call vcopy(vk,d,m)
ceee ...we do not use it
c      call mask(vk,ik,nk)
      call zsol(qd(1,1),vk,m)
ceee ...we do not use it
c      call mask(qd(1,1),ik,nk)

c      print*,qd(1,1),qd(2,1),qd(3,1),qd(4,1),qd(5,1),qd(6,1),
c     *     qd(7,1),qd(10,1)

ccc
ccc   (LU)**-1 Ax (arbitrary preconditioner defined by zsol)
ccc

      call prod(vk,x,m,sparse,irwste, indexe)
ceee ...we do not use it
c      call mask(vk,ik,nk)
      call zsol(qd(1,2),vk,m)
ceee ...we do not use it
c      call mask(qd(1,2),ik,nk)

c      print*,qd(1,2),qd(2,2),qd(3,2),qd(4,2),qd(5,2),qd(6,2),
c     *     qd(7,2),qd(10,2)

c      do i=1,m
c         if(abs(qd(i,1) - qd(i,2) ).gt. 1E-06) then
c             print*,'@@@',i, abs(qd(i,1) - qd(i,2) ),
c     *           x(i),d(i),qd(i,1),qd(i,2)
c          endif
c       enddo

      
c	do i=1,5
c	   write(*,'(a2,2i5,4e12.4)')'$$',i,m,d(i),qd(i,1),vk(i),qd(i,2)
c	   write(*,'(a2,2i5,4e12.4)')'$$',i,m,d(i),x(i),vk(i),vk(i)-d(i)
c	enddo

c        do i=1,m
c	   write(*,'(a6,i5,4e14.6)')
c	1	'MATRIX',i,x(i),vk(i),qd(i,2)
c	enddo



ccc
ccc   The initial (preconditioned) residual is now in qd(*,1)
ccc
      do 14 j=1,m
         qd(j,1)=qd(j,1)-qd(j,2)
 14   continue
      
      call xdotx(dnorm,qd(1,1),m)

ccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccc
ccc   MAIN GMRES(NVEC) LOOP       ccc
ccc   Loop to 78 for non-recursiveccc
ccc   residual.                   ccc
ccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccc
 500  continue
c     write to history file
      if(ihist.eq.1) then
         call prod(vk,x,m,sparse,irwste, indexe)
         rr=0.0D0
         do 101 j=1,m
            rr=rr+(vk(j)-d(j))**2*ibnd(j)
 101     continue
         write(22,*)it,rr
 2000    format(i6,3(1x,1PD21.14))
      endif
      
      call zero(bvec,nvec+1)
      call xdotx(sum,qd(1,1),m)
      qtq=sqrt(sum)
      amat(1,1)=qtq
      hmat(1,1)=qtq
      bvec(1)=  qtq
      ss=1.0D0/qtq
      call vmult(ss,qd(1,1),m)
 1000 format('Pass...',i6,' residual**2 = ',1pd10.3)
      if(it.eq.0) then
         if(ioption.eq.1) sum0=dnorm
         if(ioption.eq.2) sum0=sum
         if(iout.eq.1) then
            write(*,1000)it,sum
            print*,'#####################################'
            write(*,'(i5,6e12.4)') 
     *           it,sum,sum0,tol,tol2,tol**2*sum0,tol2**2
         endif
      endif
      if(sum.le.tol**2*sum0.and.sum.le.tol2**2) then
         if(ihist.eq.1) then
            
            call prod(vk,x,m,sparse,irwste, indexe)
            rr=0.0D0
            do 102 j=1,m
               rr=rr+(vk(j)-d(j))**2*ibnd(j)
 102        continue
            write(22,*)it,rr
            close(22)
         endif
         return
      endif

      kvec=nvec
      do 20 i=1,nvec
ccc
ccc   w(i):=Aq(i)
ccc   
         call prod(vk,qd(1,i),m,sparse,irwste, indexe)
ceee ...we do not use it
c         call mask(vk,ik,nk)
         call zsol(qd(1,i+1),vk,m)
ceee ...we do not use it
c         call mask(qd(1,i+1),ik,nk)
ccc     
ccc Gram-Schmidt Orthonormalise to get {q(1),q(2),..,q(nvec+1)}
ccc (It would be better to use modified  Gramm-Schmidt
ccc  if nvec large)
         do 30 j=1,i
            call xdoty(amat(j,i+1),qd(1,j),qd(1,i+1),m)
            call aminvb(qd(1,i+1),qd(1,j),amat(j,i+1),m)
 30      continue

ccc
ccc Normalise q(i+1)
ccc
         call xdotx(qtq,qd(1,i+1),m)
         qtq=sqrt(qtq)
         ss=1.0D0/qtq
         amat(i+1,i+1)=qtq
         call vmult(ss,qd(1,i+1),m)
ccc   
ccc Copy next column of H as Householder overwrites.
ccc
         call vcopy(hmat(1,i+1),amat(1,i+1),i+1)
         call house(amat(1,2),vv,bvec,i-1,mdim)
         it=it+1
         
ccc     
ccc     The 2-norm of the residual is available
ccc as abs(bvec(i+1) at each sweep without
ccc     needing to calculate r:=r-Av
ccc     
ccc

         sum=bvec(i+1)**2
         if(iout.eq.1) then
c            write(*,1000)it,sum
            write(*,'(i5,6e12.4)') 
     *           it,sum,sum0,tol,tol2,tol**2*sum0,tol2**2
         endif
         
         if(sum.le.tol**2*sum0.and.sum.le.tol2**2) then
            iflag=1
            kvec=i
            goto 5000
         endif
 20   continue

ccc
ccc   At the end of nvec cycles update r,x by completing
ccc   the QR calculation.
ccc   Solve the least squares problem
ccc   min || rk - AQr0||2
ccc   which is the same as
ccc   min || (amat(1,1),0,0,0..0)T- amat r0 ||
ccc   where amat is nvec+1 X nvec upper Hessenberg
ccc   Use fast version of Householder for this
ccc
 5000 continue
      if(idump.eq.1) then
         call hdump(hmat,kvec,mdim)
         idump=0
      endif
      call hsolve(amat(1,2),bvec,kvec+1,mdim)

ccc
ccc   x(k+1)=x(k)+v(k)
ccc   r(k+1)=r(k)-w(k)
ccc   where v(k) = Q r0
ccc   w(k) = W r0
ccc

 3333 format(A2,20(1x,f10.3))

ccc
ccc   compute hvec = H'bvec, then W bvec = QH'bvec = Q hvec
ccc
      call zero(hvec,kvec+1)
      do 6020 j=2,kvec+1
         bi=-bvec(j-1)
         call aminvb(hvec,hmat(1,j),bi,j)
 6020 continue
ccc
ccc   Restore previous residual to q(*,1)
ccc   by mutiplying by mod(r) so do i=1
ccc   as special case.
ccc
      i=1
      bi=bvec(i)
      hi=hmat(1,1)-hvec(1)
      call aminvb(x,qd(1,i),-bi,m)
      call vmult(hi,qd(1,1),m)
ccc   
ccc   Do the rest
ccc
      do 8030 i=2,kvec
         bi=bvec(i)
         hi=hvec(i)
         call aminvb(x,qd(1,i),-bi,m)
         call aminvb(qd(1,1),qd(1,i), hi,m)
 8030 continue
      call aminvb(qd(1,1),qd(1,kvec+1),hvec(kvec+1),m)
      
      if(it.ge.nit) goto 9000
ccc
ccc   Loop back to 500
ccc   to use recursive residual
ccc   otherwise loop back to 78
ccc   (exact residual requires an
ccc   extra mat/vec multiply.
      if (iflag.eq.0) then
         goto 78
      else
         goto 500
      endif
 9000 continue
c      print *,'Gmres  failed to converge in ',it,' iterations'
      if(ihist.eq.1) then
         close(22)
      endif
      return
      end

      subroutine hdump(hmat,m,mdim)
      implicit double precision (a-h,o-z)
      dimension hmat(mdim,*)
ccc
ccc   Dump the upper hessenberg matrix in
ccc   matlab format
ccc
      open(37,file='hmat',status='UNKNOWN')
      do 10 i=1,m
         write(37,1000)(hmat(i,j),j=2,m+1)
 1000    format(1000(1PE21.14,1x))
 10   continue
      close(37)
      return
      end


      subroutine mask(x,ik,nk)
      implicit double precision (a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccc
c     Zero displacements of known values    c
ccccccccccccccccccccccccccccccccccccccccccccc
      dimension x(*),ik(*)
      do 10 j=1,nk
         x(ik(j))=0.0D0
 10   continue
      return
      end
      

      subroutine xdoty(s,x,y,n)
ccccccccccccccccccccccccccccccccccccccccccccc
c     s = <x,y> real inner product of 2 vectors c
ccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      dimension x(*),y(*)
      s=0.0d0
      do 10 j=1,n
         s=x(j)*y(j)+s
 10   continue
      return
      end


      subroutine xdotx(s,x,n)
ccccccccccccccccccccccccccccccccccccccccccccc
c     s = <x,x> real inner product of  x with x c
ccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      dimension x(*)
      s=0.0d0
      do 10 j=1,n
         s=s+x(j)*x(j)
c         if(x(j)**2 .gt. 1E-04) 
c     *        print*,'X.X',j,n,x(j),x(j)**2,s
 10   continue
      return
      end
      
      

      subroutine vcopy(a,b,n)
      implicit double precision (a-h,o-z)
      dimension a(*),b(*)
ccccccccccccccccccccccccccccccccccccccccccccccc
c     Make a copy of vector b in vector a         c
ccccccccccccccccccccccccccccccccccccccccccccccc
      do 10 j=1,n
         a(j)=b(j)
 10   continue
      return
      end




      subroutine vmult(v,x,n)
      implicit double precision (a-h,o-z)
cccccccccccccccccccccccccccccccccccccccccccccccc
c     Vector X : = v*X, v scalar                   c
c     len(x)=n                                     c
cccccccccccccccccccccccccccccccccccccccccccccccc
      dimension x(*)
      do 10 j=1,n
         x(j)=v*x(j)
 10   continue
      return
      end
      
      subroutine aminvb(a,b,v,n)
      implicit double precision (a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccc
c     Vector A : = A - v*B                        c
c     v scalar, len(a)=len(b)=n                   c
ccccccccccccccccccccccccccccccccccccccccccccccc
      dimension a(*),b(*)
      J1X=N-(N/8)*8
      do 05 j=1,J1X
         a(j)=a(j)-v*b(j)
 05   continue
      do 10 j=J1X+1,n,8
         a(j)=a(j)-v*b(j)
         a(j+1)=a(j+1)-v*b(j+1)
         a(j+2)=a(j+2)-v*b(j+2)
         a(j+3)=a(j+3)-v*b(j+3)
         a(j+4)=a(j+4)-v*b(j+4)
         a(j+5)=a(j+5)-v*b(j+5)
         a(j+6)=a(j+6)-v*b(j+6)
         a(j+7)=a(j+7)-v*b(j+7)
 10   continue
      return
      end
      
      subroutine asembl(a,sparsu,indexu,irwstu,i,j)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     assemble  scalar a(i,j) in upper sparse array, j>=i c
c     a(j,i) in lower sparse array,      c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision a,sparsu(*)
      integer indexu(*),irwstu(*)
c     
c     Locate base address  and end address
c     
      istart=irwstu(i)
      iend=irwstu(i+1)-1
c     
c     check if element is available and store
c     
      do 10 k=istart,iend
         if(indexu(k).eq.j) then
            sparsu(k)=a+sparsu(k)
            return
         endif
 10   continue
      
c
c     succesful exit from loop => error in indices i,j
c     
      print *,' Index error in asembl for i,j=',i,j
      stop
      end
      
      subroutine ilu(b,sparsu,indexu,irwstu,
     $     sparsl,indexl,irwstl,n)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Overwite the incomplete lower and upper factors              c
c     on sparsl,sparsu using same indexing as A                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      dimension sparsu(*),sparsl(*)
      integer indexu(*),indexl(*),irwstu(*),irwstl(*)
      dimension b(*)
c     b is a work vector of length n
c     required for scatter and gather
c     No pivotting  !
      do 10 k=1,n-1
c     
c     Operate on U first where A=L+U
c     locate start of column k
c     
         istart=irwstl(k)
         iend=irwstl(k+1)-1
         kk=irwstu(k)
c     
c     loop down column k
c
         aa=1.0D0/sparsu(kk)
         sparsu(kk)=aa
         do 21 i=istart+1,iend
            sparsl(i)=sparsl(i)*aa
 21      continue
         do 20 i=istart+1,iend
            v=sparsl(i)
            index=indexl(i)
            jstart=irwstu(index)
            jend=irwstu(index+1)-1
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Subtract v*row(k) from each row(i)                c
c     This is the cpu intensive part                    c
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     scatter
            do 105 j=jstart,jend
               b(indexu(j))=sparsu(j)
 105        continue
c       saxpy
            do 110 j=irwstu(k)+1,irwstu(k+1)-1
               b(indexu(j))=b(indexu(j))-v*sparsu(j)
 110        continue
c     gather
            do 120 j=jstart,jend
               sparsu(j)=b(indexu(j))
 120        continue
            
 20      continue
         jend=irwstu(k+1)-1
         do 40 j=kk+1,jend
            v=sparsu(j)
            index=indexu(j)
            istart=irwstl(index)
            iend=irwstl(index+1)-1
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Subtract v*column(k) from each column j           c
c     This is the cpu intensive part                    c
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     scatter
            do 205 i=istart,iend
               b(indexl(i))=sparsl(i)
 205        continue
c     saxpy
            do 210 i=irwstl(k)+1,irwstl(k+1)-1
               b(indexl(i))=b(indexl(i))-v*sparsl(i)
 210        continue
c     gather
            do 220 i=istart,iend
               sparsl(i)=b(indexl(i))
 220        continue
            
 40      continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     17 July 1991:change factorisation to iLDU by dividing  c
c     through by the diagonal in sparsu                      c
c     5 August 1992: Store 1/diag in sparsu. This affects    c
c     routine sinv                                           c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         istart=irwstu(k)
         iend=irwstu(k+1)-1
         diag=sparsu(istart)
         do 80 j=istart+1,iend
            sparsu(j)=sparsu(j)*diag
 80      continue
 10   continue
      sparsu(irwstu(n))=1.0D0/sparsu(irwstu(n))
      return
      end
      
      subroutine sort(index,k2)
cccccccccccccccccccccccccccccccccccccccc
c     Sort the index vector into numerical c
c     order . This is not the best sort    c
c     routine , but is useful for short    c
c     initialisations                      c
cccccccccccccccccccccccccccccccccccccccc
      dimension index(*)
      do 10 j=1,k2
         do 20 i=j+1,k2
            m1=index(j)
            m2=index(i)
            if(m1.gt.m2) then
               index(j)=m2
               index(i)=m1
            endif
 20      continue
 10   continue
      return
      end

      subroutine sortin(index,irwst,n)
ccccccccccccccccccccccccccccccccccccccccc
c     Routine to sort the row index blocks  c
c     into numerical order. index(k)=k and  c
c     by definition index(k+j)>k for        c
c     j=0,1,...,irwst(k+1)-irwst(k)         c
c     There are n rows in the matrix        c
ccccccccccccccccccccccccccccccccccccccccc
      dimension index(*),irwst(*)
      do 10 k=1,n
         call sort(index(irwst(k)+1),irwst(k+1)-irwst(k)-1)
 10   continue
      return
      end
      
      subroutine prod(c,x,n,sparse,irwste, indexe)
      implicit double precision (a-h,o-z)
C     C Mat-Vec Mult Using sparse arrays
      include 'params.for'
c      common/mats/sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
c     $     spu(nmax*nzeromax),spl(nmax*nzeromax),
c     $     index(nmax*nzeromax),irwst(nmax)
      common/matsC/ spars(nmax*nzeromax),index(nmax*nzeromax),
     $     irwst(nmax+1)

      dimension c(*),x(*)
      dimension sparse(nmax*nzeromax),  indexe(nmax*nzeromax), 
     *     irwste(nmax+1)
c      call prodax(x,sparsu,index,irwst,
c     $     sparsl,index,irwst,c,n)

c     ... new version of storing of matices
c      do i=1,n
c         print*,'###',irwst(i),irwste(i),irwst(i)-irwste(i)
c      enddo

      do 10 i=1,n
         istart=irwste(i)
         iend=irwste(i+1)-1
         sum=0.0D0
         do 20 j=istart,iend
            sum=sum+x(indexe(j))*sparse(j)
cc            sum=sum+x(index(j))*spars(j)
 20      continue
         c(i)=sum
 10   continue


      return
      end
      
      
      subroutine ilufac(z,g,n)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     ILDU preconditioner                                   c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      include 'params.for'
c      common/mats/sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
c     $     spu(nmax*nzeromax),spl(nmax*nzeromax),
c     $     index(nmax*nzeromax),irwst(nmax)
      common/matsC/ spars(nmax*nzeromax),index(nmax*nzeromax),
     $     irwst(nmax+1)
      dimension z(*),g(*)


      do 10 j=1,n
         z(j)=g(j)
 10   continue
      
c      call sfwd(z,spl,index,irwst,n)
c      call sinv(z,spu,index,irwst,n)
c      call sbkd(z,spu,index,irwst,n)
      return
      end
      
      SUBROUTINE HOUSE(A,V,B,L,MDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MDIM,*),B(*),V(*)
c======================================================================c
c     Version for m X m-1 upper Hessenberg matrix                          c
c     For solution of the least-squares problem min || Ax-b||  A is M X N  c
c     M> = N, first use the routine QR to factor A. Then do the            c
c     2 steps b:= Q'b, solve for x the equations Rx = b, and overwrites thec
c     solution on b. This routine may be called any number of times        c
c     Re-ordered routine written for GMRES algorithm :                     c
c     (1) Apply Householder transformations to next column generated       c
c     from Arnoldi ,                                                   c
c     A(*,L+1):=P(L)P(L-1)..P(1)A(*,L+1)                                   c
c     Calculate new P(L+1) to eliminate A(L+1,L)                           c
c     and apply P(L+1) to vector b                                         c
c     V is a work vector.                                                  c
c======================================================================c
      DO 10 K=1,L
         HK=-V(K)*0.5D0/A(K,K)
         call xdoty(VB,A(K+1,K),a(k+1,L+1),1)
         vb=vb+hk*A(k,L+1)
         D=2.0D0*VB/V(K)
         call aminvb(a(k+1,L+1),A(k+1,k),d,1)
         a(k,L+1)=a(k,L+1)-d*hk
 10   CONTINUE
      K=L+1
      call xdotx(SS,A(k,k),2)
      S=SQRT(SS)
      IF (A(K,K).LT. 0.0) S=-S
c     Hkk = Vk/(-2rkk)
C     STORE VTV IN V
      V(K)=2.0D0*SS+2.0D0*A(K,K)*S
C     APPLY TRANSFORMATION,GET NEW A
      A(K,K)=-S
c     and apply to get new right handside
      HK=-V(K)*0.5D0/A(K,K)
      call xdoty(VB,A(K+1,K),b(k+1),1)
      vb=vb+hk*b(k)
      D=2.0D0*VB/V(K)
      call aminvb(b(k+1),A(k+1,k),d,1)
      b(k)=b(k)-d*hk
      RETURN
      END
      

      SUBROUTINE HSOLVE(A,B,M,MDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MDIM,*),B(*)
      N=M-1
c=================================================================c
c     Now solve Rx = Q'b and overwrite the solution on b
c=================================================================c
      do 40 k=n,1,-1
         b(k)=b(k)/a(k,k)
         call aminvb(b(1),a(1,k),b(k),k-1)
 40   continue
      RETURN
      END
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Split Ly=b,Dy=z,Ux=z into 3 separate routines.           c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine sfwd(b,sparsl,indexl,irwstl,n)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     sparse forward substitution                              c
c     the sparse ilu factors are contained in                  c
c     sparsl,sparsu using the sparse data-structure            c
c     on entry b contains the n-dimensional right hand vector  c
c     on exit b contains (iu)**-1(il)**-1 b                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      dimension sparsl(*),b(*)
      integer indexl(*),irwstl(*)
c
c sparse forwards using lower factor
c
      do 10 j=1,n-1
         istart=irwstl(j)
         iend=irwstl(j+1)-1
         do 20 i=istart+1,iend
            b(indexl(i))=b(indexl(i))-sparsl(i)*b(j)
 20      continue
 10   continue
      return
      end
      
      subroutine sinv(b,sparsu,indexu,irwstu,n)
      implicit double precision (a-h,o-z)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Sparse diagonal inversion Dz=y                         c
c     using upper strorage in sparsu                         c
c     05 August 1992. 1/d stored in diagonal to make         c
c     this a multiply rather than divide.                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension sparsu(*),b(*)
      integer indexu(*),irwstu(*)
c       
c       divide by diagonal factors
c
      do 10 k=1,n
         b(k)=b(k)*sparsu(irwstu(k))
 10   continue
      return
      end
      


      subroutine sbkd(b,sparsu,indexu,irwstu,n)
      implicit double precision (a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Back substitute using upper factor                   c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension sparsu(*),b(*)
      integer indexu(*),irwstu(*)
      do 30 i=n-1,1,-1
         jstart=irwstu(i)
         jend=irwstu(i+1)-1
         do 40 j=jstart+1,jend
            b(i)=b(i)-sparsu(j)*b(indexu(j))
 40      continue
 30   continue
      return
      end

      subroutine prodax(x,sparsu,indexu,irwstu,
     $     sparsl,indexl,irwstl,prod,n)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     unsymmetric Ax. To get A'x simply reverse roles of upper andc
c     lower data                                                  c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Use data structures to form prod= Ax =(L+U)x no clearing of c
c     prod required.                                              c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      dimension x(*),sparsl(*),sparsu(*),prod(*)
      integer   indexu(*),irwstu(*)
      integer   indexl(*),irwstl(*)
      call upperx(x,sparsu,indexu,irwstu,prod,n)
      call lowerx(x,sparsl,indexl,irwstl,prod,n)
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     mat-vec multiplies using the data structure    cccccc
c     Break this down into Ax = Lx + Ux where        cccccc
c     L is the strict lower triangle of A and U is   cccccc
c     the upper triangle including the diagonal.     cccccc
c     Similarly A'x = L'x + U'x where L' is upper    cccccc
c     and U' is lower triangular. To vectorise need  cccccc
c     separate rotines for upper *x and lower*x      cccccc
c     Lx = sum (i=1,n) Ci xi                         cccccc
c     i.e scalar multiply the columns of L with the  cccccc
c     components of x and add resultant vectors.     cccccc
c     whereas                                        cccccc
c     Ux = (r1.x,r2.x,...,rn.x)'                     cccccc
c     i.e take inner products of x with rows of U    cccccc
c     to get the components of Ux.                   cccccc
c     finally Ax=Lx+Ux - add the 2 partial sums.     cccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc


      subroutine lowerx(x,sparsl,indexl,irwstl,prodl,n)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     form sparse Lx from datastructure and store in prodl  c
c     A is n X n                                            c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision sparsl(*),prodl(*),x(*)
      integer indexl(*),irwstl(*)
      double precision xx
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     loop on columns of L                                  c
c     If Lx only is required then clear prodl before entry  c
c     the n'th column has no entries below the diagonal     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      do 20 j=1,n-1
c     
c     get base address of the diagonal
c     leave diagonal multiply to Ux routine
         istart=irwstl(j)
         iend= irwstl(j+1)-1
         xx=x(j)
         do  30 i=istart+1,iend
            prodl(indexl(i))=prodl(indexl(i))+xx*sparsl(i)
 30      continue
 20   continue
      return
      end
      
      subroutine upperx(x,sparsu,indexu,irwstu,produ,n)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     form sparse Ux from datastructure and store in prodl  c
c     A is n X n                                            c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision sparsu(*),produ(*),x(*)
      integer indexu(*),irwstu(*)
      double precision sum
c       
c     loop on rows of U
c     form inner products
c     
      do 10 i=1,n
         istart=irwstu(i)
         iend=irwstu(i+1)-1
         sum=0.0D0
         do 20 j=istart,iend
            sum=sum+x(indexu(j))*sparsu(j)
 20      continue
         produ(i)=sum
 10   continue
      return
      end
      
      subroutine zero(x,n)
      implicit double precision (a-h,o-z)
      dimension x(*)
      do 10 j=1,n
         x(j)=0.0D0
 10   continue
      return
      end
      
      subroutine prod1(c,x,n)
      implicit double precision (a-h,o-z)
C     C Mat-Vec Mult Using sparse arrays
      include 'params.for'
      common/matsC/ spars(nmax*nzeromax),index(nmax*nzeromax),
     $     irwst(nmax+1)
c      common/mats/sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
c     $     spu(nmax*nzeromax),spl(nmax*nzeromax),
c     $     index(nmax*nzeromax),irwst(nmax)
      dimension c(*),x(*)
c     call prodax(x,sparsu,index,irwst,
c     1    sparsl,index,irwst,c,n)
      
      do i=1,n
         c(i) = 2*x(i)
         if(i .gt. 1) c(i) =c(i) + x(i-1)
         if(i .lt. n) c(i) =c(i) + x(i+1)
      enddo
      return
      end
      
      subroutine MATRIX_DEFINE(m)
      implicit double precision (a-h,o-z)
      
      include 'params.for'
      common/mats/ sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
     $     spu(nmax*nzeromax),spl(nmax*nzeromax),
     $     index(nmax*nzeromax),irwst(nmax)
      
      do i=1,m
         istart=irwst(i)
         iend=irwst(i+1)-1
         do j=istart,iend
            spu(j) = sparsu(j)
            spl(j) = sparsl(j)
         enddo
      enddo


      return
      end


      subroutine MATRIX_DEFINE_DIAG(m)
      implicit double precision (a-h,o-z)
      include 'params.for'
      common/mats/ sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
     $     spu(nmax*nzeromax),spl(nmax*nzeromax),
     $     index(nmax*nzeromax),irwst(nmax)
      
c
c     diagonal preconditioners
c

      do i=1,m
         istart=irwst(i)
         iend=irwst(i+1)-1
         do j=istart,iend
            spu(j) = 0
            spl(j) = 0
         enddo
      enddo
      do i=1,m
         istart=irwst(i)
         iend=irwst(i+1)-1
         do j=istart,iend
            spu(istart) =  spu(istart) + sparsu(j)**2
            spu(index(j)) =  spu(index(j)) + sparsl(j)**2
         enddo
      enddo

      do i=1,m
         spu(irwst(i)) = spu(irwst(i))**0.5
      enddo

      return
      end



c     --------------------------------------------------------------
c     This subroutine sets up the preconditioner
c     --------------------------------------------------------------
      subroutine setup_precon(nnodes,bcnods)
c     --------------------------------------------------------------
c     implicit none
      implicit double precision (a-h,o-z)
      include 'params.for'
      
c     double precision  b(nodegmax),sparsu,sparsl,spu,spl
c     integer index,irwst,j,nnodes,i,
      integer bcnods(nnodes)
c     


      common/mats/
     &     sparsu(nmax*nzeromax),sparsl(nmax*nzeromax),
     &     spu(nmax*nzeromax),spl(nmax*nzeromax),
     &     index(nmax*nzeromax),irwst(nmax)


c
c
c
c       Copy array to spu, overwrite and
c       copy to spu & spl then IDLU factorise
c
      do 230 j=1,irwst(nnodes)
         spu(j)=sparsu(j)
         spl(j)=sparsl(j)
 230  continue
c
c       Overwrite boundary rows/columns with unit submatrix
c
        do 240 i=1,nnodes
          do 242 j=irwst(i),irwst(i+1)-1
            if(bcnods(index(j)).eq.0) then
              spu(j)=0.0d0
              spl(j)=0.0d0
            endif
242       continue
          if(bcnods(i).eq.0) then
c
c       Row i & column i
c
            do 250 j=irwst(i),irwst(i+1)-1
              spu(j)=0.0d0
              spl(j)=0.0d0
250         continue
c
c       Unit on the diagonal
c
            spu(irwst(i))=1.0d0
            spl(irwst(i))=1.0d0
          endif
240     continue
c
c        call ilu(b,spu,index,irwst,spl,index,irwst,nnodes)
c
        return
        end

