program main

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

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) er
  real ( kind = rk ) fp
  integer i
  character ( len = 8 ) :: pn(8) = (/ &
    'BTRIX   ', &
    'CFFT2D  ', &
    'CHOLSKY ', &
    'EMIT    ', &
    'GMTRY   ', &
    'MXM     ', &
    'VPENTA  ', &
    'Total   ' /)
  real ( kind = rk ) rt
  real ( kind = rk ) te
  real ( kind = rk ) tf
  real ( kind = rk ) tm
  real ( kind = rk ) tt

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

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'nas():'
  write ( *, '(a)' ) '  FORTRAN90 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 == 1 ) then
      call btrix_test ( er, fp, tm )
    else if ( i == 2 ) then
      call cfft2d_test ( er, fp, tm )
    else if ( i == 3 ) then
      call cholsky_test ( er, fp, tm )
    else if ( i == 4 ) then
      call emit_test ( er, fp, tm )
    else if ( i == 5 ) then
      call gmtry_test ( er, fp, tm )
    else if ( i == 6 ) then
      call mxm_test ( er, fp, tm )
    else if ( i == 7 ) then
      call vpenta_test ( er, fp, tm )
    end if

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

    te = te + er
    tf = tf + fp
    tt = tt + tm

  end do

  er = te
  fp = tf
  tm = tt
  rt = 1.0D-06 * tf / tt
  i = 8
  write ( *, '(a)' ) ' '
  write ( *, '(1x,a8,1p2e15.4,0pf12.4,f12.2)' ) &
    pn(i), er, fp, tm, rt
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'nas():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine btrix_test ( er, fp, tm )

!*****************************************************************************80
!
!! BTRIX_TEST tests BTRIX.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer jd
  integer kd
  integer ld
  integer md

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

  real ( kind = rk ) a(5,5,md,md)
  real ( kind = rk ) ans
  real ( kind = rk ) b(5,5,md,md)
  real ( kind = rk ) bx(5,5,md,md)
  real ( kind = rk ) c(5,5,md,md)
  real ( kind = rk ) er
  real ( kind = rk ) f7
  real ( kind = rk ) 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
  real ( kind = rk ) s(jd,kd,ld,5)
  real ( kind = rk ) sx(jd,kd,ld,5)
  real ( kind = rk ) t
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm

  js = 2 
  je = 29
  ls = 2 
  le = 29
  it = 20
  ans = -0.286282658663962D+00
  nb = 25 * md * md
  ns = jd * kd * ld * 5
