implicit real(a-h,o-z)
      real a(200000),af(200000),b(28000),x(28000),
     $      anrm,xnrm,resd,g(28000),p(28000),t1(28000),t2(28000)
      integer ia(28001),ja(200000),nit(3)
      real d,c2,a2,cf
      common /tomm/ a,af,ia,ja
      common /area2/ d,c2,a2,cf
       common /order/ iorflg
       data nit/1,5,10/
c
c     this is a main program for running the incomplete LU
c     conjugate gradient for solving linear systems.
c     the matrix need not be symmetric.
c
c     ilucg this flag if 1 will cause incomplete factoization
c               if 0 no factorization just straight cg
c
c     write(6,*)' start in main'
c
c     call trapov(100)
      loca = 200000
    1 continue
c     write(6,*) 'input ineum5,ineum6,irot,iorflg'
c     write(6,*)'set irot<0 to stop'
c     read(5,*) ineum5,ineum6,irot,iorflg
c     if( irot .lt. 0 ) then
c       write(6,*)' end test'
c       stop
c     end if
c     write(6,*) 'ineum5,ineum6,irot,iorflg',ineum5,ineum6,irot,iorflg
c     write(6,*)'input mesh cells for x, y, and z'
c     read(5,*)imx,imy,imz
c     write(6,*)'inner, middle, and outer fill-in integers'
c     read(5,*)ifli,iflm,iflo
cc    write(6,*)'ising (normally zero)'
cc    read(5,*)ising
      ineum5 = 1
      ineum6 = 1
      irot = 1
      iorflg = 2
      imx = 20
      imy = 20
      imz = 20
cc    if (ineum5 .ne. 1 .or. ineum6 .ne. 1) ising=0
      ising = 0
cc    imx = 5
cc    imy = 5
cc    imz = 5
      n = imx*imy*imz
c     write(6,*)   ' *********** problem type *************'
c     if( ineum5 .eq. 0 ) then
c        write(6,*)' *****       dir on bottom        *****'
c     else
c        write(6,*)' *****     neuman on bottom       *****' 
c     endif
c     if( ineum6 .eq. 0 ) then
c        write(6,*)' *****         dir on top         *****'
c     else
c        write(6,*)' *****      neuman on top         *****'
c     endif
c     if( ising .eq. 0 ) then
c        write(6,*)' *****         regular            *****'
c     else
c        write(6,*)' ***** singular if neuman on top and bottom *****'
c     endif
c     write(6,*)   ' ***** mesh cells in x dir',imx,'     *****'
c     write(6,*)   ' ***** mesh cells in y dir',imy,'     *****'
c     write(6,*)   ' ***** mesh cells in z dir',imz,'     *****'
c
c     ifli is the number of fillin stripes allowed for each of the
c          inner stripes; if zero no fillin.
c     iflm is the number of fillin stripes allowed for each of the
c          middle stripes; if zero no fillin.
c     iflo is the number of fillin stripes allowed for each of the
c          outer stripes; if zero no fillin here.
c     inum5 if 1 then neumann on bottom
c           if 0 then dirchlet on bottom 
c     inum6 if 1 then neumann on top
c           if 0 then dirchlet on top 
c     irot  if 0 then no rot flow field
c           if 1 then rot flow field
c     ising if zero regular
c           if 1 then signular if neumann top and bottom
c     imx # of mesh cells in the x direction
c     imy # of mesh cells in the y direction
c     imz # of mesh cells in the z direction
cc    ifli = 0
cc    iflm = 0
cc    iflo = 0
cc    write(6,*)' call prbtyp'
      call prbtyp(ineum5,ineum6,ising,imx,imy,imz,0)
cc    write(6,*)' call matgen'
      ti1 = second(ii)
      call matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,b)
c     write(6,*)'n,nonz,ia(n+1)'
c     write(6,*)n,nonz,ia(n+1)
cc    write(6,*)' after matgen'
      ti2 = second(ii) - ti1
      if( ia(n+1) .gt. loca ) then
         write(6,*)' not enough storage',loca,ia(n+1)
         stop
      endif
c      stop
      job = 0
      anrm = abs(a(isamax(ia(n+1)-1,a,1)))
      call scopy(n,b,1,t2,1)
      maxits = n
      tol = 1.0e-6
      write(6,101) ti2
  101 format(/20x,'time to generate matrix = ',5x,1pe10.3)
      do 888 ics = 1,3
      ti3 = second(ii)
      do 777 mm = 1, nit(ics)
      call scopy(n,t2,1,x,1)
      call scopy(n,t2,1,b,1)
      call iccglu(a,n,ia,ja,af,b,x,g,p,t1,tol,maxits,its,job,info)
c     write(6,*)'info should be zero  info = ',info
      job = 1
  777 continue
      ti4 = second(ii) - ti3
      call scopy(n,t2,1,b,1)
      call cgres(a,n,ia,ja,x,b,p)
      resd = abs(p(isamax(n,p,1)))
      xnrm = abs(x(isamax(n,x,1)))
