program main

!*****************************************************************************80
!
!! wathen_matrix_test() tests wathen_matrix().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'wathen_matrix_test():'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  test wathen_matrix().'

  call test01 ( )
  call test02 ( )
  call test05 ( )
  call test06 ( )
  call test07 ( )
  call test08 ( )
  call test10 ( )
  call test11 ( )
  call test115 ( )
  call wathen_xy_test ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'wathen_matrix_test():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine test01 ( )

!*****************************************************************************80
!
!! test01() assembles, factor and solve using wathen_ge().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    06 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer i
  integer info
  integer, allocatable :: ipvt(:)
  integer job
  integer n
  integer nx
  integer ny
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test01():'
  write ( *, '(a)' ) '  Assemble, factor and solve a Wathen system'
  write ( *, '(a)' ) '  defined by WATHEN_GE().'
  write ( *, '(a)' ) ''

  nx = 4
  ny = 4
  write ( *, '(a,i6)' ) '  Elements in X direction NX = ', nx
  write ( *, '(a,i6)' ) '  Elements in Y direction NY = ', ny
  write ( *, '(a,i6)' ) '  Number of elements = ', nx * ny
!
!  Compute the number of unknowns.
!
  call wathen_order ( nx, ny, n )
  write ( *, '(a,i6)' ) '  Number of nodes N = ', n
!
!  Set up a random solution X.
!
  allocate ( x1(1:n) )
  call random_number ( harvest = x1(1:n) )
!
!  Compute the matrix.
!
  allocate ( a(1:n,1:n) )
  call wathen_ge ( nx, ny, n, a )
!
!  Compute the corresponding right hand side B.
!
  allocate ( b(1:n) )
  b = matmul ( a, x1 )
!
!  Solve the linear system.
!
  allocate ( ipvt(1:n) )
  call dgefa ( a, n, n, ipvt, info )

  allocate ( x2(1:n) )
  do i = 1, n
    x2(i) = b(i)
  end do
  job = 0

  call dgesl ( a, n, n, ipvt, x2, job )
!
!  Compute the maximum solution error.
!
  e = maxval ( abs ( x1 - x2 ) )
  write ( *, '(a,g14.6)' ) '  Maximum solution error is ', e 
!
!  Free memory.
!
  deallocate ( a )
  deallocate ( b )
  deallocate ( ipvt )
  deallocate ( x1 )
  deallocate ( x2 )

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! test02() assembles, factors and solves using WATHEN_GB().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer info
  integer, allocatable :: ipvt(:)
  integer job
  integer lda
  integer md
  integer ml
  integer mu
  integer n
  integer nx
  integer ny
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test02():'
  write ( *, '(a)' ) '  Assemble, factor and solve a Wathen system'
  write ( *, '(a)' ) '  using WATHEN_GB().'
  write ( *, '(a)' ) ''

  nx = 4
  ny = 4
  write ( *, '(a,i6)' ) '  Elements in X direction NX = ', nx
  write ( *, '(a,i6)' ) '  Elements in Y direction NY = ', ny
  write ( *, '(a,i6)' ) '  Number of elements = ', nx * ny
!
!  Compute the number of unknowns.
!
  call wathen_order ( nx, ny, n )
  write ( *, '(a,i6)' ) '  Number of nodes N = ', n
!
!  Compute the bandwidth.
!
  call wathen_bandwidth ( nx, ny, ml, md, mu )
  write ( *, '(a,i6)' ) '  Lower bandwidth ML = ', ml
  write ( *, '(a,i6)' ) '  Upper bandwidth MU = ', mu
!
!  Set up a random solution X1.
!
  allocate ( x1(1:n) )
  call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
  allocate ( a(1:2*ml+mu+1,1:n) )
  call wathen_gb ( nx, ny, n, a )
!
!  Compute the corresponding right hand side B.
!
  allocate ( b(1:n) )
  call mv_gb ( n, n, ml, mu, a, x1, b )
