program main

c*********************************************************************72
c
cc nas() is the NAS kernel benchmark program.
c
c  Discussion:
c
c    This is a version of the NAS kernel benchmark program,
c    whose original version was created by David Bailey,
c    dated 17 December 1984.
c
c    Each of the tests begins by filling arrays with pseudorandom values
c    generated by the recursion:
c      x(n+1) = 5^7 * x(n)  (mod 2^30)
c    This recursion will generate 2^28 (approx. 268 million) numbers
c    before repeating.  For this scheme to work properly, the hardware
c    multiply operation must be correct to 47 bits of precision.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      double precision er(8)
      double precision fp(8)
      integer i
      character*8 pn(8)
      double precision rt(8)
      double precision te
      double precision tf
      double precision tm(8)
      double precision tt

      data pn / 
     &  'MXM', 
     &  'CFFT2D', 
     &  'CHOLSKY', 
     &  'BTRIX',
     &  'GMTRY',
     &  'EMIT',
     &  'VPENTA',
     &  'Total'/

      te = 0.0D+00
      tf = 0.0D+00
      tt = 0.0D+00

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NAS:'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '                The NAS kernel benchmark program'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) ' Program        Error          FP Ops'
     &  //'        Seconds     MFLOPS'
      write ( *, '(a)' ) ' '

      do i = 1, 7

        if ( i .eq. 1 ) then
          call mxmtst ( er(1), fp(1), tm(1) )
        else if ( i .eq. 2 ) then
          call ffttst ( er(2), fp(2), tm(2) )
        else if ( i .eq. 3 ) then
          call chotst ( er(3), fp(3), tm(3) )
        else if ( i .eq. 4 ) then
          call btrtst ( er(4), fp(4), tm(4) )
        else if ( i .eq. 5 ) then
          call gmttst ( er(5), fp(5), tm(5) )
        else if ( i .eq. 6 ) then
          call emitst ( er(6), fp(6), tm(6) )
        else if ( i .eq. 7 ) then
          call vpetst ( er(7), fp(7), tm(7) )
        end if

        rt(i) = 1.0D-06 * fp(i) / tm(i)
        write ( *, '(1x,a8,1p2e15.4,0pf12.4,f12.2)' ) 
     &    pn(i), er(i), fp(i), tm(i), rt(i)

        te = te + er(i)
        tf = tf + fp(i)
        tt = tt + tm(i)

      end do

      er(8) = te
      fp(8) = tf
      tm(8) = tt
      rt(8) = 1.0D-06 * tf / tt
      i = 8
      write ( *, '(a)' ) ' '
      write ( *, '(1x,a8,1p2e15.4,0pf12.4,f12.2)' ) 
     &  pn(i), er(i), fp(i), tm(i), rt(i)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NAS:'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      function cptime ( )

c*********************************************************************72
c
cc CPTIME returns a reading of the CPU or elapsed time clock.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      double precision cptime

      cptime = second ( )

      return
      end
      subroutine copy ( n, a, b )

c*********************************************************************72
c
cc COPY copies a double precision vector.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer n

      double precision a(n)
      double precision b(n)
      integer i

      do i = 1, n
        b(i) = a(i)
      end do

      return
      end
      subroutine mxmtst ( er, fp, tm )

c*********************************************************************72
c
cc MXMTST is the test program for MXM.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer l
      integer m
      integer n

      parameter ( l = 256 )
      parameter ( m = 128 )
      parameter ( n = 64 )

      double precision a(l,m)
      double precision ans
      double precision b(m,n)
      double precision c(l,n)
      double precision cptime
      double precision er
      double precision f7
      double precision fp
      integer i
      integer ii
      integer it
      integer j
      double precision t
      double precision t30
      double precision tm

      it = 100
      ans = 35.2026179738722D+00
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t = f7 / t30

      do j = 1, m
        do i = 1, l
          t = dmod ( f7 * t, 1.0D+00 )
          a(i,j) = t
        end do
      end do

      do j = 1, n
        do i = 1, m
          t = dmod ( f7 * t, 1.0D+00 )
          b(i,j) = t
        end do
      end do
c
c  Timing.
c
      tm = cptime ( )

      do ii = 1, it
        call mxm ( a, b, c, l, m, n )
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( c(19,19) - ans ) / ans )
      fp = dble ( 2 * it * l * m * n )

      return
      end
      subroutine mxm ( a, b, c, l, m, n )

c*********************************************************************72
c
cc MXM computes the matrix product C = A * B.
c
c  Discussion:
c
c    The function uses 4-way unrolled loops to carry out matrix multiplication.
c
c    M must be a multiple of 4.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer l
      integer m
      integer n

      double precision a(l,m)
      double precision b(m,n)
      double precision c(l,n)
      integer i
      integer j
      integer k

      do k = 1, n
        do i = 1, l
          c(i,k) = 0.0D+00
        end do
      end do

      do j = 1, m, 4
        do k = 1, n
          do i = 1, l
            c(i,k) = c(i,k) 
     &        + a(i,j)   * b(j,k)
     &        + a(i,j+1) * b(j+1,k) 
     &        + a(i,j+2) * b(j+2,k)
     &        + a(i,j+3) * b(j+3,k)
          end do
        end do
      end do

      return
      end
      subroutine ffttst ( er, fp, tm )

c*********************************************************************72
c
cc FFTTST is the test program for CFFT2D1 and CFFTD2.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer m
      integer m1
      integer n

      parameter ( m = 128 )
      parameter ( m1 = 128 )
      parameter ( n = 256 )

      double precision ans
      double precision cptime
      double precision er
      double precision f7
      double precision fp
      integer i
      integer ip(2*n)
      integer it
      integer j
      integer k
      double precision rmn
      double precision t1
      double precision t2
      double precision t30
      double precision tm
      double complex w1(m)
      double complex w2(n)
      double complex x(m1,n)

      it = 100
      ans = 0.894799941219277D+00
      rmn = 1.0D+00 / dble ( m * n )
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t2 = f7 / t30

      do j = 1, n
        do i = 1, m
          t1 = dmod ( f7 * t2, 1.0D+00 )
          t2 = dmod ( f7 * t1, 1.0D+00 )
          x(i,j) = dcmplx ( t1, t2 )
        end do
      end do

      call cfft2d1 ( 0, m, m1, n, x, w1, ip )
      call cfft2d2 ( 0, m, m1, n, x, w2, ip )