c     write(6,*)'anrm,xnrm,resd',anrm,xnrm,resd
      sres = resd/(anrm*xnrm)
      nonzr = ia(n+1) - 1
c     ti5 = ti2 + ti4
c     write(6,100) sres,ti5,ti2,ti4,nonzr,its
c 100 format(20x,'scaled residual = ',13x,1pe10.3,
c    1      /20x,'total time = ',18x,1pe10.3,
c    2      /20x,'time to generate matrix = ',5x,1pe10.3,
c    3      /20x,'time for cg iterations = ',6x,1pe10.3,
c    4      /20x,'number of nonzero elements = ',2x,i10,
c    5      /20x,'number of cg iterations = ',5x,i10)
      write(6,102) nit(ics),its,sres,ti4
  102 format(20x,'results for ',i3,' sets of ',i4,' cg iterations ',
     1      'for each set',
     2      //20x,'scaled residual = ',13x,1pe10.3,
     3      /20x,'time for cg iterations = ',6x,1pe10.3,//)
  888 continue
c     go to 1
      stop
      end
       subroutine matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,f)
      implicit real(a-h,o-z)
      integer ia(1),ja(1),ifli,iflo
      real a(1)
      integer bwi,bwo
      real f(28000),acof0(28000),acof1(28000),acof2(28000),
     # acof3(28000),
     # acof4(28000),acof5(28000),acof6(28000),gl(1000),gu(1000)
       real cvx(30,30,30),cvy(30,30,30),cvz(30,30,30)
     #  ,omg(30,30,30)
       common/coef/acof0,acof1,acof2,acof3,acof4,acof5,acof6
       common /order/ iorflg
c      ineum5=1 neuman on bot;=0 dir on bot.
c      ineum6=1 neuman on top;=0 dir on top.
c      irot = 1, rotational flow field
c      irot = 0, non-rotational flow field
      flx=1.e0
      fly=1.e0
      flz=1.e0
c      write(6,*)' call prbtyp'
       call prbtyp(ineum5,ineum6,ising,imx,jmx,kmx,1)
c      write(6,*)' out of prbtyp in matgen'
      id=1
      idp1=id+1
       iptflg=0
       go to (3,5,7),iorflg
 3     continue
c      write(6,*)' continue 3'
       bwo=imx*jmx
       bwi=imx
       ni=1
       nj=bwi
       nk=bwo
       nc=-bwo-bwi
       nib=1
       njb=imx
       ncb=-imx
c       (ijk) (5310246)
       go to 9
 5     continue
c      write(6,*)' continue 5'
c       (kij) (3150624)
       bwo=kmx*imx
       bwi=kmx
       ni=bwi
       nj=bwo
       nk=1
       nc=-bwo-bwi
       nib=1
       njb=imx
       ncb=-imx
       go to 9
 7     continue
c      write(6,*)' continue 7'
c       (jki) (1530462)
       bwo=jmx*kmx
       bwi=jmx
       ni=bwo
       nj=1
       nk=bwi
       nc=-bwo-bwi
       nib=jmx
       njb=1
       ncb=-jmx
 9     continue
c      write(6,*)' continue 9'
c      write(6,*)'imx,jmx,kmx',imx,jmx,kmx
      dx=flx/float(imx)
      dy=fly/float(jmx)
      dz=flz/float(kmx)
      rdx2=1.e0/dx**2
      rdy2=1.e0/dy**2
      rdz2=1.e0/dz**2
      mx=imx*jmx*kmx
      imxm1=imx-1
      jmxm1=jmx-1
      kmxm1=kmx-1
      nonz=3*mx-2+2*(mx-bwi+ifli)+2*(mx-bwo+iflo)
      nonzt=nonz
c      write(6,*)' call inits'
      call inits(a,mx,ia,ja)
      if(ifli .ge. bwi-1) go to 220
      if(iflo .ge. bwo-bwi) go to 220
      do 10 m=1,mx
      f(m)=0.e0
  10  continue
c      write(6,*)' continue 10'
      do 11 i=1,imx
      do 11 j=1,jmx
      do 11 k=1,kmx
      xi=(float(i)-0.5e0)*dx
      yj=(float(j)-0.5e0)*dy
      zk=(float(k)-0.5e0)*dz
      facx = 1.0e0
      facy = 1.0e0
      if( irot .eq. 0 ) go to 12
      facx = xi - .05e0
      facy = yj - .05e0
   12 continue
      dum=32.0e0*xi*(1.e0-xi)*yj*(1.e0-yj)*zk
      dum1=2.0e0*xi*yj*zk**2
       dum=25.e0*dum
       dum1=2.0e0*dum1
       cvx(i,j,k)=dum*facx
       cvy(i,j,k)=dum*facy
       cvz(i,j,k)=dum1
       omg(i,j,k)=0.e0
       dum2 = (xi*xi)*yj*zk
       m = i*ni + j*nj + k*nk + nc
       f(m) = dum2
