c 
c Example for using ROWMAP with left preconditioning.
c 
c (For references see file "rowmap_pc.f".)
c 
c An ODE system is generated from the following two-dimensional 
c nonlinear diffusion equation ("NILIDI")
c
c     du           [  d^2 u     d^2 u  ]
c     --  = exp(u) [  -----  +  -----  ] + u*(18 exp(u) - 1),
c     dt           [  dx^2      dy^2   ]
c
c     t in [0,1], (x,y) in [0,pi/3]^2
c
c with Dirichlet boundary conditions and initial values chosen
c to yield the exact solution
c
c     u(t,x,y) = exp(-t) * sin(3x) * sin(3y).
c
c The PDE system is treated by central differences on a uniform
c mesh with 69 x 69 inner grid points.
c
c ROWMAP uses BLAS Level 1 and 2 and the preconditioner below
c uses in addition LAPACK for solving linear equations. Therefore,
c you have to link these libraries. 
c 
c link:
c       rowmap driver lapack blas
c  
c
c Version of August 25, 1997.
c
c-----!----------------------------------------------------------------- 
c
      program driver
      implicit none
      integer n,liwork,lwork,nx,ipar(1),idid,ipre
      parameter(liwork=100, lwork=500000, nx=5000)
      real*8 u(nx),t,tend,atol,rtol,work(lwork),hs,rpar(1)
      external fex,jacv,solout,fdt,pfac,psol
      integer iwork(liwork),itol,ifcn,ijacv,ifdt,iout,i

c  set the initial values
      call init(n,u,t,tend,ifcn,rpar,ipar)
c  tolerances
      atol=1d-2
      rtol=1d-2
c  atol and rtol are scalars.
      itol=0
c  use default initial step size
      hs=0d0
c  use finite differences for jacobian-vector products
      ijacv=0
c  problem is autonomous, don't call a subroutine computing the
c  derivative f_t:
      ifdt=0
c  don't call solout
      iout=0 
c  use preconditioning
      ipre=1
c  clear work and iwork, i.e. use the default values for optional input.
      do 100 i=1,20
        iwork(i)=0
        work(i)=0d0
 100  continue
      rpar(1)=0d0
      ipar(1)=0

c
c Call to ROWMAP.
c 
      call rowmap(n,fex,ifcn,t,u,tend,hs,rtol,atol,itol,
     1                  jacv,ijacv,fdt,ifdt,solout,iout,
     2                  pfac,psol,ipre,work,
     3                  lwork,iwork,liwork,rpar,ipar,idid)


      write (*,9990) iwork(5),iwork(6),iwork(7),iwork(8),
     1                float(iwork(8))/float(iwork(5))
      write (*,9992) iwork(9),iwork(10)

 9990 format (/ 'Statistics: '/,
     1 'number of steps     =',i5,5x,'number of rejected steps  =',i5/,
     2 'number of f evals.  =',i5,5x,'number of jac-vector op.  =',i5/,
     3 'average krylov dim. =',f7.2)
 9992 format(
     1 'minimum for lwork   =',i7,3x,'minimum for liwork        =',i5/)

      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   S O L O U T
c-----!-----------------------------------------------------------------
c This subroutine is for output of the numerical solution. It applies
c the dense output subroutine "ROWCON" at times t=0.1,0.2,..,1.0;
c ROWMAP calls this subroutine solout only if "iout" is set equal to 1 above.
c
      subroutine solout(n,told,tnew,uold,unew,fold,fnew,ucon,
     1                  intr,rpar,ipar)
      implicit none
      integer n,intr,ipar(*)
      real*8 told,tnew,uold(n),unew(n),fold(n),fnew(n),ucon(n)
      real*8 rpar(*),s
      s=rpar(1)
      if(s.eq.0d0) write (6,9990)  
      if (s.ge.tnew) return