!
!  Solve the linear system.
!
  lda = 2 * ml + mu + 1
  allocate ( ipvt(1:n) )
  call dgbfa ( a, lda, n, ml, mu, ipvt, info )

  allocate ( x2(1:n) )
  x2(1:n) = b(1:n)
  job = 0
  call dgbsl ( a, lda, n, ml, mu, ipvt, x2, job )
!
!  Compute the maximum solution error.
!
  e = maxval ( abs ( x1 - x2 ) );
  write ( *, '(a,g14.6)' ) '  Maximum solution error is ', e
!
!  Free memory.
!
  deallocate ( a )
  deallocate ( b )
  deallocate ( ipvt )
  deallocate ( x1 )
  deallocate ( x2 )

  return
end
subroutine test05 ( )

!*****************************************************************************80
!
!! test05() measures the storage needed for the Wathen system.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  integer bd1
  integer bd2
  integer bl1
  integer bl2
  integer bu1
  integer bu2
  integer bw1
  integer bw2
  integer n
  integer nx
  integer ny
  integer storage_gb
  integer storage_ge
  integer storage_sparse
  integer test

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test05():'
  write ( *, '(a)' ) '  For various problem sizes and storage schemes,'
  write ( *, '(a)' ) '  measure the storage used for the Wathen system.'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '                                   Predicted  Observed'
  write ( *, '(a)' ) '                              GE        Band      Band      Band    Sparse'
  write ( *, '(a)' ) '    NX  Elements   Nodes   storage     width     width   storage   storage'
  write ( *, '(a)' ) ''

  nx = 1
  ny = 1

  do test = 1, 6
!
!  Compute the number of unknowns.
!
    call wathen_order ( nx, ny, n )
!
!  Predict the bandwidth.
!
    call wathen_bandwidth ( nx, ny, bl1, bd1, bu1 )
    bw1 = bl1 + bd1 + bu1
!
!  Compute the matrix.
!
    allocate ( a(1:n,1:n) )
    call wathen_ge ( nx, ny, n, a )

    storage_ge = n * n

    call bandwidth ( n, n, a, bw2, bl2, bd2, bu2 )
    storage_gb = ( 2 * bl2 + 1 + bu2 ) * n

    call nonzeros ( n, n, a, storage_sparse )
!
!  Report.
!
    write ( *, '(2x,i4,6x,i4,2x,i6,2x,i8,2x,i8,2x,i8,2x,i8,2x,i8)' )  &
      nx, nx * ny, n, storage_ge, bw1, bw2, storage_gb, storage_sparse
!
!  Ready for next iteration.
!
    nx = nx * 2
    ny = ny * 2

    deallocate ( a )

  end do

  return
end
subroutine test06 ( )

!*****************************************************************************80
!
!! test06() times WATHEN_GE() assembly and solution.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer i
  integer info
  integer, allocatable :: ipvt(:)
  integer job
  integer n
  integer nx
  integer ny
  integer storage_ge
  real ( kind = rk ) t0
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  integer test
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test06():'
  write ( *, '(a)' ) '  For various problem sizes,'
  write ( *, '(a)' ) '  time the assembly and factorization of a Wathen system'
  write ( *, '(a)' ) '  using the WATHEN_GE() function.'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '    NX  Elements   Nodes   Storage    Assembly      Factor      Error'
  write ( *, '(a)' ) ''

  nx = 1
  ny = 1

  do test = 1, 6
!
!  Compute the number of unknowns.
!
    call wathen_order ( nx, ny, n )
    storage_ge = n * n
!
!  Set up a random solution X1.
!
    allocate ( x1(1:n) )
    call random_number ( harvest = x1 )
!
!  Compute the matrix, and measure the storage required.
!
    allocate ( a(1:n,1:n) )

    call cpu_time ( t0 )
    call wathen_ge ( nx, ny, n, a )
    call cpu_time ( t1 )
    t1 = t1 - t0