11    continue
c      write(6,*)' continue 11'
      do 15 j=1,jmx
      do 15 i=1,imx
      mb=i*nib+j*njb+ncb
      gl(mb)=1.e0
      gu(mb)=2.0e0
  15  continue
c      write(6,*)' continue 15'
      do 20 m=1,mx
      acof0(m)=0.e0
      acof1(m)=0.e0
      acof2(m)=0.e0
      acof3(m)=0.e0
      acof4(m)=0.e0
      acof5(m)=0.e0
      acof6(m)=0.e0
  20  continue
c      write(6,*)' continue 20'
      if(imx .lt. 3) go to 45
      if(jmx .lt. 3) go to 45
      if(kmx .lt. 3) go to 45
      do 30 j=2,jmxm1
      do 30 i=2,imxm1
      do 30 k=2,kmxm1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=2.0e0*(rdx2+rdy2+rdz2)+omg(i,j,k)
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
  30  continue
c      write(6,*)' continue 30'
  45  continue
c      write(6,*)' continue 45'
      if((jmx.lt.3).or.(kmx.lt.3)) go to 55
      do 50 j=2,jmxm1
      do 50 k=2,kmxm1
      i=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+2.0e0*(rdy2+rdz2)+omg(i,j,k)-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      i=imx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+2.0e0*(rdy2+rdz2)+omg(i,j,k) +0.5e0*cvx(i,j,k)/dx
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
  50  continue
c      write(6,*)' continue 50'
  55  continue
c      write(6,*)' continue 55'
      if((imx.lt.3).or.(kmx.lt.3)) go to 65
      do 60 i=2,imxm1
      do 60 k=2,kmxm1
      j=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdy2+ 2.0e0*(rdx2+rdz2)+omg(i,j,k) -0.5e0*cvy(i,j,k)/dy
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      j=jmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdy2+2.0e0*(rdx2+rdz2)+omg(i,j,k)+0.5e0*cvy(i,j,k)/dy
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
  60  continue
c      write(6,*)' continue 60'
  65  continue
c      write(6,*)' continue 65'
      if((imx.lt.3).or.(jmx.lt.3)) go to 75
      do 70 i=2,imxm1
      do 70 j=2,jmxm1
      k=1
      m=(j-1)*bwo+(i-1)*bwi +k
      acof0(m)=2.0e0*(rdx2+rdy2)+3.0e0*rdz2 +omg(i,j,k)
     #  +cvz(i,j,k)/(2.0e0*dz)
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 71
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
   71 continue
c      write(6,*)' continue 71'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=2.0e0*(rdx2+rdy2)+3.0e0*rdz2+omg(i,j,k)
     #  -cvz(i,j,k)/(2.0e0*dz)
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 72
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
   72 continue
c      write(6,*)' continue 72'
  70  continue
c      write(6,*)' continue 70'
  75  continue
c      write(6,*)' continue 75'
      if(kmx.lt.3) go to 85
      do 80 k=2,kmxm1
      i=1
      j=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      j=jmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      i=imx
      m=(j-1)*bwo+(i-1)*bwi +k
      acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+ 0.5e0*cvz(i,j,k)/dz
      j=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
  80  continue
c      write(6,*)' continue 80'
  85  continue
c      write(6,*)' continue 85'
      if(imx.lt.3)go to 95
      do 90 i=2,imxm1
      k=1
      j=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvy(i,j,k)/dy+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 91
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
   91 continue
c      write(6,*)' continue 91'
      j=jmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=2.0e0*rdx2 +rdy2 +3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvy(i,j,k)/dy+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 92
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+ cvz(i,j,k)/dz)*gl(mb)
   92 continue
c      write(6,*)' continue 92'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvy(i,j,k)/dy-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 93
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
   93 continue
c      write(6,*)' continue 93'
      j=1
      m=(j-1)*bwo+(i-1)*bwi +k
      acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvy(i,j,k)/dy-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 94
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
   94 continue
c      write(6,*)' continue 94'
  90  continue
c      write(6,*)' continue 90'
  95  continue
c      write(6,*)' continue 95'
      if(jmx.lt.3)go to 105
      do 100 j=2,jmxm1
      k=1
      i=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+ 2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
     # -0.5e0*cvx(i,j,k)/dx+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 101
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  101 continue
c      write(6,*)' continue 101'
      i=imx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
     # +0.5e0*cvx(i,j,k)/dx+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 102
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  102 continue
c      write(6,*)' continue 102'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
     # +0.5e0*cvx(i,j,k)/dx-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 103
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  103 continue
c      write(6,*)' continue 103'
      i=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
     # -0.5e0*cvx(i,j,k)/dx-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 104
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  104 continue
c      write(6,*)' continue 104'
 100  continue
c      write(6,*)' continue 100'
 105  continue