c Call for dense output.
1     call rowcon(n,s,told,tnew,uold,unew,fold,fnew,ucon)
      write (6,9991) s,ucon(1)
      s=s+1d-1
      if (s.lt.tnew) goto 1
c Save last evaluation point:
      rpar(1)=s
 9990 format ('Solution u(t,x,y) at grid point (1,1) at time ' )
 9991 format ('t=',f5.1,2x,'is',2x,d12.6) 
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   F D T 
c-----!-----------------------------------------------------------------
c A dummy only, the problem is autonomous.
c
      subroutine fdt(n,t,u,ft,rpar,ipar) 
      real*8 t,u(n),ft(n),rpar(*)
      integer ipar(*)
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   F E X 
c-----!-----------------------------------------------------------------
c The semidiscrete problem.
c
      subroutine fex(n,t,u,fw,rpar,ipar)
      implicit none
      integer n,m, i,j,ipar(*)
      parameter (m=69)
      real*8 t,u(m,m),fw(m,m), pi,x1,x2,fk2,eu,dx,d2,rpar(*)
      data pi/ 3.141592654D0/
      fk2 = (m+1)*(m+1)*9D0/(pi*pi)
      dx = 1D0/(m+1)
      do 100 i=1,m
       x1 = i*dx
       do 100 j=1,m
       x2 = j*dx
       d2 = -4D0*u(i,j)
       if (i.gt.1) d2 = d2+u(i-1,j)
       if (i.lt.m) d2 = d2+u(i+1,j)
       if (j.gt.1) d2 = d2+u(i,j-1)
       if (j.lt.m) d2 = d2+u(i,j+1)
       eu = dexp(u(i,j))
       fw(i,j) = eu*fk2*d2+u(i,j)*(18D0*eu-1D0)
  100 continue
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   J A C V 
c-----!-----------------------------------------------------------------
c This subroutine computes z:=A*y with A=Jacobian of "fex".
c It is only called if "ijacv" is set equal 1 above.
c
      subroutine jacv(n,t,u,y,z,rpar,ipar)
      implicit none
      integer n,m, i,j,ipar(*)
      parameter (m=69)
      real*8 t,u(m,m),y(m,m),z(m,m), pi,fk2,eu,d2,du,rpar(*)
      data pi/3.14159265358979D0/
      fk2 = (m+1)*(m+1)*9D0/(pi*pi)
      do 100 i=1,m
       do 100 j=1,m
       eu = dexp(u(i,j))
       d2 = -4D0*y(i,j)
       du = -4D0*u(i,j)
       if (i.gt.1) then
        d2 = d2+y(i-1,j)
        du = du+u(i-1,j)
       endif
       if (i.lt.m) then
        d2 = d2+y(i+1,j)
        du = du+u(i+1,j)
       endif
       if (j.gt.1) then
        d2 = d2+y(i,j-1)
        du = du+u(i,j-1)
       endif
       if (j.lt.m) then
        d2 = d2+y(i,j+1)
        du = du+u(i,j+1)
       endif
       z(i,j) = eu*fk2*(d2+du*y(i,j))+(18D0*eu*(1D0+u(i,j))-1D0)*y(i,j)
  100 continue
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   I N I T 
c-----!-----------------------------------------------------------------
c
c
      subroutine init(n,u,t0,te,ifcn,rpar,ipar)
      implicit none
      integer n,m, i,j,ifcn,ipar(*)
      parameter (m=69)
      real*8 u(m,m),t0,te, x1,x2,dx,pi,rpar(*)
      data pi/3.14159265358979D0/
c The problem is autonomous.
      ifcn=0
c Number of equations.
      n = m*m
c Initial values.
      t0 = 0D0
      dx = pi/(m+1)
      do 100 i=1,m
       x1 = i*dx
       do 100 j=1,m
       x2 = j*dx
  100  u(i,j) = dsin(x1)*dsin(x2)