!
!  Compute the corresponding right hand side B.
!
    allocate ( b(1:n) )
    b = matmul ( a, x1 )
!
!  Solve the system.
!
    allocate ( ipvt(1:n) )
    allocate ( x2(1:n) )
    do i = 1, n
      x2(i) = b(i)
    end do
    job = 0

    call cpu_time ( t0 )
    call dgefa ( a, n, n, ipvt, info )
    call dgesl ( a, n, n, ipvt, x2, job )
    call cpu_time ( t2 )
    t2 = t2 - t0
!
!  Compute the maximum solution error.
!
    e = maxval ( abs ( x1 - x2 ) )
!
!  Report.
!
    write ( *, '(2x,i4,6x,i4,2x,i6,2x,i8,2x,e10.2,2x,e10.2,2x,e10.2)' ) &
      nx, nx * ny, n, storage_ge, t1, t2, e
!
!  Ready for next iteration.
!
    nx = nx * 2
    ny = ny * 2
!
!  Free memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( ipvt )
    deallocate ( x1 )
    deallocate ( x2 )

  end do

  return
end
subroutine test07 ( )

!*****************************************************************************80
!
!!  test07() times WATHEN_GB() assembly and solution.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer info
  integer, allocatable :: ipvt(:)
  integer job
  integer lda
  integer md
  integer ml
  integer mu
  integer n
  integer nx
  integer ny
  integer storage_gb
  real ( kind = rk ) t0
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  integer test
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test07()'
  write ( *, '(a)' ) '  For various problem sizes,'
  write ( *, '(a)' ) '  time the assembly and factorization of a Wathen system'
  write ( *, '(a)' ) '  using the WATHEN_GB() function.'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '    NX  Elements   Nodes   Storage    Assembly      Factor      Error'
  write ( *, '(a)' ) ''

  nx = 1
  ny = 1

  do test = 1, 6
!
!  Compute the number of unknowns.
!
    call wathen_order ( nx, ny, n )
!
!  Compute the bandwidth.
!
    call wathen_bandwidth ( nx, ny, ml, md, mu )
    storage_gb = ( 2 * ml + mu + 1 ) * n
!
!  Set up a random solution X1.
!
    allocate ( x1(1:n) )
    call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
    allocate ( a(1:2*ml+mu+1,1:n) )
    call cpu_time ( t0 )
    call wathen_gb ( nx, ny, n, a )
    call cpu_time ( t1 )
    t1 = t1 - t0
!
!  Compute the corresponding right hand side B.
!
    allocate ( b(1:n) )
    call mv_gb ( n, n, ml, mu, a, x1, b )
!
!  Solve the system.
!
    lda = 2 * ml + mu + 1
    allocate ( ipvt(1:n) )
    allocate ( x2(1:n) )
    x2(1:n) = b(1:n)
    job = 0

    call cpu_time ( t0 )
    call dgbfa ( a, lda, n, ml, mu, ipvt, info )
    call dgbsl ( a, lda, n, ml, mu, ipvt, x2, job )
    call cpu_time ( t2 )
    t2 = t2 - t0
!
!  Compute the maximum solution error.
!
    e = maxval ( abs ( x1 - x2 ) )
!
!  Report.
!
    write ( *, '(2x,i4,6x,i4,2x,i6,2x,i8,2x,e10.2,2x,e10.2,2x,e10.2)' ) &
      nx, nx * ny, n, storage_gb, t1, t2, e
!
!  Ready for next iteration.
!
    nx = nx * 2
    ny = ny * 2
!
!  Free memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( ipvt )
    deallocate ( x1 )
    deallocate ( x2 )

  end do

  return
end
subroutine test08 ( )