c      write(6,*)' continue 105'
      i=1
      j=1
      k=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
     # +0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dx
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 106
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  106 continue
c      write(6,*)' continue 106'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
     # -0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6. eq. 1) go to 107
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  107 continue
c      write(6,*)' continue 107'
      i=1
      j=jmx
      k=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
     # +0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 108
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  108 continue
c      write(6,*)' continue 108'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  -0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
     # -0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 109
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  109 continue
c      write(6,*)' continue 109'
      i=imx
      j=1
      k=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
     # +0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 111
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  111 continue
c      write(6,*)' continue 111'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
     # -0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 112
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  112 continue
c      write(6,*)' continue 112'
      i=imx
      j=jmx
      k=1
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
     # +0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
      if (ineum5 .eq. 1) go to 113
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
  113 continue
c      write(6,*)' continue 113'
      k=kmx
      m=i*ni+j*nj+k*nk+nc
      acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
     #  +0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
     # -0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
      acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
      acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
      acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
      if (ineum6 .eq. 1) go to 114
      mb=i*nib+j*njb+ncb
      f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
  114 continue
c      write(6,*)' continue 114'
       if(ineum5*ineum6 .eq. 0) go to 117
       if( ising .eq. 1 ) go to 117
       i=1
       j=1
       k=1
       m=i*ni+j*nj+k*nk+nc
       acof2(m)=0.e0
       acof4(m)=0.e0
       acof6(m)=0.e0
       f(m)=0.e0
       i=2
       m=i*ni+j*nj+k*nk+nc
       acof1(m)=0.e0
       i=1
       j=2
       m=i*ni+j*nj+k*nk+nc
       acof3(m)=0.e0
       j=1
       k=2
       m=i*ni+j*nj+k*nk+nc
       acof5(m)=0.e0
 117    continue
c      write(6,*)' continue 117'
c  load coef
c      write(6,*)'load matrix'
      do 150 i=1+bwo,mx
             j=i-bwo
             t=acof3(i)
             call put(t,a,mx,ia,ja,i,j)
 150   continue
c      write(6,*)' continue 150'
       do 155 i=1+bwi,mx
              j=i-bwi
              t=acof1(i)
              call put(t,a,mx,ia,ja,i,j)
 155   continue
c      write(6,*)' continue 155'
      do 160 i=2,mx
         j=i-1
         t=acof5(i)
         call put(t,a,mx,ia,ja,i,j)
 160   continue
c      write(6,*)' continue 160'
       do 165 i=1,mx
        j=i
        t=acof0(i)
       call put(t,a,mx,ia,ja,i,j)
 165   continue
c      write(6,*)' continue 165'
       do 170 i=1,mx-1
          j=i+1
          t=acof6(i)
          call put(t,a,mx,ia,ja,i,j)
 170   continue
c      write(6,*)' continue 170'
       do 175 i=1,mx-bwi
          j=i+bwi
          t=acof2(i)
          call put(t,a,mx,ia,ja,i,j)
 175   continue
c      write(6,*)' continue 175'
       do 180 i=1,mx-bwo
          j=i+bwo
          t=acof4(i)
          call put(t,a,mx,ia,ja,i,j)
 180   continue
c      write(6,*)' continue 180'
c   load zeros for fill-in stripes
       if(iflo .eq. 0) go to 195
       do 192 k = 1,iflo
       do 185 i=1+bwo-k,mx
          j=i-bwo+k
          t=0.0e0
          call put(t,a,mx,ia,ja,i,j)
 185   continue
c      write(6,*)' continue 185'
      do 190 i=1,mx-bwo+k
         j=i+bwo-k
         t=0.0e0
         call put(t,a,mx,ia,ja,i,j)
 190   continue
 192   continue
c      write(6,*)' continue 190'
 195   continue
c      write(6,*)' continue 195'
       if(ifli .eq. 0) go to 210
       do 208 k = 1,ifli
       do 200 i=1+bwi-k,mx
       j=i-bwi+k
       t=0.0e0
       call put(t,a,mx,ia,ja,i,j)
 200   continue
c      write(6,*)' continue 200'
       do 205 i=1,mx-bwi+k
          j=i+bwi-k
          t=0.0e0
          call put(t,a,mx,ia,ja,i,j)
 205   continue
 208   continue
c      write(6,*)' continue 205'
 210   continue
c      write(6,*)' continue 210'
       if( iflm .eq. 0 ) go to 219
       do 217 k = 1,iflm
          do 212 i = 1+bwi+k,mx
             j = i - bwi - k
             call put(0.0e0,a,mx,ia,ja,i,j)
  212     continue
          do 213 i = 1,mx-bwi-k
             j = i + bwi + k
             call put(0.0e0,a,mx,ia,ja,i,j)
  213     continue
  217 continue
  219 continue
c   have now loaded the a,ia,ja arrays
c      write(6,*) 'formula,actual',nonzt,ia(mx+1)
      go to 215
 220  write(6,*) 'ifli or iflo too large'
 215  continue
