program main

c*********************************************************************72
c
cc MAIN generates a Markov chain matrix to test eigenvalue routines.
c
c  The matrix A produced by this routine can also be used to generate
c  a related singular matrix, I-A.
c  The matrix models  simple random walk on a triangular grid.
c  see additional comments in subroutine.
c
c   will create a matrix in the HARWELL/BOEING format and put it in
c   the file markov.mat
c
      implicit none

      integer nmax
      parameter ( nmax = 5000 )
      integer nzmax
      parameter ( nzmax= 4 * nmax )

      double precision a(nzmax)
      integer ia(nmax+1)
      integer ifmt
      integer ios
      integer iout
      integer ja(nzmax)
      integer job
      character * ( 8 ) key
      integer m
      integer n
      character * ( 72 ) title
      character * ( 3 ) type
      double precision rhs(1)

      call timestamp ( )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSEKIT_PRB07'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Generate a Markov chain matrix to test the'
      write ( *, '(a)' ) '  eigenvalue routine.'

      open ( unit = 11, file = 'markov.mat', status = 'replace', 
     &  err = 99 )
c
c  Set grid size - will not accept too large grids.
c
      m = 5

      if ( 2 * nmax < m * ( m + 1 ) ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSEKIT_PRB07 - Fatal error!'
        write ( *, '(a)' ) 
     &    '  M is too large - unable to produce matrix.'
        stop
      end if
c
c  Call the matrix generator.
c
      call markgen ( m, n, a, ja, ia )
c
c  Store result in file.
c
      title = ' Test matrix from SPARSKIT - Markov chain model'
      key = 'randwk01'
      type = 'rua'
      iout = 11
      job = 2
      ifmt = 10

      call prtmt ( n, n, a, ja, ia, rhs, 'NN', title, key, 
     &  type, ifmt, job, iout )

      close ( unit = iout )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '  The matrix has been stored in Harwell/Boeing format'
      write ( *, '(a)' ) '  in the file "markov.mat".'

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSEKIT_PRB07'
      write ( *, '(a)' ) '  Normal end of execution.'

      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop

99    continue

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSEKIT_PRB07 - Fatal error!'
      write ( *, '(a)' ) '  Could not open the file.'
      write ( *, '(a)' ) '  Abnormal end of execution.'

      write ( *, '(a)' ) ' '
      call timestamp ( )
      stop

      end
      function afun ( x, y, z )

c*********************************************************************72
c
cc AFUN
c
      implicit none

      double precision afun
      double precision x
      double precision y
      double precision z

      afun = -1.0D+00

      return
      end
      function bfun ( x, y, z )

c*********************************************************************72
c
cc BFUN
c
      implicit none

      double precision bfun
      double precision x
      double precision y
      double precision z

      bfun = -1.0D+00

      return
      end
      function cfun ( x, y, z )

c*********************************************************************72
c
cc CFUN
c
      implicit none

      double precision cfun
      double precision x
      double precision y
      double precision z

      cfun = -1.0D+00

      return
      end
      function dfun ( x, y, z )

c*********************************************************************72
c
cc DFUN
c
      implicit none

      double precision dfun
      double precision x
      double precision y
      double precision z

      dfun = 0.0D+00

      return
      end
      function efun ( x, y, z )

c*********************************************************************72
c
cc EFUN
c
      implicit none

      double precision efun
      double precision x
      double precision y
      double precision z

      efun = 0.0D+00

      return
      end
      function ffun ( x, y, z )

c*********************************************************************72
c
cc FFUN
c
      implicit none

      double precision ffun
      double precision x
      double precision y
      double precision z

      ffun = 0.0D+00

      return
      end
      function gfun ( x, y, z )

c*********************************************************************72
c
cc GFUN
c
      implicit none

      double precision gfun
      double precision x
      double precision y
      double precision z

      gfun = 0.0D+00

      return
      end
      subroutine afunbl ( nfree, x, y, z, coeff )

c*********************************************************************72
c
cc AFUNBL ???
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine bfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine cfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine dfunbl (nfree,x,y,z,coeff)
 
c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine efunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine ffunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine gfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end