c
c  Timing.
c
      tm = cptime ( )

      do k = 1, it

        do j = 1, n
          do i = 1, m
            x(i,j) = rmn * x(i,j)
          end do
        end do

        call cfft2d1 ( 1, m, m1, n, x, w1, ip )
        call cfft2d2 ( 1, m, m1, n, x, w2, ip )
        call cfft2d2 ( -1, m, m1, n, x, w2, ip )
        call cfft2d1 ( -1, m, m1, n, x, w1, ip )

      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( dble ( x(19,19) ) - ans ) / ans )
      fp = dble ( it * m * n ) * ( 2.0D+00  
     &  + 10.0D+00 * dlog ( dble ( m * n ) ) / dlog ( 2.0D+00 ) )

      return
      end
      subroutine cfft2d1 ( is, m, m1, n, x, w, ip )

c*********************************************************************72
c
cc CFFT2D1 performs complex radix 2 FFT's on the first dimension.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer m
      integer m1
      integer n

      double complex ct
      double complex cx
      integer i
      integer i1
      integer i2
      integer ii
      integer im
      integer ip(2,m)
      integer is
      integer j
      integer k
      integer l
      integer m2
      double precision pi
      double precision t
      double complex w(m)
      double complex x(m1,n)

      pi = 3.141592653589793D+00
c
c  If IS = 0 then initialize only.
c
      m2 = m / 2
      if ( is .eq. 0 ) then
        do i = 1, m2
          t = 2.0D+00 * pi * dble ( i - 1 ) / dble ( m )
          w(i) = dcmplx ( dcos ( t ), dsin ( t ) )
        end do
        return
      end if
c
c  Perform forward or backward FFT's according to IS = 1 or -1.
c
      do i = 1, m
        ip(1,i) = i
      end do

      l = 1
      i1 = 1

120   continue

      i2 = 3 - i1

      do j = l, m2, l

        cx = w(j-l+1)
        if ( is .lt. 0 ) then
          cx = dconjg ( cx )
        end if

        do i = j - l + 1, j
          ii = ip(i1,i)
          ip(i2,i+j-l) = ii
          im = ip(i1,i+m2)
          ip(i2,i+j) = im
          do k = 1, n
            ct = x(ii,k) - x(im,k)
            x(ii,k) = x(ii,k) + x(im,k)
            x(im,k) = ct * cx
          end do
        end do

      end do

      l = 2 * l
      i1 = i2

      if ( l .le. m2 ) then
        go to 120
      end if

      do i = 1, m
        ii = ip(i1,i)
        if ( i .lt. ii ) then
          do k = 1, n
            ct = x(i,k)
            x(i,k) = x(ii,k)
            x(ii,k) = ct
          end do
        end if
      end do

      return
      end
      subroutine cfft2d2 ( is, m, m1, n, x, w, ip )

c*********************************************************************72
c
cc CFFT2D2 performs complex radix 2 FFT's on the second dimension.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer m1
      integer n

      double complex ct
      double complex cx
      integer i
      integer i1
      integer i2
      integer ii
      integer im
      integer ip(2,n)
      integer is
      integer j
      integer k
      integer l
      integer m
      integer n2
      double precision pi
      double precision t
      double complex w(n)
      double complex x(m1,n)

      pi = 3.141592653589793D+00
c
c  If IS = 0, then initialize only.
c
      n2 = n / 2
      if ( is .eq. 0 ) then
        do i = 1, n2
          t = 2.0D+00 * pi * dble ( i - 1 ) / dble ( n )
          w(i) = dcmplx ( dcos ( t ), dsin ( t ) )
        end do
        return
      end if
c
c  Perform forward or backward FFT's according to IS = 1 or -1.
c
      do i = 1, n
        ip(1,i) = i
      end do

      l = 1
      i1 = 1

120   continue

      i2 = 3 - i1

      do j = l, n2, l

        cx = w(j-l+1)
        if ( is .lt. 0 ) then
          cx = dconjg ( cx )
        end if

        do i = j - l + 1, j
          ii = ip(i1,i)
          ip(i2,i+j-l) = ii
          im = ip(i1,i+n2)
          ip(i2,i+j) = im
          do k = 1, m
            ct = x(k,ii) - x(k,im)
            x(k,ii) = x(k,ii) + x(k,im)
            x(k,im) = ct * cx
          end do
        end do
      end do

      l = 2 * l
      i1 = i2

      if ( l .le. n2 ) then
        go to 120
      end if

      do i = 1, n
        ii = ip(i1,i)
        if ( i .lt. ii ) then
          do k = 1, m
            ct = x(k,i)
            x(k,i) = x(k,ii)
            x(k,ii) = ct
          end do
        end if
      end do

      return
      end
      subroutine chotst ( er, fp, tm )

c*********************************************************************72
c
cc CHOTST is the test program for CHOLSKY.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer ida
      integer m
      integer n
      integer nmat
      integer nrhs

      parameter ( ida = 250 )
      parameter ( m = 4 )
      parameter ( n = 40 )
      parameter ( nmat = 250 )
      parameter ( nrhs = 3 )

      double precision a(0:ida,-m:0,0:n)
      double precision ans
      double precision ax(0:ida,-m:0,0:n)
      double precision b(0:nrhs,0:nmat,0:n)
      double precision bx(0:nrhs,0:nmat,0:n)
      double precision cptime
      double precision er
      double precision f7
      double precision fp
      integer i
      integer it
      integer j
      integer k
      integer la
      integer lb
      double precision t
      double precision t30
      double precision tm

      it = 200
      ans = 5177.88531774562D+00
      la = ( ida + 1 ) * ( m + 1 ) * ( n + 1 )
      lb = ( nrhs + 1 ) * ( nmat + 1 ) * ( n + 1 )
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t = f7 / t30

      do k = 0, n
        do j = - m, 0
          do i = 0, ida
            t = dmod ( f7 * t, 1.0D+00 )
            ax(i,j,k) = t
          end do
        end do
      end do

      do k = 0, n
        do j = 0, nmat
          do i = 0, nrhs
            t = dmod ( f7 * t, 1.0D+00 )
            bx(i,j,k) = t
          end do
        end do
      end do