c      write(6,*)' continue 215'
      return
      end
      subroutine prbtyp(ineum5,ineum6,ising,imx,imy,imz,iflg)
       integer ineum5,ineum6,ising,imx,imy,imz,iflg
       integer n5,n6,isg,ix,iy,iz
c  ineum5=0,dir on bot; =1,neuman on bot
c  ineum6=0,dir on top; =1,neuman on top
c  ising =0,regular; =1 singular if neuman on top & bot
c  imx = # mesh cells in x dir
c  imy = # mesh cells in y dir
c  imz = # mesh cells in z dir
c  iflg =0 sets local values;=1 returns local values
c
c
      if(iflg .eq. 0) go to 10
c      write(6,*)' in prbtyp ix,iy,iz recall',ix,iy,iz
      ineum5=n5
      ineum6=n6
      ising=isg
      imx=ix
      imy=iy
      imz=iz
c      write(6,*)' in prbtyp recall imx,imy,imz',imx,imy,imz
      go to 20
 10   continue
c      write(6,*)' continue #'
c      write(6,*)' in prbtyp init imx,imy,imz',imx,imy,imz
      n5=ineum5
      n6=ineum6
      isg=ising
      ix=imx
      iy=imy
      iz=imz
c      write(6,*)   ' *********** problem type *************'
c      if( ineum5 .eq. 0 ) then
c         write(6,*)' *****       dir on bottom        *****'
c      else
c         write(6,*)' *****     neuman on bottom       *****' 
c      endif
c      if( ineum6 .eq. 0 ) then
c         write(6,*)' *****         dir on top         *****'
c      else
c         write(6,*)' *****      neuman on top         *****'
c      endif
c      if( ising .eq. 0 ) then
c         write(6,*)' *****         regular            *****'
c      else
c         write(6,*)' ***** singular if neuman on top and bottom *****'
c      endif
c      write(6,*)   ' ***** mesh cells in x dir',imx,'     *****'
c      write(6,*)   ' ***** mesh cells in y dir',imy,'     *****'
c      write(6,*)   ' ***** mesh cells in z dir',imz,'     *****'
 20   continue
c      write(6,*)' continue #'
      return
      end
      subroutine iccglu(a,n,ia,ja,af,b,x,r,p,s,tol,maxits,its,job,info)
c
      integer n,ia(1),ja(1),maxits,its,job,info
      real a(1),af(1),x(1),r(1),p(1),s(1),b(1),tol
c
c     this routine performs preconditioned conjugate gradient on a
c     sparse matrix. the preconditioner is an incomplete lu of
c     the matrix.
c
c     on entry:
c
c        a  real()
c           contains the elements of the matrix.
c
c        n  integer
c           is the order of the matrix.
c
c        ia integer(n+1)
c           contains pointers to the start of the rows in the arrays
c           a and ja.
c
c        ja integer()
c           contains the column location of the corresponding elements
c           in the array a.
c
c        b real (n)
c           contains the right hand side.
c
c        x real (n)
c           contains an estimate of the solution, the closer it
c           is to the true solution the faster tthe method will
c           converge.
c
c        tol real
c           is the accuracy desired in the solution.
c
c        maxits integer
c           is the maximun number of iterations to be taken
c           by the routine.
c
c        job integer
c           is a flag to signal if incomplete factorization
c           aready exits in array af.
c           job = 0  perform incomplete factorization
c           job = 1  skip incomplete factorization
c
c     on output
c
c        af real ()
c           contains the incomplete factorization of the matrix
c           contained in the array a.
c
c        x real (n)
c          contains the solution.
c
c        r,p,s real (n)
c          these are scratch work arrays.
c
c        its integer
c          contains the number of iterations need to converge.
c
c        info integer
c          signals if normal termination.
c          info = 0 method converged in its iterations
c          info = 1 method converged, but exit occurred because
c                   residual norm was less than sqrt(rtxmin).
c          info = -999 method didnot converge in maxits iterations.
c
c
c     the algorithm has the following form.
c
c        form incomplete factors l and u
c        x(0) <- initial estiate
c        r(0) <- b - a*x(0)
c        r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0)
c        p(0) <- r(0)
c        i    <- 0
c        while r(i) > tol do
c           s      <- inv(u)*inv(l)*a*p(i)
c           a(i)   <- trans(r(i))*r(i)/(trans(s)*s)
c           x(i+1) <- x(i) + a(i)*p(i)
c           r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s
c           b(i)   <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i))
c           p(i+1) <- r(i+1) + b(i)*p(i)
c           i      <- i + 1
c        end
c
      real ai,bi,rowold,rownew,xnrm,anrm
      real sdot
      real rtxmax,rtxmin
      common /gear14/ rtxmax,rtxmin
      data rtxmin/1.0e-16/
      info = 0
      anrm = abs(a(isamax(ia(n+1)-1,a,1)))