!*****************************************************************************80
!
!! test08() times WATHEN_GE() and WATHEN_GB().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer i
  integer info
  integer, allocatable :: ipvt(:)
  integer job
  integer lda
  integer md
  integer ml
  integer mu
  integer n
  integer nx
  integer ny
  integer storage_gb
  integer storage_ge
  real ( kind = rk ) t0
  real ( kind = rk ) t1
  real ( kind = rk ) t2
  integer test
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test08()'
  write ( *, '(a)' ) '  For various problem sizes,'
  write ( *, '(a)' ) '  time the assembly and factorization of a Wathen system'
  write ( *, '(a)' ) '  WATHEN_GE() and WATHEN_GB()'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '                   NX  Elements   Nodes   Storage    Assembly      Factor      Error'

  nx = 1
  ny = 1

  do test = 1, 6
!
!  Compute the number of unknowns.
!
    call wathen_order ( nx, ny, n )
    storage_ge = n * n
!
!  Set up a random solution X1.
!
    allocate ( x1(1:n) )
    call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
    allocate ( a(1:n,1:n) )
    call cpu_time ( t0 )
    call wathen_ge ( nx, ny, n, a )
    call cpu_time ( t1 )
    t1 = t1 - t0
!
!  Compute the corresponding right hand side B.
!
    allocate ( b(1:n) )
    b = matmul ( a, x1 )
!
!  Solve the system.
!
    allocate ( ipvt(1:n) )
    allocate ( x2(1:n) )
    do i = 1, n
      x2(i) = b(i)
    end do
    job = 0

    call cpu_time ( t0 )
    call dgefa ( a, n, n, ipvt, info )
    call dgesl ( a, n, n, ipvt, x2, job )
    call cpu_time ( t2 )
    t2 = t2 - t0
!
!  Compute the maximum solution error.
!
    e = maxval ( abs ( x1 - x2 ) )
!
!  Report.
!
    write ( *, '(a)' ) ''
    write ( *, '(2x,a,2x,i4,6x,i4,2x,i6,2x,i8,2x,e10.2,2x,e10.2,2x,e10.2)' ) &
      'WATHEN_GE    ', nx, nx * ny, n, storage_ge, t1, t2, e
!
!  Free memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( ipvt )
    deallocate ( x1 )
    deallocate ( x2 )
!
!  Compute the bandwidth.
!
    call wathen_bandwidth ( nx, ny, ml, md, mu )
    storage_gb = ( 2 * ml + mu + 1 ) * n
!
!  Set up a random solution X1.
!
    allocate ( x1(1:n) )
    call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
    allocate ( a(1:2*ml+mu+1,1:n) )
    call cpu_time ( t0 )
    call wathen_gb ( nx, ny, n, a )
    call cpu_time ( t1 )
    t1 = t1 - t0
!
!  Compute the corresponding right hand side B.
!
    allocate ( b(1:n) )
    call mv_gb ( n, n, ml, mu, a, x1, b )
!
!  Solve the system.
!
    lda = 2 * ml + mu + 1
    allocate ( ipvt(1:n) )
    allocate ( x2(1:n) )
    x2(1:n) = b(1:n)
    job = 0

    call cpu_time ( t0 )
    call dgbfa ( a, lda, n, ml, mu, ipvt, info )
    call dgbsl ( a, lda, n, ml, mu, ipvt, x2, job )
    call cpu_time ( t2 )
    t2 = t2 - t0
!
!  Compute the maximum solution error.
!
    e = maxval ( abs ( x1 - x2 ) )
!
!  Report.
!
    write ( *, '(2x,a,2x,i4,6x,i4,2x,i6,2x,i8,2x,e10.2,2x,e10.2,2x,e10.2)' ) &
      'WATHEN_GB    ', nx, nx * ny, n, storage_gb, t1, t2, e
!
!  Free memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( ipvt )
    deallocate ( x1 )
    deallocate ( x2 )
!
!  Ready for next iteration.
!
    nx = nx * 2
    ny = ny * 2

  end do

  return
end
subroutine test10 ( )