c
c  Timing.
c
      tm = cptime ( )

      do j = 1, it
        call copy ( la, ax, a )
        call copy ( lb, bx, b )
        call cholsky ( ida, nmat, m, n, a, nrhs, ida, b )
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( b(1,19,19) - ans ) / ans )
      fp = dble ( it * ( nmat + 1 ) * 4403 )

      return
      end
      subroutine cholsky ( ida, nmat, m, n, a, nrhs, idb, b )

c*********************************************************************72
c
cc CHOLSKY carries out Cholesky decomposition and back substitution.
c
c  Discussion:
c
c    The Cholesky decomposition is performed on a set of input matrices
c    which are provided as a single three-dimensional array.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer ida
      integer idb
      integer m
      integer n
      integer nrhs

      double precision a(0:ida,-m:0, 0:n)
      double precision b(0:nrhs,0:idb,0:n)
      double precision eps
      double precision epss(0:256)
      integer i
      integer i0
      integer j
      integer jj
      integer k
      integer l
      integer nmat

      eps = 1.0D-13
c
c  Cholesky decomposition.
c
      do j = 0, n

        i0 = max ( -m, -j )
c
c  Off diagonal elements.
c
        do i = i0, -1
          do jj = i0 - i, -1
            do l = 0, nmat
              a(l,i,j) = a(l,i,j) - a(l,jj,i+j) * a(l,i+jj,j)
            end do
          end do
          do l = 0, nmat
            a(l,i,j) = a(l,i,j) * a(l,0,i+j)
          end do
        end do
c
c  Store inverse of diagonal elements.
c
        do l = 0, nmat
          epss(l) = eps * a(l,0,j)
        end do

        do jj = i0, -1
          do l = 0, nmat
            a(l,0,j) = a(l,0,j) - a(l,jj,j) ** 2
          end do
        end do

        do l = 0, nmat
          a(l,0,j) = 1.0D+00 / dsqrt ( dabs ( epss(l) + a(l,0,j) ) )
        end do

      end do
c
c  Solution.
c
      do i = 0, nrhs

        do k = 0, n

          do l = 0, nmat
            b(i,l,k) = b(i,l,k) * a(l,0,k)
          end do

          do jj = 1, min ( m, n - k )
            do l = 0, nmat
              b(i,l,k+jj) = b(i,l,k+jj) - a(l,-jj,k+jj) * b(i,l,k)
            end do
          end do

        end do

        do k = n, 0, -1

          do l = 0, nmat
            b(i,l,k) = b(i,l,k) * a(l,0,k)
          end do

          do jj = 1, min ( m, k )
            do l = 0, nmat
              b(i,l,k-jj) = b(i,l,k-jj) - a(l,-jj,k) * b(i,l,k)
            end do
          end do

        end do

      end do

      return
      end
      subroutine btrtst ( er, fp, tm )

c*********************************************************************72
c
cc BTRTST is the test program for BTRIX.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer jd
      integer kd
      integer ld
      integer md

      parameter ( jd = 30 )
      parameter ( kd = 30 )
      parameter ( ld = 30 )
      parameter ( md = 30 )

      double precision a(5,5,md,md)
      double precision ans
      double precision b(5,5,md,md)
      double precision bx(5,5,md,md)
      double precision c(5,5,md,md)
      double precision cptime
      double precision er
      double precision f7
      double precision fp
      integer i
      integer ii
      integer it
      integer j
      integer je
      integer js
      integer k
      integer l
      integer le
      integer ls
      integer nb
      integer ns
      double precision s(jd,kd,ld,5)
      double precision sx(jd,kd,ld,5)
      double precision t
      double precision t30
      double precision tm

      js = 2 
      je = 29
      ls = 2 
      le = 29
      it = 20
      ans = -0.286282658663962D+00
      nb = 25 * md * md
      ns = jd * kd * ld * 5
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t = f7 / t30

      do l = 1, md
        do k = 1, md
          do j = 1, 5
            do i = 1, 5
              t = dmod ( f7 * t, 1.0D+00 )
              a(i,j,k,l) = t
              t = dmod ( f7 * t, 1.0D+00 )
              bx(i,j,k,l) = t
              t = dmod ( f7 * t, 1.0D+00 )
              c(i,j,k,l) = t
            end do
          end do
        end do
      end do

      do l = 1, 5
        do k = 1, ld
          do j = 1, kd
            do i = 1, jd
              t = dmod ( f7 * t, 1.0D+00 )
              sx(i,j,k,l) = t
            end do
          end do
        end do
      end do
c
c  Timing.
c
      tm = cptime ( )

      do ii = 1, it
        call copy ( ns, sx, s )
        do k = 1, kd
          call copy ( nb, bx, b )
          call btrix ( js, je, ls, le, k, jd, kd, ld, md,
     &      a, b, c, s )
        end do
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( s(19,19,19,1) - ans ) / ans )
      fp = dble ( it * md * ( le - 1 ) * 19165 )

      return
      end
      subroutine btrix ( js, je, ls, le, k, jd, kd, ld, md,
     &  a, b, c, s )