!
!  Random initialization.
!
  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 = mod ( f7 * t, 1.0D+00 )
          a(i,j,k,l) = t
          t = mod ( f7 * t, 1.0D+00 )
          bx(i,j,k,l) = t
          t = mod ( 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 = mod ( f7 * t, 1.0D+00 )
          sx(i,j,k,l) = t
        end do
      end do
    end do
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

  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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( 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 )

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

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer jd
  integer kd
  integer ld
  integer md

  real ( kind = rk ) a(5,5,md,md)
  real ( kind = rk ) b(5,5,md,md)
  real ( kind = rk ) c(5,5,md,md)
  real ( kind = rk ) c1
  real ( kind = rk ) c2
  real ( kind = rk ) c3
  real ( kind = rk ) c4
  real ( kind = rk ) c5
  real ( kind = rk ) d1
  real ( kind = rk ) d2
  real ( kind = rk ) d3
  real ( kind = rk ) d4
  real ( kind = rk ) d5
  integer j
  integer je
  integer jem1
  integer js
  integer k
  integer l
  integer le
  integer ls
  real ( kind = rk ) l11(md)
  real ( kind = rk ) l21(md)
  real ( kind = rk ) l31(md)
  real ( kind = rk ) l41(md)
  real ( kind = rk ) l51(md)
  real ( kind = rk ) l22(md)
  real ( kind = rk ) l32(md)
  real ( kind = rk ) l33(md)
  real ( kind = rk ) l42(md)
  real ( kind = rk ) l43(md)
  real ( kind = rk ) l44(md)
  real ( kind = rk ) l52(md)
  real ( kind = rk ) l53(md)
  real ( kind = rk ) l54(md)
  real ( kind = rk ) l55(md)
  integer m
  integer n
  real ( kind = rk ) s(jd,kd,ld,5)
  real ( kind = rk ) u12(md)
  real ( kind = rk ) u13(md)
  real ( kind = rk ) u14(md)
  real ( kind = rk ) u15(md)
  real ( kind = rk ) u23(md)
  real ( kind = rk ) u24(md)
  real ( kind = rk ) u25(md)
  real ( kind = rk ) u34(md)
  real ( kind = rk ) u35(md)
  real ( kind = rk ) u45(md)
!
!  Part 1.  Forward block sweep.
!
  do j = js, je
!
!  Step 1.  Construct L(I) in B.
!
    if ( j /= 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
!
!  Step 2.  Compute L inverse.
!
!  A.  Decompose L(I) into L and U.
!
    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
!
!  Step 3.  Solve for intermediate vector.
!
!  A.  Construct the right hand side.
!
    if ( j /= 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
!
!  B. Intermediate vector.
!
!  Forward substitution.
!
    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)
!
!  Backward substitution.
!
      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
!
!  Step 4.  Construct U(I) = inverse(L(I))*C(I+1) by columns and store in B.
!
    if ( j /= je ) then

      do n = 1, 5
        do l = ls, le
!
!  Forward substitution.
!
          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)
!
!  Backward substitution.
!
          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
!
!  Part 2.  Backward block sweep.
!
  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 copy ( n, a, b )

!*****************************************************************************80
!
!! COPY copies a real ( kind = rk ) vector.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a(n)
  real ( kind = rk ) b(n)
  integer i

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

  return
end
subroutine mxm_test ( er, fp, tm )

!*****************************************************************************80
!
!! MXM_TEST tests MXM.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer l
  integer m
  integer n

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

  real ( kind = rk ) a(l,m)
  real ( kind = rk ) ans
  real ( kind = rk ) b(m,n)
  real ( kind = rk ) c(l,n)
  real ( kind = rk ) er
  real ( kind = rk ) f7
  real ( kind = rk ) fp
  integer i
  integer ii
  integer it
  integer j
  real ( kind = rk ) t
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm

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

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

  do j = 1, n
    do i = 1, m
      t = mod ( f7 * t, 1.0D+00 )
      b(i,j) = t
    end do
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( c(19,19) - ans ) / ans )
  fp = dble ( 2 * it * l * m * n )

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

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

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer l
  integer m
  integer n

  real ( kind = rk ) a(l,m)
  real ( kind = rk ) b(m,n)
  real ( kind = rk ) 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 cfft2d_test ( er, fp, tm )

!*****************************************************************************80
!
!! CFFT2D_TEST is the test program for CFFT2D1 and CFFTD2.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 July 2015
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer m1
  integer n

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

  real ( kind = rk ) ans
  real ( kind = rk ) er
  real ( kind = rk ) f7
  real ( kind = rk ) fp
  integer i
  integer ip(2*max(m,n))
  integer it
  integer j
  integer k
  real ( kind = rk ) rmn
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm
  complex ( kind = ck ) w1(m)
  complex ( kind = ck ) w2(n)
  complex ( kind = ck ) x(m1,n)

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

  do j = 1, n
    do i = 1, m
      t1 = mod ( f7 * t2, 1.0D+00 )
      t2 = mod ( 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 )
