program main

c*********************************************************************72
c
cc biharmonic_fd2d_test() tests biharmonic_fd2d().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    31 May 2021
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'biharmonic_fd2d_test()'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test biharmonic_fd2d().'

      call test1 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'biharmonic_fd2d_test()'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop 0
      end
      subroutine test1 ( )

c*********************************************************************72
c
cc test1() solves the biharmonic equation for example #1.
c
c  Discussion:
c
c    0 <= x <= 1, 0 <= y <= 1
c    uxxxx + 2 uxxyy + uyyyy = 8
c    u = 0 on boundaries
c    derivative conditions determined from exact solution.
c
c    u  = x * ( 1 - x ) * y * ( 1 - y )
c    ux = ( 1 - 2 x )   * y * ( 1 - y )
c    uy = x * ( 1 - x ) * ( 1 - 2 y )
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer m
      parameter ( m = 15 )
      integer n
      parameter ( n = 15 )

      integer lw
c     parameter ( lw = max ( 7 * n, 3 * m ) + 2 * ( n + m ) + 19 )

      parameter ( lw = max ( 3 * m, 4 * n ) 
     &  + 4 * n + 2 * m + ( n + 1 )**2 / 2 + 19 )

      double precision alpha
      double precision bda(n)
      double precision bdb(n)
      double precision bdc(m)
      double precision bdd(m)
      double precision beta
      double precision dudx
      double precision dudy
      double precision e
      double precision f(m+2,n+2)
      integer i
      integer idf
      integer iflag
      integer itcg
      integer j
      double precision tol
      double precision u
      double precision w(lw)
      double precision x
      double precision xmax
      double precision xmin
      double precision y
      double precision ymax
      double precision ymin
c
c  Rectangle limits.
c    xmin <= x <= xmax
c    ymin <= y <= ymax
c
      xmin = 0.0D+00
      xmax = 1.0D+00
      ymin = 0.0D+00
      ymax = 1.0D+00
c
c  alpha and beta are coefficients in the equation.
c
      alpha = 0.0D+00
      beta = 0.0D+00
c
c  BDA, BDB, BDC, BDD are the derivative values at the boundaries.
c
c  Here we want -dudx
c
      x = xmin
      do j = 1, n
        y = ymin + j * ( ymax - ymin ) / ( n + 1 )
        dudx = ( 1.0 - 2.0 * x ) * y * ( 1.0 - y )
        bda(j) = - dudx
      end do
c
c  Here we want dudx
c
      x = xmax
      do j = 1, n
        y = ymin + j * ( ymax - ymin ) / ( n + 1 )
        dudx = ( 1.0 - 2.0 * x ) * y * ( 1.0 - y )
        bdb(j) = dudx
      end do
c
c  Here we want -dudy
c
      y = ymin
      do i = 1, m
        x = xmin + i * ( xmax - xmin ) / ( m + 1 )
        dudy = x * ( 1.0 - x ) * ( 1.0 - 2.0 * y )
        bdc(i) = - dudy
      end do
c
c  Here we want dudy
c
      y = ymax
      do i = 1, m
        x = xmin + i * ( xmax - xmin ) / ( m + 1 )
        dudy = x * ( 1.0 - x ) * ( 1.0 - 2.0 * y )
        bdd(i) = dudy
      end do
c
c  Typical extremely confusing mapping between matrix coordinates
c  and spatial coordinates.
c
c  In "MATRIX WORLD" (I,J):
c
c    I=1   X=XMIN   F(1,1)   .....  F(1,N+2)
c          ...      .......  .....  .....
c    I=M+2 X=XMAX   F(M+2,1) .....  F(M+2,N+2)
c                   J=1 Y=YMIN      J=N+2 Y=YMAX
c
c  In "CARTESIAN WORLD (X,Y):
c
c    I=M+2 X=XMAX   F(M+2,1) .....  F(M+2,N+2)
c          ...      .......  .....  .....
c    I=1   X=XMIN   F(1,1)   .....  F(1,N+2)
c                   J=1 Y=YMIN      J=N+2 Y=YMAX
c

c
c  Interior entries of F are the right hand side.
c
      do j = 2, n + 1
        y = ( ( n + 2 - j ) * ymin + ( j - 1 ) * ymax ) / ( n + 1 )
        do i = 2, m + 1
          x = ( ( m + 2 - i ) * xmin + ( i - 1 ) * xmax ) / ( m + 1 )
          f(i,j) = 8.0D+00
        end do
      end do
c
c  Boundary entries of F are the boundary values of U.
c
      i = 1
      x = xmin
      do j = 1, n + 2
        y = ( ( n + 2 - j ) * ymin + ( j - 1 ) * ymax ) / ( n + 1 )
        f(i,j) = 0.0D+00
      end do

      i = m
      x = xmax
      do j = 1, n + 2
        y = ( ( n + 2 - j ) * ymin + ( j - 1 ) * ymax ) / ( n + 1 )
        f(i,j) = 0.0D+00
      end do

      j = 1
      y = ymin
      do i = 1, m + 2
        x = ( ( m + 2 - i ) * xmin + ( i - 1 ) * xmax ) / ( m + 1 )
        f(i,j) = 0.0D+00
      end do

      j = n
      y = ymax
      do i = 1, m + 2
        x = ( ( m + 2 - i ) * xmin + ( i - 1 ) * xmax ) / ( m + 1 )
        f(i,j) = 0.0D+00
      end do
c
c  IDF is the row dimension of F.
c
      idf = m + 2
c
c  TOL is the error tolerance for the conjugate gradient iteration.
c
      tol = 1.0D-8
c
c  IFLAG chooses the solution method.
c
      iflag = 3

      call dbihar ( xmin, xmax, m, bda, bdb, bdc, bdd, ymin, ymax, 
     &  n, f, idf, alpha, beta, iflag, tol, itcg, w, lw )

      write ( *, '(a)' ) ''
      write ( *, '(a,i8)' ) '  Return flag = ', iflag 
      write ( *, '(a,i8)' ) '  Conjugate gradient iterations = ', itcg 
c
c  Compute RMS norm of error.
c
      e = 0.0D+00
      do j = 1, n + 2
        y = ( ( n + 2 - j ) * ymin + ( j - 1 ) * ymax ) / ( n + 1 )
        do i = 1, m + 2
          x = ( ( m + 2 - i ) * xmin + ( i - 1 ) * xmax ) / ( m + 1 )
          u = x * ( 1.0 - x ) * y * ( 1.0 - y )
          e = e + ( u - f(i,j) )**2
          write ( *, * ) i, j, x, y, u, f(i,j)
        end do
      end do
      e = e / m / n
      e = sqrt ( e )
      write ( *, '(a,g14.6)' ) '  RMS error ||Exact-Computed|| = ', e

      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 June 2014
c
c  Author:
c
c    John Burkardt
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, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end