c*********************************************************************72
c
cc BTRIX is a block tridiagonal solver in one direction.
c
c  Discussion:
c
c    The array has four dimensions.  The routine solves along the
c    "J" index.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer jd
      integer kd
      integer ld
      integer md

      double precision a(5,5,md,md)
      double precision b(5,5,md,md)
      double precision c(5,5,md,md)
      double precision c1
      double precision c2
      double precision c3
      double precision c4
      double precision c5
      double precision d1
      double precision d2
      double precision d3
      double precision d4
      double precision d5
      integer j
      integer je
      integer jem1
      integer js
      integer k
      integer l
      integer le
      integer ls
      double precision l11(md)
      double precision l21(md)
      double precision l31(md)
      double precision l41(md)
      double precision l51(md)
      double precision l22(md)
      double precision l32(md)
      double precision l33(md)
      double precision l42(md)
      double precision l43(md)
      double precision l44(md)
      double precision l52(md)
      double precision l53(md)
      double precision l54(md)
      double precision l55(md)
      integer m
      integer n
      double precision s(jd,kd,ld,5)
      double precision u12(md)
      double precision u13(md)
      double precision u14(md)
      double precision u15(md)
      double precision u23(md)
      double precision u24(md)
      double precision u25(md)
      double precision u34(md)
      double precision u35(md)
      double precision u45(md)
c
c  Part 1.  Forward block sweep.
c
      do j = js, je
c
c  Step 1.  Construct L(I) in B.
c
        if ( j .ne. js ) then
          do m = 1, 5
            do n = 1, 5
              do l = ls, le
                b(m,n,j,l) = b(m,n,j,l) 
     &            - a(m,1,j,l) * b(1,n,j-1,l)
     &            - a(m,2,j,l) * b(2,n,j-1,l) 
     &            - a(m,3,j,l) * b(3,n,j-1,l)
     &            - a(m,4,j,l) * b(4,n,j-1,l) 
     &            - a(m,5,j,l) * b(5,n,j-1,l)
              end do
            end do
          end do
        end if
c
c  Step 2.  Compute L inverse.
c
c  A.  Decompose L(I) into L and U.
c
        do l = ls, le
          l11(l) = 1.0D+00 / b(1,1,j,l)
          u12(l) = b(1,2,j,l) * l11(l)
          u13(l) = b(1,3,j,l) * l11(l)
          u14(l) = b(1,4,j,l) * l11(l)
          u15(l) = b(1,5,j,l) * l11(l)
          l21(l) = b(2,1,j,l)
          l22(l) = 1.0D+00 / ( b(2,2,j,l) - l21(l) * u12(l) )
          u23(l) = ( b(2,3,j,l) - l21(l) * u13(l) ) * l22(l)
          u24(l) = ( b(2,4,j,l) - l21(l) * u14(l) ) * l22(l)
          u25(l) = ( b(2,5,j,l) - l21(l) * u15(l) ) * l22(l)
          l31(l) = b(3,1,j,l)
          l32(l) = b(3,2,j,l) - l31(l) * u12(l)
          l33(l) = 1.0D+00 / 
     &      ( b(3,3,j,l) - l31(l) * u13(l) - l32(l) * u23(l) )
          u34(l)  = ( b(3,4,j,l) - l31(l) * u14(l) 
     &      - l32(l) * u24(l) ) * l33(l)
          u35(l)  = ( b(3,5,j,l) - l31(l) * u15(l) 
     &      - l32(l) * u25(l) ) * l33(l)
        end do

        do l = ls, le
          l41(l) = b(4,1,j,l)
          l42(l) = b(4,2,j,l) - l41(l) * u12(l)
          l43(l) = b(4,3,j,l) - l41(l) * u13(l) - l42(l) * u23(l)
          l44(l) = 1.0D+00 / ( b(4,4,j,l) - l41(l) * u14(l) 
     &      - l42(l) * u24(l) - l43(l) * u34(l) )
          u45(l) = ( b(4,5,j,l) - l41(l) * u15(l) - l42(l) * u25(l)
     &      - l43(l) * u35(l) ) * l44(l)
          l51(l) = b(5,1,j,l)
          l52(l) = b(5,2,j,l) - l51(l) * u12(l)
          l53(l) = b(5,3,j,l) - l51(l) * u13(l) - l52(l) * u23(l)
          l54(l) = b(5,4,j,l) - l51(l) * u14(l) - l52(l) * u24(l)
     &      - l53(l) * u34(l)
          l55(l) = 1.0D+00 / 
     &      ( b(5,5,j,l) - l51(l) * u15(l) - l52(l) * u25(l)
     &      - l53(l) * u35(l) - l54(l) * u45(l) )
        end do
c
c  Step 3.  Solve for intermediate vector.
c
c  A.  Construct the right hand side.
c
        if ( j .ne. js ) then
          do m = 1, 5
            do l = ls, le
              s(j,k,l,m) = s(j,k,l,m) 
     &          - a(m,1,j,l) * s(j-1,k,l,1)
     &          - a(m,2,j,l) * s(j-1,k,l,2) 
     &          - a(m,3,j,l) * s(j-1,k,l,3)
     &          - a(m,4,j,l) * s(j-1,k,l,4) 
     &          - a(m,5,j,l) * s(j-1,k,l,5)
            end do
          end do
        end if
c
c  B. intermediate vector.
c
c  Forward substitution.
c
        do l = ls, le
          d1 = s(j,k,l,1) * l11(l)
          d2 = ( s(j,k,l,2) - l21(l) * d1 ) * l22(l)
          d3 = ( s(j,k,l,3) - l31(l) * d1 - l32(l) * d2 ) * l33(l)
          d4 = ( s(j,k,l,4) - l41(l) * d1 - l42(l) * d2 
     &      - l43(l) * d3 ) * l44(l)
          d5 = ( s(j,k,l,5) - l51(l) * d1 - l52(l) * d2 - l53(l) * d3
     &      - l54(l) * d4 ) * l55(l)
c
c  Backward substitution.
c
          s(j,k,l,5) = d5
          s(j,k,l,4) = d4 - u45(l) * d5
          s(j,k,l,3) = d3 - u34(l) * s(j,k,l,4) - u35(l) * d5
          s(j,k,l,2) = d2 - u23(l) * s(j,k,l,3) - u24(l) * s(j,k,l,4)
     &      - u25(l) * d5
          s(j,k,l,1) = d1 - u12(l) * s(j,k,l,2) - u13(l) * s(j,k,l,3)
     &      - u14(l) * s(j,k,l,4) - u15(l) * d5
        end do