!
!  Timing.
!
  call cpu_time ( time1 )

  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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( dble ( x(19,19) ) - ans ) / ans )
  fp = dble ( it * m * n ) * ( 2.0D+00  &
    + 10.0D+00 * log ( dble ( m * n ) ) / log ( 2.0D+00 ) )

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

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

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer m1
  integer n

  complex ( kind = ck ) ct
  complex ( kind = ck ) 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
  real ( kind = rk ) pi
  real ( kind = rk ) t
  complex ( kind = ck ) w(m)
  complex ( kind = ck ) x(m1,n)

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

  l = 1
  i1 = 1

  do

    i2 = 3 - i1

    do j = l, m2, l

      cx = w(j-l+1)
      if ( is < 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 ( m2 < l ) then
      exit
    end if

  end do

  do i = 1, m
    ii = ip(i1,i)
    if ( i < 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 )

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

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m1
  integer n

  complex ( kind = ck ) ct
  complex ( kind = ck ) 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
  real ( kind = rk ) pi
  real ( kind = rk ) t
  complex ( kind = ck ) w(n)
  complex ( kind = ck ) x(m1,n)

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

  l = 1
  i1 = 1

  do

    i2 = 3 - i1

    do j = l, n2, l

      cx = w(j-l+1)
      if ( is < 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 ( n2 < l ) then
      exit
    end if

  end do

  do i = 1, n
    ii = ip(i1,i)
    if ( i < 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 cholsky_test ( er, fp, tm )

!*****************************************************************************80
!
!! CHOLSKY_TEST tests CHOLSKY.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  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 )

  real ( kind = rk ) a(0:ida,-m:0,0:n)
  real ( kind = rk ) ans
  real ( kind = rk ) ax(0:ida,-m:0,0:n)
  real ( kind = rk ) b(0:nrhs,0:nmat,0:n)
  real ( kind = rk ) bx(0:nrhs,0:nmat,0:n)
  real ( kind = rk ) er
  real ( kind = rk ) f7
  real ( kind = rk ) fp
  integer i
  integer it
  integer j
  integer k
  integer la
  integer lb
  real ( kind = rk ) t
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm

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

  do k = 0, n
    do j = - m, 0
      do i = 0, ida
        t = mod ( 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 = mod ( f7 * t, 1.0D+00 )
        bx(i,j,k) = t
      end do
    end do
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

  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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( b(1,19,19) - ans ) / ans )
  fp = dble ( it * ( nmat + 1 ) * 4403 )

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

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

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer ida
  integer idb
  integer m
  integer n
  integer nrhs

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

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

    i0 = max ( -m, -j )
!
!  Off diagonal elements.
!
    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
!
!  Store inverse of diagonal elements.
!
    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 / sqrt ( abs ( epss(l) + a(l,0,j) ) )
    end do

  end do
!
!  Solution.
!
  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 gmtry_test ( er, fp, tm )

!*****************************************************************************80
!
!! GMTRY_TEST tests GMTRY.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nb
  integer nw

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

  real ( kind = rk ) ans
  real ( kind = rk ) er
  real ( kind = rk ) f7
  real ( kind = rk ) fp
  integer i
  integer it
  integer j
  integer lw
  integer nwall(nb)
  complex ( kind = ck ) proj(nw,nb)
  real ( kind = rk ) rmatrx(nw*nb,nw*nb)
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm
  complex ( kind = ck ) wall(nw,nb)
  real ( kind = rk ) xmax(nb)
  complex ( kind = ck ) zcr(nw,nb)

  it = 2
  ans = -2.57754233214174D+00
  lw = 2 * nw * nb
!
!  Random initialization.
!
  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 = mod ( f7 * t2, 1.0D+00 )
      t2 = mod ( f7 * t1, 1.0D+00 )
      wall(i,j) = dcmplx ( t1, t2 )
    end do
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( 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 )

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

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nb
  integer nw

  real ( kind = rk ) arcl
  real ( kind = rk ) 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)
  real ( kind = rk ) period
  real ( kind = rk ) pi
  real ( kind = rk ) pidp
  complex ( kind = ck ) proj(nw,nb)
  real ( kind = rk ) r0
  real ( kind = rk ) rmatrx(nw*nb,nw*nb)
  real ( kind = rk ) sig2
  real ( kind = rk ) sigma
  complex ( kind = ck ) wall(nw,nb)
  real ( kind = rk ) xmax(nb)
  real ( kind = rk ) ylimit
  real ( kind = rk ) ymax
  real ( kind = rk ) ymin
  complex ( kind = ck ) z1
  complex ( kind = ck ) zcr(nw,nb)
  complex ( kind = ck ) zi
  complex ( kind = ck ) zz

  pi = 3.141592653589793D+00
  period = 3.0D+00
!
!  Compute arclength.
!
  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 + abs ( wall(k,l) - wall(1+mod(k,nwall(l)),l) )
    end do
  end do
!
!  Compute core radius.
!
  r0 = 0.5D+00 * arcl / dble ( matdim )
  sigma = r0 / 2.0D+00
!
!  Define creation points.
!
  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 / abs ( zz ) ) * zz
    end do