c
c     form incomplete factors l and u
c
      if( job .ne. 0 ) go to 5
         call scopy(ia(n+1)-1,a,1,af,1)
         call lu(af,n,ia,ja)
    5 continue
c
c     r(0) <- b - a*x(0)
c
      call cgres(a,n,ia,ja,x,b,r)
c
c     r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0)
c
c     inv(u)*inv(l)*r(0)
c
      call scopy(n,r,1,s,1)
      call ssol(af,n,ia,ja,s,1)
      call ssol(af,n,ia,ja,s,2)
c
c     invtrans(l)*invtrans(u)*above
c
      call ssol(af,n,ia,ja,s,-2)
      call ssol(af,n,ia,ja,s,-1)
c
c     trans(a)*above
c
      call mmult(a,n,ia,ja,r,s,-1)
c
c     p(0) <- r(0)
c
      call scopy(n,r,1,p,1)
      rowold = sdot(n,r,1,r,1)
      i = 0
c
c     while r(i) > tol do
c
      ai = 1.0d0
   10 continue
      xnrm = abs(x(isamax(n,x,1)))
cc    write(6,*) ' iter residual xnrm ai',i,sqrt(rowold)/(anrm*xnrm),
cc   $            xnrm,ai
      if( sqrt(rowold) .le. tol*(anrm*xnrm) ) go to 12
cc    if (rowold .le. rtxmin) go to 25
   13 continue
c
c        s      <- inv(u)*inv(l)*a*p(i)
c
         call mmult(a,n,ia,ja,s,p,1)
         call ssol(af,n,ia,ja,s,1)
         call ssol(af,n,ia,ja,s,2)
c
c        a(i)   <- trans(r(i))*r(i)/(trans(s)*s)
c
         ai = rowold/sdot(n,s,1,s,1)
c
c        x(i+1) <- x(i) + a(i)*p(i)
c
         call saxpy(n,ai,p,1,x,1)
c
c
c        r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s
c
         call ssol(af,n,ia,ja,s,-2)
         call ssol(af,n,ia,ja,s,-1)
         call mmult(a,n,ia,ja,b,s,-1)
         call saxpy(n,-ai,b,1,r,1)
c
c        b(i)   <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i))
c
         rownew = sdot(n,r,1,r,1)
         bi = rownew/rowold
c
c        p(i+1) <- r(i+1) + b(i)*p(i)
c
         call saxpy2(n,bi,p,1,r,1)
         rowold = rownew
c
c        i      <- i + 1
c
         i = i + 1
         if( i .gt. maxits ) go to 20
         go to 10
   12 continue
   15 continue
      its = i
      return
   20 continue
      info = -999
      its = maxits
      return
c  25 continue
c     info = 1
c     its = i
c     return
      end
      subroutine axpy(a,n,ia,ja,i,js,je,t,y)
c
      integer n,ia(1),ja(1),i,js,je
      real a(1),y(1)
      real t