c
c  Step 4.  Construct U(I) = inverse(L(I))*C(I+1) by columns and store in B.
c
        if ( j .ne. je ) then

          do n = 1, 5
            do l = ls, le
c
c  Forward substitution.
c
              c1 = c(1,n,j,l) * l11(l)
              c2 = ( c(2,n,j,l) - l21(l) * c1 ) * l22(l)
              c3 = ( c(3,n,j,l) - l31(l) * c1 - l32(l) * c2 ) * l33(l)
              c4 = ( c(4,n,j,l) - l41(l) * c1 - l42(l) * c2 
     &          - l43(l) * c3 ) * l44(l)
              c5 = ( c(5,n,j,l) - l51(l) * c1 - l52(l) * c2 
     &          - l53(l) * c3 - l54(l) * c4 ) * l55(l)
c
c  Backward substitution.
c
              b(5,n,j,l) = c5
              b(4,n,j,l) = c4 - u45(l) * c5
              b(3,n,j,l) = c3 - u34(l) * b(4,n,j,l) - u35(l) * c5
              b(2,n,j,l) = c2 - u23(l) * b(3,n,j,l) 
     &          - u24(l) * b(4,n,j,l) - u25(l) * c5
              b(1,n,j,l) = c1 - u12(l) * b(2,n,j,l) 
     &          - u13(l) * b(3,n,j,l) - u14(l) * b(4,n,j,l) 
     &          - u15(l) * c5
            end do
          end do

        end if

      end do
c
c  Part 2.  Backward block sweep.
c
      jem1 = je - 1

      do j = jem1, js, -1
        do m = 1, 5
          do l = ls, le
            s(j,k,l,m) = s(j,k,l,m) 
     &        - b(m,1,j,l) * s(j+1,k,l,1)
     &        - b(m,2,j,l) * s(j+1,k,l,2)
     &        - b(m,3,j,l) * s(j+1,k,l,3)
     &        - b(m,4,j,l) * s(j+1,k,l,4) 
     &        - b(m,5,j,l) * s(j+1,k,l,5)
          end do
        end do
      end do

      return
      end
      subroutine gmttst ( er, fp, tm )

c*********************************************************************72
c
cc GMTTST is the test program for GMTRY.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer nb
      integer nw

      parameter ( nb = 5 )
      parameter ( nw = 100 )

      double precision ans
      double precision cptime
      double precision er
      double precision f7
      double precision fp
      integer i
      integer it
      integer j
      integer lw
      integer nwall(nb)
      double complex proj(nw,nb)
      double precision rmatrx(nw*nb,nw*nb)
      double precision t1
      double precision t2
      double precision t30
      double precision tm
      double complex wall(nw,nb)
      double precision xmax(nb)
      double complex zcr(nw,nb)

      it = 2
      ans = -2.57754233214174D+00
      lw = 2 * nw * nb
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t2 = f7 / t30

      do j = 1, nb
        nwall(j) = nw
      end do

      do j = 1, nb
        do i = 1, nw
          t1 = dmod ( f7 * t2, 1.0D+00 )
          t2 = dmod ( f7 * t1, 1.0D+00 )
          wall(i,j) = dcmplx ( t1, t2 )
        end do
      end do
c
c  Timing.
c
      tm = cptime ( )

      do i = 1, it
        call gmtry ( nb, nw, nwall, proj, rmatrx, wall, xmax, zcr )
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( rmatrx(19,19) - ans ) / ans )
      fp = dble ( it ) * ( dble ( 120 * ( nb * nw ) ** 2 )
     &  +  0.666D+00 * dble ( nb * nw ) ** 3 )

      return
      end
      subroutine gmtry ( nb, nw, nwall, proj, rmatrx, wall, xmax, zcr )

c*********************************************************************72
c
cc GMTRY computes solid-related arrays.
c
c  Discussion:
c
c    This function was extracted from a vortex method program.
c    It sets up arrays needed for the computation, and performs
c    Gauss elimination on the matrix of wall-influence coefficients.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer nb
      integer nw

      double precision arcl
      double precision dum
      integer i
      integer i0
      integer j
      integer j0
      integer k
      integer kp
      integer kron
      integer ks
      integer l
      integer l1
      integer l2
      integer ls
      integer matdim
      integer nwall(nb)
      double precision period
      double precision pi
      double precision pidp
      double complex proj(nw,nb)
      double precision r0
      double precision rmatrx(nw*nb,nw*nb)
      double precision sig2
      double precision sigma
      double complex wall(nw,nb)
      double precision xmax(nb)
      double precision ylimit
      double precision ymax
      double precision ymin
      double complex z1
      double complex zcr(nw,nb)
      double complex zi
      double complex zz

      pi = 3.141592653589793D+00
      period = 3.0D+00
c
c  Compute arclength.
c
      matdim = 0
      arcl = 0.0D+00
      ymin = 1.0D+30
      ymax = -1.0D+30
      pidp = pi / period

      do l = 1, nb
        matdim = matdim + nwall(l)
        do k = 1, nwall(l)
          arcl = arcl + cdabs ( wall(k,l) - wall(1+mod(k,nwall(l)),l) )
        end do
      end do
c
c  Compute core radius.
c
      r0 = 0.5D+00 * arcl / dble ( matdim )
      sigma = r0 / 2.0D+00
c
c  Define creation points.
c
      do l = 1, nb

        do k = 1, nwall(l)
          zz = wall(1+mod(k+nwall(l)-2,nwall(l)),l)
     &      - wall(1+mod(k,nwall(l)),l)
          zcr(k,l) = wall(k,l) 
     &      + dcmplx ( 0.0D+00, r0 / cdabs ( zz ) ) * zz
        end do