!
!  Check that wall and creation points are not crossed due to
!  too sharp a concave kink or an error in defining the body.
!  Also find highest, lowest and right-most point.
!
    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 < 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
!
!  The "main period" will be between ylimit and ylimit + period.
!
  ylimit = ( ymin - period + ymax ) / 2.0D+00
!
!  Project creation points into main period.  This is technical.
!
  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
!
!  Compute matrix.
!
  sig2 = ( 2.0D+00 * pidp * sigma ) ** 2
  i0 = 0

  do l1 = 1, nb

    j0 = 0

    do l2 = 1, nb

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

      do j = 1, nwall(l2)
        rmatrx(i0+1,j0+j) = kron
        z1 = exp ( ( 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 = exp ( ( wall(i,l1) - zcr(j,l2) ) * pidp )
          zz = zi - 1.0D+00 / zi
          rmatrx(i0+i,j0+j) = -0.25D+00 / pi * log ( dum / &
            ( sig2 + dreal ( zz ) ** 2 + dimag ( zz ) ** 2 ) )
        end do
      end do

      j0 = j0 + nwall(l2)

    end do

    i0 = i0 + nwall(l1)

  end do
!
!  Gauss elimination.
!
  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 emit_test ( er, fp, tm )

!*****************************************************************************80
!
!! EMIT_TEST tests EMIT.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nb
  integer nv
  integer nvm
  integer nw

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

  real ( kind = rk ) ans
  real ( kind = rk ) cp(nw,nb)
  real ( kind = rk ) dpds(nw,nb)
  real ( kind = rk ) er
  complex ( kind = ck ) expmz(nvm)
  complex ( kind = ck ) expz(nvm)
  real ( kind = rk ) f7
  complex ( kind = ck ) force(nb)
  real ( kind = rk ) fp
  real ( kind = rk ) gamma(nvm)
  integer i
  integer it
  integer j
  integer nwall(nb)
  real ( kind = rk ) ps(nvm)
  real ( kind = rk ) psi(nw)
  complex ( kind = ck ) refpt(nb)
  real ( kind = rk ) rhs(nw*nb)
  real ( kind = rk ) rmatrx(nw*nb,nw*nb)
  real ( kind = rk ) rmom(nb)
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm
  complex ( kind = ck ) wall(nw,nb)
  complex ( kind = ck ) z(nvm)
  complex ( kind = ck ) 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 = mod ( f7 * t2, 1.0D+00 )
      t2 = mod ( f7 * t1, 1.0D+00 )
      wall(i,j) = dcmplx ( t1, t2 )
      t1 = mod ( f7 * t2, 1.0D+00 )
      t2 = mod ( 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 = mod ( 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 = mod ( f7 * t2, 1.0D+00 )
    t2 = mod ( f7 * t1, 1.0D+00 )
    z(i) = dcmplx ( t1, t2 )
    t2 = mod ( f7 * t2, 1.0D+00 )
    gamma(i) = t2
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

  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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( 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 )

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

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )
  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nb
  integer nvm
  integer nw

  real ( kind = rk ) chord
  real ( kind = rk ) cp(nw,nb)
  real ( kind = rk ) cpm
  real ( kind = rk ) cupst
  real ( kind = rk ) delt
  real ( kind = rk ) dpds(nw,nb)
  complex ( kind = ck ) dum3
  complex ( kind = ck ) expmwk
  complex ( kind = ck ) expmz(nvm)
  complex ( kind = ck ) expwkl
  complex ( kind = ck ) expz(nvm)
  complex ( kind = ck ) force(nb)
  real ( kind = rk ) gamma(nvm)
  integer i
  integer i0
  integer j
  integer k
  integer l
  integer matdim
  integer nolld
  integer nv
  integer nwall(nb)
  real ( kind = rk ) period
  real ( kind = rk ) pi
  real ( kind = rk ) pidp
  real ( kind = rk ) ps(nvm)
  real ( kind = rk ) psi(nw)
  complex ( kind = ck ) refpt(nb)
  real ( kind = rk ) rhs(nw*nb)
  real ( kind = rk ) rmatrx(nw*nb,nw*nb)
  real ( kind = rk ) rmom(nb)
  real ( kind = rk ) sig2
  real ( kind = rk ) sps
  real ( kind = rk ) u0
  complex ( kind = ck ) uupstr
  complex ( kind = ck ) wall(nw,nb)
  complex ( kind = ck ) z(nvm)
  complex ( kind = ck ) 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 )
!
!  Store exp(z(i)) and exp(-z(i)) to reduce work in inner loop.
!
!  Note that the NV used here is a variable, whereas the NV in the
!  calling program is a constant.  They are separate quantities.
!
  nv = 1000
  pidp = pi / period

  do i = 1, nv
    expz(i) = exp ( 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 = exp ( 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) * log ( 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
!
!  Compute the right-hand side.
!
    do k = 1, nwall(l)
      rhs(i0+k) = psi(k) - psi(1)
    end do

    i0 = i0 + nwall(l)

  end do
!
!  Solve the system.
!
  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
!
!  Create new vortices.
!
  nolld = nv
  i0 = 0

  do l = 1, nb

    do k = 1, nwall(l)
!
!  Put the new vortex at the end of the array.
!
      nv = nv + 1
      z(nv) = zcr(k,l)
      gamma(nv) = rhs(i0+k)
!
!  Record the gain of linear and angular momentum.
!
      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
!
!  Filter and integrate pressure gradient to get pressure.
!
    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
!
!  Normalize the pressure.
!
    do k = 1, nwall(l)
      cp(k,l) = cp(k,l) - cpm
    end do
!
!  Finish computing force and moment, as time rate of change of linear
!  and angular momentum.
!
    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 vpenta_test ( er, fp, tm )

!*****************************************************************************80
!
!! VPENTA_TEST tests VPENTA.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 November 2010
!
!  Author:
!
!    Original FORTRAN77 version by David Bailey.
!    FORTRAN90 version by John Burkardt.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  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 )
  
  real ( kind = rk ) a(nja,njb)
  real ( kind = rk ) ans
  real ( kind = rk ) b(nja,njb)
  real ( kind = rk ) c(nja,njb)
  real ( kind = rk ) d(nja,njb)
  real ( kind = rk ) e(nja,njb)
  real ( kind = rk ) er
  real ( kind = rk ) f(nja,njb,3)
  real ( kind = rk ) f7
  real ( kind = rk ) fp
  real ( kind = rk ) fx(nja,njb,3)
  integer i
  integer it
  integer j
  integer k
  integer lf
  real ( kind = rk ) t
  real ( kind = rk ) t30
  real ( kind = rk ) time1
  real ( kind = rk ) time2
  real ( kind = rk ) tm
  real ( kind = rk ) x(nja,njb)
  real ( kind = rk ) y(nja,njb)

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

  do j = kl, ku
    do i = jl, ju
      t = mod ( f7 * t, 1.0D+00 )
      a(i,j) = t
      t = mod ( f7 * t, 1.0D+00 )
      b(i,j) = t
      t = mod ( f7 * t, 1.0D+00 )
      c(i,j) = t
      t = mod ( f7 * t, 1.0D+00 )
      d(i,j) = t
      t = mod ( f7 * t, 1.0D+00 )
      e(i,j) = t
      do k = 1, 3
        t = mod ( f7 * t, 1.0D+00 )
        fx(i,j,k) = t
      end do
    end do
  end do
!
!  Timing.
!
  call cpu_time ( time1 )

  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

  call cpu_time ( time2 )

  tm = time2 - time1
!
!  Results.
!
  er = abs ( ( 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 )

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

  integer, parameter :: rk = kind ( 1.0D+00 )

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

  real ( kind = rk ) a(nja,njb)
  real ( kind = rk ) b(nja,njb)
  real ( kind = rk ) c(nja,njb)
  real ( kind = rk ) d(nja,njb)
  real ( kind = rk ) e(nja,njb)
  real ( kind = rk ) f(nja,njb,3)
  integer j
  integer jx
  integer k
  real ( kind = rk ) rld
  real ( kind = rk ) rld1
  real ( kind = rk ) rld2
  real ( kind = rk ) rldi
  real ( kind = rk ) x(nja,njb)
  real ( kind = rk ) y(nja,njb)
!
!  Start forward generation process and sweep.
!
  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
!
!  Back sweep solution.
!
  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 ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

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

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

  return
end