c Endpoint of integration.
      te = 1d0
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   P F A C
c-----!-----------------------------------------------------------------
c For preconditioning we use the preconditioner P=P_x*P'_y obtained by
c splitting the diffusion term with respect to x and y as 
c described in reference [4] in "rowmap_pc.f". P_x is a tridiagonal 
c matrix, corresponding to the disrectized operator u_{xx} and the
c half of the reaction term. P_y is defined analogously, the permutation
c P'_y is also tridiagonal.
c Due to the symmetry of the problem it is sufficent to compute
c the LU-decomposition of P_x only. 
c
      subroutine pfac(n,t,y,fw,hg,rpar,ipar)
      implicit none
      integer n,ldd,ipar(*),ml,mu
      integer m,i
      parameter (m=69)

      integer gi(m,m)
      real*8 g(4,m,m),tv(m*m)
      common /cpfac/g,tv,gi

      real*8 t,y(m,m),fw(n),rpar(*),hg,pi
      real*8 u,fk2,eu,up,um,t1,t2,t3
      integer j,info

      data pi /3.141592654D0/
      fk2 = (m+1)*(m+1)*9D0/(pi*pi)
      ml=1
      mu=1
      ldd=4

      call dcopy(ldd*n,0d0,0,g,1)

      do j=1,m
      do i=1,m
        u=y(i,j)
        eu=dexp(u)
        t1=fk2*eu
        um=0d0
        up=0d0
        if (i.gt.1) then
            um=y(i-1,j)
            g(4,i-1,j)=-hg*t1
        end if
        if (i.lt.m) then
            up=y(i+1,j)
            g(2,i+1,j)=-hg*t1
        end if
        t2=t1*(um-2d0*u+up-2d0)
        t3=0.5d0*(18d0*eu-1d0)+0.5d0*u*(18d0*eu)
        g(3,i,j)=1d0-hg*(t2+t3)
      end do
      g(4,m,j)=0d0
      g(2,1,j)=0d0
      call dgbtrf(m,m,1,1,g(1,1,j),ldd,gi(1,j),info)
      if (info.ne.0) then
        write (*,*) 'Error in "pfac": '
        write (*,*) '"dgbtrf" returned  info=',info
        stop
      end if
      end do
      return
      end


c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   P S O L 
c-----!-----------------------------------------------------------------
c Here we solve the linear systems for the preconditioning
c   (P_x) (P'_y) z =  b.
c With the permutation S, we get P_x = P_y := (S-inverse) P'_y S.
c The LU-decomposition of P_x was already computed by "pfac"
c and is stored the array "g". We compute z from b in situ. 
c
      subroutine psol(n,b,rpar,ipar)
      implicit none
      integer n,ipar(*),info,ldd,ml,mu
      integer m
      parameter (m=69)
      integer i,j
      real*8 rpar(*),b(n)

      integer gi(m,m)
      real*8 g(4,m,m),tv(m*m)
      common /cpfac/g,tv,gi

      ml=1
      mu=1
      ldd=4
      info=0
    
      do j=1,m

      call 
     +  dgbtrs('n',m,ml,mu,1,g(1,1,j),ldd,gi(1,j),b(1+(j-1)*m),m,info) 
      if (info.ne.0) then
        write (*,*) 'Error in "psol":'
        write (*,*) '"dgbtrs" returned  info=',info
        write (*,*) 'info=',info
        stop
      end if

      end do

      do 100 i=1,m
      do 100 j=1,m
        tv(i+m*(j-1))=b(m*(i-1)+j)
 100  continue
      do j=1,m
      call 
     +  dgbtrs('n',m,ml,mu,1,g(1,1,j),ldd,gi(1,j),tv(1+(j-1)*m),m,info)
      if (info.ne.0) then
        write (*,*) 'Error in "psol":'
        write (*,*) '"dgbtrs" returned  info=',info
        write (*,*) 'info=',info
        stop
      end if
      end do

      do 110 i=1,m
      do 110 j=1,m
        b(m*(i-1)+j)=tv(i+m*(j-1))
 110  continue
      return
      end