c
c  Check that wall and creation points are not crossed due to
c  too sharp a concave kink or an error in defining the body.
c  Also find highest, lowest and right-most point.
c
        xmax(l) = dreal ( zcr(1,l) )
        ls = 0

        do k = 1, nwall(l)

          ymin = min ( ymin, dimag ( zcr(k,l) ) )
          ymax = max ( ymax, dimag ( zcr(k,l) ) )
          xmax(l) = max ( xmax(l), dreal ( zcr(k,l) ) )
          kp = 1 + mod ( k, nwall(l) )

          if ( 0.0D+00 .lt. dreal ( ( zcr(kp,l) - zcr(k,l) ) *
     &        dconjg ( wall(kp,l) - wall(k,l) ) ) ) then
            ls = l
            ks = k
          end if

        end do

      end do
c
c  The "main period" will be between ylimit and ylimit + period.
c
      ylimit = ( ymin - period + ymax ) / 2.0D+00
c
c  Project creation points into main period.  This is technical.
c
      do l = 1, nb
        do k = 1, nwall(l)
          proj(k,l) = zcr(k,l) 
     &      - dcmplx ( 0.0D+00, period *
     &      ( int ( 5.0D+00 + ( dimag ( zcr(k,l) ) - ylimit ) / period ) 
     &      - 5.0D+00 ) )
        end do
      end do
c
c  Compute matrix.
c
      sig2 = ( 2.0D+00 * pidp * sigma ) ** 2
      i0 = 0

      do l1 = 1, nb

        j0 = 0

        do l2 = 1, nb

          if ( l1 .eq. l2 ) then
            kron = 1
          else
            kron = 0
          end if

          do j = 1, nwall(l2)
            rmatrx(i0+1,j0+j) = kron
            z1 = cdexp ( ( wall(1,l1) - zcr(j,l2) ) * pidp )
            z1 = z1 - 1.0D+00 / z1
            dum = sig2 + dreal ( z1 )**2 + dimag ( z1 )**2
            do i = 2, nwall(l1)
              zi = cdexp ( ( wall(i,l1) - zcr(j,l2) ) * pidp )
              zz = zi - 1.0D+00 / zi
              rmatrx(i0+i,j0+j) = -0.25D+00 / pi * dlog ( dum /
     &          ( sig2 + dreal ( zz ) ** 2 + dimag ( zz ) ** 2 ) )
            end do
          end do

          j0 = j0 + nwall(l2)

        end do

        i0 = i0 + nwall(l1)

      end do
c
c  Gauss elimination.
c
      do i = 1, matdim
        rmatrx(i,i) = 1.0D+00 / rmatrx(i,i)
        do j = i + 1, matdim
          rmatrx(j,i) = rmatrx(j,i) * rmatrx(i,i)
          do k = i + 1, matdim
            rmatrx(j,k) = rmatrx(j,k) - rmatrx(j,i) * rmatrx(i,k)
          end do
        end do
      end do

      return
      end
      subroutine emitst ( er, fp, tm )

c*********************************************************************72
c
cc EMITST is the test program for EMIT.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer nb
      integer nv
      integer nvm
      integer nw

      parameter ( nb = 5 )
      parameter ( nv = 1000 )
      parameter ( nvm = 1500 )
      parameter ( nw = 100 )

      double precision ans
      double precision cp(nw,nb)
      double precision cptime
      double precision dpds(nw,nb)
      double precision er
      double complex expmz(nvm)
      double complex expz(nvm)
      double precision f7
      double complex force(nb)
      double precision fp
      double precision gamma(nvm)
      integer i
      integer it
      integer j
      integer nwall(nb)
      double precision ps(nvm)
      double precision psi(nw)
      double complex refpt(nb)
      double precision rhs(nw*nb)
      double precision rmatrx(nw*nb,nw*nb)
      double precision rmom(nb)
      double precision t1
      double precision t2
      double precision t30
      double precision tm
      double complex wall(nw,nb)
      double complex z(nvm)
      double complex zcr(nw,nb)

      it = 10
      ans = 6.0088546832072D+00
!
!  Random initialization.
!
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t2 = f7 / t30

      do j = 1, nb
        nwall(j) = nw
        refpt(j) = 0.0D+00
        force(j) = 0.0D+00
        rmom(j) = 0.0D+00
        do i = 1, nw
          t1 = dmod ( f7 * t2, 1.0D+00 )
          t2 = dmod ( f7 * t1, 1.0D+00 )
          wall(i,j) = dcmplx ( t1, t2 )
          t1 = dmod ( f7 * t2, 1.0D+00 )
          t2 = dmod ( f7 * t1, 1.0D+00 )
          zcr(i,j) = dcmplx ( t1, t2 )
          dpds(i,j) = 0.0D+00
        end do
      end do

      do j = 1, nw * nb
        rmatrx(j,j) = 1.0D+00
        do i = 1, j - 1
          t2 = dmod ( f7 * t2, 1.0D+00 )
          rmatrx(i,j) = 0.001D+00 * t2
          rmatrx(j,i) = 0.001D+00 * t2
        end do
      end do

      do i = 1, nvm
        t1 = dmod ( f7 * t2, 1.0D+00 )
        t2 = dmod ( f7 * t1, 1.0D+00 )
        z(i) = dcmplx ( t1, t2 )
        t2 = dmod ( f7 * t2, 1.0D+00 )
        gamma(i) = t2
      end do
c
c  Timing.
c
      tm = cptime ( )

      do i = 1, it
        call emit ( nb, nvm, nw, cp, dpds, expmz, expz, force, 
     &    gamma, nwall, ps, psi, refpt, rhs, rmatrx, rmom, wall, 
     &    z, zcr )
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( rhs(19) - ans ) / ans )
      fp = dble ( it * ( 56 * nv + nb * nw 
     &  * ( 97 + 44 * nv + 2 * nb * nw ) ) )

      return
      end
      subroutine emit ( nb, nvm, nw, cp, dpds, expmz, expz, force, 
     &  gamma, nwall, ps, psi, refpt, rhs, rmatrx, rmom, wall, z, zcr )