!*****************************************************************************80
!
!! test10() assembles, factors and solves using WATHEN_GE() and CG_GE().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer n
  integer nx
  integer ny
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test10()'
  write ( *, '(a)' ) '  Assemble, factor and solve a Wathen system'
  write ( *, '(a)' ) '  defined by WATHEN_GE() and CG_GE().'
  write ( *, '(a)' ) ''

  nx = 1
  ny = 1
  write ( *, '(a,i6)' ) '  Elements in X direction NX = ', nx
  write ( *, '(a,i6)' ) '  Elements in Y direction NY = ', ny
  write ( *, '(a,i6)' ) '  Number of elements = ', nx * ny
!
!  Compute the number of unknowns.
!
  call wathen_order ( nx, ny, n )
  write ( *, '(a,i6)' ) '  Number of nodes N = ', n
!
!  Set up a random solution X.
!
  allocate ( x1(1:n) )
  call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
  allocate ( a(1:n,1:n) )
  call wathen_ge ( nx, ny, n, a )
!
!  Compute the corresponding right hand side B.
!
  allocate ( b(1:n) )
  b = matmul ( a, x1 )
!
!  Solve the linear system.
!
  allocate ( x2(1:n) )
  x2(1:n) = 1.0D+00
  call cg_ge ( n, a, b, x2 )
!
!  Compute the maximum solution error.
!
  e = maxval ( abs ( x1 - x2 ) )
  write ( *, '(a,g14.6)' ) '  Maximum solution error is ', e 
!
!  Free memory.
!
  deallocate ( a )
  deallocate ( b )
  deallocate ( x1 )
  deallocate ( x2 )

  return
end
subroutine test11 ( )

!*****************************************************************************80
!
!! test11() assembles, factors and solves using WATHEN_ST() + CG_ST().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:)
  real ( kind = rk ), allocatable :: b(:)
  integer, allocatable :: col(:)
  real ( kind = rk ) e
  integer n
  integer nx
  integer ny
  integer nz_num
  integer, allocatable :: row(:)
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test11():'
  write ( *, '(a)' ) '  Assemble, factor and solve a Wathen system'
  write ( *, '(a)' ) '  defined by WATHEN_ST() and CG_ST().'
  write ( *, '(a)' ) ''

  nx = 1
  ny = 1
  write ( *, '(a,i6)' ) '  Elements in X direction NX = ', nx
  write ( *, '(a,i6)' ) '  Elements in Y direction NY = ', ny
  write ( *, '(a,i6)' ) '  Number of elements = ', nx * ny
!
!  Compute the number of unknowns.
!
  call wathen_order ( nx, ny, n )
  write ( *, '(a,i6)' ) '  Number of nodes N = ', n
!
!  Set up a random solution X1.
!
  allocate ( x1(1:n) )
  call random_number ( harvest = x1 )
!
!  Compute the matrix size.
!
  call wathen_st_size ( nx, ny, nz_num )
  write ( *, '(a,i6)' ) '  Number of nonzeros NZ_NUM = ', nz_num
!
!  Compute the matrix.
!
  allocate ( row(1:nz_num) )
  allocate ( col(1:nz_num) )
  allocate ( a(1:nz_num) )
  call wathen_st ( nx, ny, nz_num, row, col, a )
!
!  Compute the corresponding right hand side B.
!
  allocate ( b(1:n) )
  call mv_st ( n, n, nz_num, row, col, a, x1, b )
!
!  Solve the linear system.
!
  allocate ( x2(1:n) )
  x2(1:n) = 1.0D+00
  call cg_st ( n, nz_num, row, col, a, b, x2 )
!
!  Compute the maximum solution error.
!
  e = maxval ( abs ( x1 - x2 ) )
  write ( *, '(a,g14.6)' ) '  Maximum solution error is ', e 
!
!  Free memory.
!
  deallocate ( a )
  deallocate ( b )
  deallocate ( col )
  deallocate ( row )
  deallocate ( x1 )
  deallocate ( x2 )

  return