c
c     this routine computes an axpy for a row of a sparse matrix
c     with a vector.
c     an axpy is a multiple of a row of the matrix is added to the
c     vector y.
c
      is = ia(i)
      ie = ia(i+1) - 1
      do 10 ir = is,ie
         j = ja(ir)
         if( j .gt. je ) go to 20
         if( j .lt. js ) go to 10
         y(j) = y(j) + t*a(ir)
   10 continue
   20 continue
      return
      end
      subroutine saxpy2(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real dx(1),dy(1),da
      integer i,incx,incy,m,mp1,n  
c
      if(n.le.0)return
      if (da .eq. 0.0e0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dx(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dx(i) = dy(i) + da*dx(i)
        dx(i + 1) = dy(i + 1) + da*dx(i + 1)
        dx(i + 2) = dy(i + 2) + da*dx(i + 2)
        dx(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      real function dot(a,n,ia,ja,i,js,je,b)
c
      integer n,ia(1),ja(1),i,js,je
      real a(1),b(1)
      real t
c
c     this routine computes an inner product for a row of a sparse matri
c     with a vector.
c
      t = 0.0e0
      is = ia(i)
      ie = ia(i+1) - 1
      do 10 ir = is,ie
         j = ja(ir)
         if( j .gt. je ) go to 20
         if( j .lt. js ) go to 10
         t = t + a(ir)*b(j)
   10 continue
   20 continue
      dot = t
      return
      end
      subroutine inits(a,n,ia,ja)
c
      integer n,ia(1),ja(1)
      real a(1)
c
c     this routine initializes storage for the sparse
c     matrix arrays.
c     note: the matrix is initialized to have zeroes
c     on the diagonal.
c
      do 10 i = 1,n
         a(i) = 0.0e0
         ja(i) = i
         ia(i) = i
   10 continue
      ia(n+1) = n+1
      return
      end
      subroutine insert(t,a,n,ia,ja,i,j,k)
c
      integer n,ia(1),ja(1),i,j,k
      integer ip1,np1
      real t,a(1)
c
c     this routine rearranges the elements in arrays a and ja
c     and updates array ia for the new element.
c
cc    write(6,1000)i,j,ia(i),ia(i+1),t
1000  format('  *** from insert i,j,ia(i),ia(i+1),t ',4i5,1x,e15.8)
      l1 = k
      l2 = ia(n+1) - 1
      do 10 lb = l1,l2
         l = l2 - lb + l1
         a(l+1) = a(l)
         ja(l+1) = ja(l)
   10 continue
      a(k) = t
      ja(k) = j
      ip1 = i + 1
      np1 = n + 1
      do 20 l = ip1,np1
         ia(l) = ia(l) + 1
   20 continue
      return
      end
      logical function locate(a,n,ia,ja,i,j,k)
c
      integer n,ia(1),ja(1),i,j
      real a(1)
c
c     this routine will locate the i,j-th element in the
c     sparse matrix structure.
c
      is = ia(i)
      ie = ia(i+1) - 1
c
c     row i is from is to ie in array a.
c
      do 10 l = is,ie
         if( j .gt. ja(l) ) go to 10
         if( j .ne. ja(l) ) go to 5
            locate = .true.
            k = l
            go to 20
    5    continue
         if( j .ge. ja(l) ) go to 10
            locate = .false.
            k = 0
         go to 20
   10 continue
c
c     get here if should be at the end of a row.
c
      locate = .false.
      k = 0
   20 continue
      return
      end
      subroutine lu(a,n,ia,ja)
c
      logical locate
      integer n,ia(1),ja(1)
      integer ikp1,kp1
      real a(1)
c
c     this subroutine does incomplete gaussian elimenation
c     on a sparse matrix. the matrix is stored in a sparse
c     data structure.
c     note: no pivoting is done
c           and the factorization is incomplete,
c           i.e., only where there exists a storage location
c           will the operation take place.
c
      do 40 k = 1,n
         if( .not. locate(a,n,ia,ja,k,k,l2) ) go to 50
         if( a(l2) .eq. 0.0e0 ) go to 50
         kp1 = k + 1
         do 30 i = kp1,n
            if( .not. locate(a,n,ia,ja,i,k,ik) ) go to 25
               is = ik
               ie = ia(i+1) - 1
               kj = l2
               ke = ia(k+1) - 1
               a(ik) = a(ik)/a(kj)
               if( kj .eq. ke ) go to 30
               kj = kj + 1
               ikp1 = ik + 1
               do 20 j = ikp1,ie
   10             continue
                     if( kj .gt. ke ) go to 30
                     if( ja(kj) .ge. ja(j) ) go to 15
                        kj = kj + 1
                        go to 10
   15                continue
                  if( ja(kj) .gt. ja(j) ) go to 20
                  a(j) = a(j) - a(ik)*a(kj)
   20          continue
   25       continue
   30    continue
   40 continue
      return
   50 continue
cc    write(6,*)' value of k and a(k,k)',k,a(l2)
cc    write(6,*)' error zero on diagonal'
cc    write(6,*)' matrix probably specified incorrectly or'
cc    write(6,*)' incomplete factorization produces a singular matrix'
      return
      end
      subroutine mmult(a,n,ia,ja,b,x,job)
c
      integer n,ia(1),ja(1),job
      real a(1),b(1),x(1)
c
c     this routine performs matrix vector multiple
c     for a sparse matrix structure
c
c     job has the following input meanings:
c         job = 1  matrix vector multiple
c             = -1 matrix transpose vector multiple
c             = 2  unit lower matrix vector multiple
c             = -2 unit lower matrix transpose vector multiple
c             = 3  upper matrix vector multiple
c             = -3 upper matrix transpose vector multiple
      real dot
c
c     a*x
c
      if( job .ne. 1 ) go to 15
         do 10 i = 1,n
            b(i) = dot(a,n,ia,ja,i,1,n,x)
   10    continue
c
c     trans(a)*x
c
   15 continue
      if( job .ne. -1 ) go to 35
         do 20 i = 1,n
            b(i) = 0.0e0
   20    continue
         do 30 i = 1,n
            call axpy(a,n,ia,ja,i,1,n,x(i),b)
   30    continue
c
c     l*x   when l is unit lower
c
   35 continue
      if( job .ne. 2 ) go to 45
         do 40 i = 1,n
            b(i) = x(i) + dot(a,n,ia,ja,i,1,i-1,x)
   40    continue
c
c     trans(l)*x   when l is unit lower
c
   45 continue
      if( job .ne. -2 ) go to 55
         do 49 i = 1,n
            b(i) = x(i)
   49    continue
         do 50 i = 1,n
            call axpy(a,n,ia,ja,i,i,n,x(i),b)
   50    continue
c
c     u*x
c
   55 continue
      if( job .ne. 3 ) go to 65
         do 60 i = 1,n
            b(i) = dot(a,n,ia,ja,i,i,n,x)
   60    continue
c
c     trans(u)*x
c
   65 continue
      if( job .ne. -3 ) go to 85
         do 70 i = 1,n
            b(i) = 0.0e0
   70    continue
         do 80 i = 1,n
            call axpy(a,n,ia,ja,i,1,i,x(i),b)
   80    continue
   85 continue
      return
      end
      subroutine put(t,a,n,ia,ja,i,j)
c
      integer n,ia(1),ja(1),i,j
      real t,a(1)
c
c     this routine will insert an element into the sparse matrix
c     structure.
c
      is = ia(i)
      ie = ia(i+1) - 1
cc    write(6,100)i,j,ia(i),ia(i+1)
100   format('  *** from put i,j,ia(i),ia(i+1) ',4i7)
c
c     row i is from is to ie in array a.
c
      do 10 k = is,ie
         if( j .gt. ja(k) ) go to 10
         if( j .ne. ja(k) ) go to 5
            a(k) = t
            go to 20
    5    continue
         if( j .ge. ja(k) ) go to 12
            call insert(t,a,n,ia,ja,i,j,k)
            go to 20
   12    continue
   10 continue
c
c     get here if should be at the end of a row.
c
      k = ie + 1
      call insert(t,a,n,ia,ja,i,j,k)
      go to 30
   20 continue
   30 continue
      return
      end
      subroutine cgres(a,n,ia,ja,x,b,r)
c
      integer n,ia(1),ja(1)
      real a(1),x(1),b(1),r(1)
c
c     this routine computes a residual for a*x=b where
c     a is in a sparse structure
c
      real dot
      do 10 i = 1,n
         r(i) = b(i) - dot(a,n,ia,ja,i,1,n,x)
   10 continue
      return
      end
      subroutine ssol(a,n,ia,ja,b,job)
c
      integer n,ia(1),ja(1),job
      real a(1),b(1)
c
c     this routine solves a system of equations based on a sparse
c     matrix date structure.
c     the array b contains the right hand side on input
c         and on output has the solution
c
c     job has the value 1 if l*x = b is to be solved.
c     job has the value -1 if trans(l)*x = b is to be solved.
c     job has the value 2 if u*x = b is to be solved.
c     job has the value -2 if trans(u)*x = b is to be solved.
c
      logical locate
      real t
      real dot
c
c     job = 1  solve  l*x = b
c
      if( job .ne. 1 ) go to 15
c
c        solve l*y=b
c
         do 10 i = 2,n
            b(i) = b(i) - dot(a,n,ia,ja,i,1,i-1,b)
   10    continue
   15 continue
      if( job .ne. 2 ) go to 25
c
c        solve u*x=y
c
         do 20 ib = 1,n
            i = n - ib + 1
            t = dot(a,n,ia,ja,i,i+1,n,b)
            if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30
            b(i) = (b(i) - t)/a(k)
   20    continue
c
c     job = -2  solve  trans(u)*x = b
c
   25 continue
      if( job .ne. -2 ) go to 35
c
c        solve trans(u)*y=b
c
         do 21 i = 1,n
            if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30
            b(i) = b(i)/a(k)
            call axpy(a,n,ia,ja,i,i+1,n,-b(i),b)
   21    continue
c
c        solve trans(l)*x=y
c
   35 continue
      if( job .ne. -1 ) go to 45
         do 22 ib = 2,n
            i = n - ib + 2
            call axpy(a,n,ia,ja,i,1,i-1,-b(i),b)
   22    continue
   45 continue
      return
   30 continue
cc    write(6,*)' error no diagonal element: from solve'
      return
      end
      subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),sa
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sy(i) + sa*sx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        sy(i) = sy(i) + sa*sx(i)
        sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
        sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
        sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
   50 continue
      return
      end
      subroutine scopy(n,sx,incx,sy,incy)
c
c     copies a vector, x, to a vector, y.
c     uses unrolled loops for increments equal to 1.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1)
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        sy(i) = sx(i)
        sy(i + 1) = sx(i + 1)
        sy(i + 2) = sx(i + 2)
        sy(i + 3) = sx(i + 3)
        sy(i + 4) = sx(i + 4)
        sy(i + 5) = sx(i + 5)
        sy(i + 6) = sx(i + 6)
   50 continue
      return
      end
      real function sdot(n,sx,incx,sy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp = stemp + sx(ix)*sy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      sdot = stemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + sx(i)*sy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
     *   sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
   50 continue
   60 sdot = stemp
      return
      end
      integer function isamax(n,sx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),smax
      integer i,incx,ix,n
c
      isamax = 0
      if( n .lt. 1 ) return
      isamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      smax = abs(sx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(abs(sx(ix)).le.smax) go to 5
         isamax = i
         smax = abs(sx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 smax = abs(sx(1))
      do 30 i = 2,n
         if(abs(sx(i)).le.smax) go to 30
         isamax = i
         smax = abs(sx(i))
   30 continue
      return
      end