c*********************************************************************72
c
cc EMIT creates new vortices according to certain boundary conditions.
c
c  Discussion:
c
c    This function was extracted from a vortex method program.
c    It emits new vortices to satisfy the boundary condition.
c    It also finishes computing pressure, forces, and other quantities.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer nb
      integer nvm
      integer nw

      double precision chord
      double precision cp(nw,nb)
      double precision cpm
      double precision cupst
      double precision delt
      double precision dpds(nw,nb)
      double complex dum3
      double complex expmwk
      double complex expmz(nvm)
      double complex expwkl
      double complex expz(nvm)
      double complex force(nb)
      double precision gamma(nvm)
      integer i
      integer i0
      integer j
      integer k
      integer l
      integer matdim
      integer nolld
      integer nv
      integer nwall(nb)
      double precision period
      double precision pi
      double precision pidp
      double precision ps(nvm)
      double precision psi(nw)
      double complex refpt(nb)
      double precision rhs(nw*nb)
      double precision rmatrx(nw*nb,nw*nb)
      double precision rmom(nb)
      double precision sig2
      double precision sps
      double precision u0
      double complex uupstr
      double complex wall(nw,nb)
      double complex z(nvm)
      double complex zcr(nw,nb)

      period = 3.0D+00
      sig2 = 3.0D+00
      u0 = 4.0D+00
      matdim = 500 
      delt = 1.0D+00
      chord = 5.0D+00 
      pi = 3.141592653589793D+00
      uupstr = dcmplx ( 3.0D+00, 4.0D+00 )
c
c  Store exp(z(i)) and exp(-z(i)) to reduce work in inner loop.
c
c  Note that the NV used here is a variable, whereas the NV in the
c  calling program is a constant.  They are separate quantities.
c
      nv = 1000
      pidp = pi / period

      do i = 1, nv
        expz(i) = cdexp ( z(i) * pidp )
        expmz(i) = 1.0D+00 / expz(i)
      end do

      i0 = 0
      cupst = ( dreal ( uupstr ) ) ** 2 + ( dimag ( uupstr ) ) ** 2

      do l = 1, nb

        do k = 1, nwall(l)

          expwkl = cdexp ( wall(k,l) * pidp )
          expmwk = 1.0D+00 / expwkl

          sps = 0.0D+00
          do i = 1, nv
            dum3 = expz(i) * expmwk - expwkl * expmz(i)
            ps(i) = gamma(i) * dlog ( dreal ( dum3 ) ** 2 +
     &        dimag ( dum3 ) ** 2 + sig2 )
            sps = sps + ps(i)
          end do

          psi(k) = dimag ( wall(k,l) 
     &      * dconjg ( uupstr + dcmplx ( 0.0D+00, u0 ) ) )
     &      - sps * 0.25D+00 / pi
        end do
c
c  Compute the right-hand side.
c
        do k = 1, nwall(l)
          rhs(i0+k) = psi(k) - psi(1)
        end do

        i0 = i0 + nwall(l)

      end do
c
c  Solve the system.
c
      do i = 1, matdim
        do j = i + 1, matdim
          rhs(j) = rhs(j) - rmatrx(j,i) * rhs(i)
        end do
      end do

      do i = matdim, 1, -1
        rhs(i) = rmatrx(i,i) * rhs(i)
        do j = 1, i - 1
          rhs(j) = rhs(j) - rmatrx(j,i) * rhs(i)
        end do
      end do
c
c  Create new vortices.
c
      nolld = nv
      i0 = 0

      do l = 1, nb

        do k = 1, nwall(l)
c
c  Put the new vortex at the end of the array.
c
          nv = nv + 1
          z(nv) = zcr(k,l)
          gamma(nv) = rhs(i0+k)
c
c  Record the gain of linear and angular momentum.
c
          force(l) = force(l) + gamma(nv) * z(nv)
          rmom(l) = rmom(l) 
     &      + gamma(nv) * ( dreal ( z(nv) - refpt(l) ) ** 2
     &      + dimag ( z(nv) - refpt(l) ) ** 2 )
          dpds(k,l) = dpds(k,l) - gamma(nv)

        end do
c
c  Filter and integrate pressure gradient to get pressure.
c
        cp(1,l) = 0.0D+00
        cpm = -1.0D+30

        do k = 2, nwall(l)
          cp(k,l) = cp(k-1,l) + ( 3.0D+00 * ( dpds(k,l) + dpds(k-1,l) )
     &      + dpds(1+mod(k,nwall(l)),l)
     &      + dpds(1+mod(k+nwall(l)-3,nwall(l)),l) )
     &      / ( 4.0D+00 * delt * cupst )
          cpm = max ( cpm, cp(k,l) )
        end do
c
c  Normalize the pressure.
c
        do k = 1, nwall(l)
          cp(k,l) = cp(k,l) - cpm
        end do
c
c  Finish computing force and moment, as time rate of change of linear
c  and angular momentum.
c
        force(l) = force(l) 
     &    * dcmplx ( 0.0D+00, 2.0D+00 / ( delt * chord * cupst ) )

        rmom(l) = rmom(l) * 2.0D+00 / ( delt * chord ** 2 * cupst )

        i0 = i0 + nwall(l)

      end do

      return
      end
      subroutine vpetst ( er, fp, tm )

c*********************************************************************72
c
cc VPETST is the test program for VPENTA.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer jl
      integer ju
      integer kl
      integer ku
      integer nja
      integer njb

      parameter ( jl = 1 )
      parameter ( ju = 128 )
      parameter ( kl = 1 )
      parameter ( ku = 128 )
      parameter ( nja = 128 )
      parameter ( njb = 128 )
      
      double precision a(nja,njb)
      double precision ans
      double precision b(nja,njb)
      double precision c(nja,njb)
      double precision cptime
      double precision d(nja,njb)
      double precision e(nja,njb)
      double precision er
      double precision f(nja,njb,3)
      double precision f7
      double precision fp
      double precision fx(nja,njb,3)
      integer i
      integer it
      integer j
      integer k
      integer lf
      double precision t
      double precision t30
      double precision tm
      double precision x(nja,njb)
      double precision y(nja,njb)

      it = 400
      ans = -0.354649411858726D+00
      lf = nja * njb * 3