end
subroutine test115 ( )

!*****************************************************************************80
!
!! test115() assembles, factors and solves using WATHEN_GB() and CG_GB().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 June 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  real ( kind = rk ), allocatable :: a(:,:)
  real ( kind = rk ), allocatable :: b(:)
  real ( kind = rk ) e
  integer md
  integer ml
  integer mu
  integer n
  integer nx
  integer ny
  real ( kind = rk ), allocatable :: x1(:)
  real ( kind = rk ), allocatable :: x2(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'test115():'
  write ( *, '(a)' ) '  Assemble, factor and solve a Wathen system'
  write ( *, '(a)' ) '  using WATHEN_GB() and CG_GB().'
  write ( *, '(a)' ) ''

  nx = 4
  ny = 4
  write ( *, '(a,i6)' ) '  Elements in X direction NX = ', nx
  write ( *, '(a,i6)' ) '  Elements in Y direction NY = ', ny
  write ( *, '(a,i6)' ) '  Number of elements = ', nx * ny
!
!  Compute the number of unknowns.
!
  call wathen_order ( nx, ny, n )
  write ( *, '(a,i6)' ) '  Number of nodes N = ', n
!
!  Compute the bandwidth.
!
  call wathen_bandwidth ( nx, ny, ml, md, mu )
  write ( *, '(a,i6)' ) '  Lower bandwidth ML = ', ml
  write ( *, '(a,i6)' ) '  Upper bandwidth MU = ', mu
!
!  Set up a random solution X1.
!
  allocate ( x1(1:n) )
  call random_number ( harvest = x1 )
!
!  Compute the matrix.
!
  allocate ( a(1:2*ml+mu+1,1:n) )
  call wathen_gb ( nx, ny, n, a )
!
!  Compute the corresponding right hand side B.
!
  allocate ( b(1:n) )
  call mv_gb ( n, n, ml, mu, a, x1, b )
!
!  Solve the linear system.
!
  allocate ( x2(1:n) )
  x2(1:n) = 1.0D+00
  call cg_gb ( n, ml, mu, a, b, x2 )
!
!  Compute the maximum solution error.
!
  e = maxval ( abs ( x1 - x2 ) );
  write ( *, '(a,g14.6)' ) '  Maximum solution error is ', e
!
!  Free memory.
!
  deallocate ( a )
  deallocate ( b )
  deallocate ( x1 )
  deallocate ( x2 )

  return
end
subroutine wathen_xy_test ( )

!*****************************************************************************80
!
!! wathen_xy_test() tests wathen_xy().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 February 2020
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  integer i
  integer j
  integer k
  integer n
  integer nx
  integer ny
  real ( kind = rk ), allocatable :: x(:)
  real ( kind = rk ), allocatable :: y(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'wathen_xy_test():'
  write ( *, '(a)' ) '  wathen_xy() returns the (X,Y) coordinates of nodes in the'
  write ( *, '(a)' ) '  Wathen finite element system.'

  nx = 3
  ny = 3
  call wathen_order ( nx, ny, n )

  allocate ( x(1:n) )
  allocate ( y(1:n) )
  call wathen_xy ( nx, ny, n, x, y )

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '   k   i   j         x          y'
  write ( *, '(a)' ) ''

  k = 0
  do j = 1, 2 * ny + 1

    if ( mod ( j, 2 ) == 1 ) then

      do i = 1, 2 * nx + 1
        k = k + 1
        write ( *, '(2x,i2,2x,i2,2x,i2,2x,f8.4,2x,f8.4)' ) k, i, j, x(k), y(k)
      end do

    else

      do i = 1, nx + 1
        k = k + 1
        write ( *, '(2x,i2,2x,i2,2x,i2,2x,f8.4,2x,f8.4)' ) k, i, j, x(k), y(k)
      end do

    end if

  end do
 
  deallocate ( x )
  deallocate ( y )

  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
!
!  Parameters:
!
!    None
!
  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