c
c  Random initialization.
c
      f7 = 78125.0D+00
      t30 = 1073741824.0D+00
      t = f7 / t30

      do j = kl, ku
        do i = jl, ju
          t = dmod ( f7 * t, 1.0D+00 )
          a(i,j) = t
          t = dmod ( f7 * t, 1.0D+00 )
          b(i,j) = t
          t = dmod ( f7 * t, 1.0D+00 )
          c(i,j) = t
          t = dmod ( f7 * t, 1.0D+00 )
          d(i,j) = t
          t = dmod ( f7 * t, 1.0D+00 )
          e(i,j) = t
          do k = 1, 3
            t = dmod ( f7 * t, 1.0D+00 )
            fx(i,j,k) = t
          end do
        end do
      end do
c
c  Timing.
c
      tm = cptime ( )

      do i = 1, it
        call copy ( lf, fx, f )
        call vpenta ( jl, ju, kl, ku, nja, njb, a, b, c, d, e, f,
     &    x, y )
      end do

      tm = cptime ( ) - tm
c
c  Results.
c
      er = dabs ( ( f(19,19,1) - ans ) / ans )
      fp = dble ( it * ku * ( 40 * ku - 53 ) )

      return
      end
      subroutine vpenta ( jl, ju, kl, ku, nja, njb, a, b, c, d, e, f,
     &  x, y )

c*********************************************************************72
c
cc VPENTA inverts 3 pentadiagonal systems simultaneously.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 November 2010
c
c  Author:
c
c    Original FORTRAN77 version by David Bailey.
c    This version by John Burkardt.
c
      implicit none

      integer jl
      integer ju
      integer kl
      integer ku
      integer nja
      integer njb

      double precision a(nja,njb)
      double precision b(nja,njb)
      double precision c(nja,njb)
      double precision d(nja,njb)
      double precision e(nja,njb)
      double precision f(nja,njb,3)
      integer j
      integer jx
      integer k
      double precision rld
      double precision rld1
      double precision rld2
      double precision rldi
      double precision x(nja,njb)
      double precision y(nja,njb)
c
c  Start forward generation process and sweep.
c
      j = jl
      do k = kl, ku
        rld = c(j,k)
        rldi = 1.0D+00 / rld
        f(j,k,1) = f(j,k,1) * rldi
        f(j,k,2) = f(j,k,2) * rldi
        f(j,k,3) = f(j,k,3) * rldi
        x(j,k) = d(j,k) * rldi
        y(j,k) = e(j,k) * rldi
      end do

      j = jl + 1
      do k = kl, ku
        rld1 = b(j,k)
        rld = c(j,k) - rld1 * x(j-1,k)
        rldi = 1.0D+00 / rld
        f(j,k,1) = ( f(j,k,1) - rld1 * f(j-1,k,1) ) * rldi
        f(j,k,2) = ( f(j,k,2) - rld1 * f(j-1,k,2) ) * rldi
        f(j,k,3) = ( f(j,k,3) - rld1 * f(j-1,k,3) ) * rldi
        x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi
        y(j,k) = e(j,k) * rldi
      end do

      do j = jl + 2, ju - 2
        do k = kl, ku
          rld2 = a(j,k)
          rld1 = b(j,k) - rld2 * x(j-2,k)
          rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) )
          rldi = 1.0D+00 / rld
          f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) 
     &      - rld1 * f(j-1,k,1) ) * rldi
          f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) 
     &      - rld1 * f(j-1,k,2) ) * rldi
          f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) 
     &      - rld1 * f(j-1,k,3) ) * rldi
          x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi
          y(j,k) = e(j,k) * rldi
        end do
      end do

      j = ju - 1
      do k = kl, ku
        rld2 = a(j,k)
        rld1 = b(j,k) - rld2 * x(j-2,k)
        rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) )
        rldi = 1.0D+00 / rld
        f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) - rld1 * f(j-1,k,1) ) 
     &    * rldi
        f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) - rld1 * f(j-1,k,2) ) 
     &    * rldi
        f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) - rld1 * f(j-1,k,3) ) 
     &    * rldi
        x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi
      end do

      j = ju
      do k = kl, ku
        rld2 = a(j,k)
        rld1 = b(j,k) - rld2 * x(j-2,k)
        rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) )
        rldi = 1.0D+00 / rld
        f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) - rld1 * f(j-1,k,1) ) 
     &    * rldi
        f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) - rld1 * f(j-1,k,2) ) 
     &    * rldi
        f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) - rld1 * f(j-1,k,3) ) 
     &    * rldi
      end do
c
c  Back sweep solution.
c
      do k = kl, ku
        f(ju,k,1) = f(ju,k,1)
        f(ju,k,2) = f(ju,k,2)
        f(ju,k,3) = f(ju,k,3)
        f(ju-1,k,1) = f(ju-1,k,1) - x(ju-1,k) * f(ju,k,1)
        f(ju-1,k,2) = f(ju-1,k,2) - x(ju-1,k) * f(ju,k,2)
        f(ju-1,k,3) = f(ju-1,k,3) - x(ju-1,k) * f(ju,k,3)
      end do

      do j = 2, ju - jl
        jx = ju - j
        do k = kl, ku
          f(jx,k,1) = f(jx,k,1) - x(jx,k) * f(jx+1,k,1) -
     &      y(jx,k) * f(jx+2,k,1)
          f(jx,k,2) = f(jx,k,2) - x(jx,k) * f(jx+1,k,2) -
     &      y(jx,k) * f(jx+2,k,2)
          f(jx,k,3) = f(jx,k,3) - x(jx,k) * f(jx+1,k,3) -
     &      y(jx,k) * f(jx+2,k,3)
        end do
      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end