! pde46h_eq.f90  homogeneous Biharmonic equation in six dimensions
! solve Uxxxx(x,y,z,t,u,v) + Uyyyy(x,y,z,t,u,v) + Uzzzz(x,y,z,t,u,v) + 
!       Utttt(x,y,z,t,u,v) + Uuuuu(x,y,z,t,u,v) + Uwwww(x,y,z,t,u,v) +
!       2 Uxx(x,y,z,t,u,v) + 2 Uyy(x,y,z,t,u,v) + 2 Uzz(x,y,z,t,u,v) +
!       2 Utt(x,y,z,t,u,v) + 2 Uuu(x,y,z,t,u,v) + 2 Uvv(x,y,z,t,u,v) +
!       6 U(x,y,z,t,u,v) = 0
!
! boundary conditions computed using U(x,y,z,t,u,v)
! analytic solution is U(x,y,z,t,u,v) = sin(x+y+z+t+u+v)
!

module global
  implicit none
  ! problem parameters
  integer, parameter :: nx=5
  integer, parameter :: ny=5
  integer, parameter :: nz=5
  integer, parameter :: nt=5
  integer, parameter :: nu=5
  integer, parameter :: nv=5
  integer, parameter :: nxyztuv=(nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)
  double precision :: xmin = 0.0
  double precision :: xmax = 0.5
  double precision :: ymin = 0.0
  double precision :: ymax = 0.5
  double precision :: zmin = 0.0
  double precision :: zmax = 0.5
  double precision :: tmin = 0.0
  double precision :: tmax = 0.5
  double precision :: umin = 0.0
  double precision :: umax = 0.5
  double precision :: vmin = 0.0
  double precision :: vmax = 0.5
  double precision :: hx, hy, hz, ht, hu, hv
  double precision, dimension(nxyztuv,nxyztuv) :: A ! stiffness matrix
     ! A(s(i,ii,i3,i4,i5,i6),s(j,jj,j3,j4,j5,j6))
  double precision, dimension(nx) :: xg ! x grid coordinates
  double precision, dimension(ny) :: yg ! y
  double precision, dimension(nz) :: zg ! z
  double precision, dimension(nt) :: tg ! t
  double precision, dimension(nu) :: ug ! u
  double precision, dimension(nv) :: vg ! v
  double precision, dimension(nxyztuv) :: R ! RHS
  double precision, dimension(nxyztuv) :: Uk ! unknown solution
  double precision, dimension(nxyztuv) :: Ua ! known solution for checking

  interface
    function uana(x,y,z,t,u,v) result(val) ! analytic solution and boundary
      double precision, intent(in) :: x,y,z,t,u,v
      double precision :: val
    end function uana
    subroutine simeq(n, sz, A, Y, X)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz), intent(in) :: Y
      double precision, dimension(sz), intent(out) :: X
    end subroutine simeq
    subroutine rderiv1(order, npoints, point, hx, cx)
      integer, intent(in) :: order, npoints, point
      double precision, intent(in) :: hx
      double precision, dimension(npoints), intent(out) :: cx
    end subroutine rderiv1
  end interface
end module global

module formain
  interface
    subroutine buildA()
    end subroutine buildA
    subroutine check_soln(UU) ! for specific PDE
      double precision, dimension(*), intent(in) :: UU
    end subroutine check_soln
    subroutine write_soln(UU) ! for specific PDE
      double precision, dimension(*), intent(in) :: UU
    end subroutine write_soln

  end interface
end module formain

module uses
  interface
    function s(i, ii, i3, i4, i5, i6) result(val)
      integer, intent(in) :: i, ii, i3, i4, i5, i6
      integer :: val
    end function s
  end interface
end module uses

! functions called by program
  function uana(x,y,z,t,u,v) result(val) ! analytic solution and boundary
    implicit none
    double precision, intent(in) :: x,y,z,t,u,v
    double precision :: val

    val =  sin(x+y+z+t+u+v)
  end function uana

  function s(i, ii, i3, i4, i5, i6) result(val)
    use global
    implicit none
    integer, intent(in) :: i, ii, i3, i4, i5, i6
    integer :: val

    val = (i-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2) + &
          (ii-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2) + &
          (i3-2)*(nz-2)*(nt-2)*(nu-2) + &
          (i4-2)*(nt-2)*(nu-2) + &
          (i5-2)*(nu-2) + i6-1
  end function s

  subroutine buildA() ! dpends on specific PDE
    use global
    use uses
    implicit none
    integer :: i, j, ii, jj, i3, j3, i4, j4, i5, j5, i6, j6, k
    double precision :: val = 0
    double precision, dimension(nx) :: cxx, cxxxx
    double precision, dimension(ny) :: cyy, cyyyy
    double precision, dimension(nz) :: czz, czzzz
    double precision, dimension(nt) :: ctt, ctttt
    double precision, dimension(nu) :: cuu, cuuuu
    double precision, dimension(nv) :: cvv, cvvvv

    ! compute entries=A matrix, inside boundary
    print *, 'compute A matrix and Y RHS'


    print *, 'check rderiv1(4,nv,1,hv,cvvvv)'
    call rderiv1(4,nv,1,hv,cvvvv)
    do i6=1,nv 
      print *, 'cvvvv(',i6,')=',cvvvv(i6)
    end do
    print *, 'check rderiv1(4,nv,nv,hv,cvvvv)'
    call rderiv1(4,nv,nv-1,hv,cvvvv)
    do i6=1,nv
      print *, 'cvvvv(',i6,')=',cvvvv(i6)
    end do
    print *, ' '

    val = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do i3=2,nz-1
          do i4=2,nt-1
            do i5=2,nu-1
              do i6=2,nv-1
        k = s(i,ii,i3,i4,i5,i6) ! row index

        ! do each term

        ! do d^4U/dx^4
        call rderiv1(4,nx,i,hx,cxxxx)
        do j=1,nx
          val = cxxxx(j)
          if(j.eq.1 .or. j.eq.nx) then
            R(k) = R(k)-val*uana(xg(j),yg(ii),zg(i3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(j,ii,i3,i4,i5,i6))=A(k,s(j,ii,i3,i4,i5,i6)) + val
          end if
        end do ! j

        ! do d^4U/dy^4
        call rderiv1(4,ny,ii,hy,cyyyy)
        do jj=1,ny
          val = cyyyy(jj)
          if(jj.eq.1 .or. jj.eq.ny) then
            R(k) = R(k)-val*uana(xg(i),yg(jj),zg(i3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(i,jj,i3,i4,i5,i6))=A(k,s(i,jj,i3,i4,i5,i6)) + val
          end if
        end do ! jj

        ! do d^4U/dz^4
        call rderiv1(4,nz,i3,hz,czzzz)
        do j3=1,nz
          val = czzzz(j3)
          if(j3.eq.1 .or. j3.eq.nz) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(j3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(i,ii,j3,i4,i5,i6))=A(k,s(i,ii,j3,i4,i5,i6)) + val
          end if
        end do ! j3

        ! do d^4U/dt^4
        call rderiv1(4,nt,i4,ht,ctttt)
        do j4=1,nt
          val = ctttt(j4)
          if(j4.eq.1 .or. j4.eq.nt) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(j4),ug(i5),vg(i6))
          else
            A(k,s(i,ii,i3,j4,i5,i6))=A(k,s(i,ii,i3,j4,i5,i6)) + val
          end if
        end do ! j4

        ! do d^4U/du^4
        call rderiv1(4,nu,i5,hu,cuuuu)
        do j5=1,nu
          val = cuuuu(j5)
          if(j5.eq.1 .or. j5.eq.nu) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(j5),vg(i6))
          else
            A(k,s(i,ii,i3,i4,j5,i6))=A(k,s(i,ii,i3,i4,j5,i6)) + val
          end if
        end do ! j5

        ! do d^4U/dv^4
        call rderiv1(4,nv,i6,hv,cvvvv)
        do j6=1,nv
          val = cvvvv(j6)
          if(j6.eq.1 .or. j6.eq.nv) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(i5),vg(j6))
          else
            A(k,s(i,ii,i3,i4,i5,j6))=A(k,s(i,ii,i3,i4,i5,j6)) + val
          end if
        end do ! j6

        ! do 2 d^2u/dx^2
        call rderiv1(2,nx,i,hx,cxx)
        do j=1,nx
          val = 2.0*cxx(j)
          if(j.eq.1 .or. j.eq.nx) then
            R(k) = R(k)-val*uana(xg(j),yg(ii),zg(i3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(j,ii,i3,i4,i5,i6))=A(k,s(j,ii,i3,i4,i5,i6)) + val
          end if
        end do ! j

        ! do 2 d^2u/dy^2
        call rderiv1(2,ny,ii,hy,cyy)
        do jj=1,ny
          val = 2.0*cyy(jj)
          if(jj.eq.1 .or. jj.eq.ny) then
            R(k) = R(k)-val*uana(xg(i),yg(jj),zg(i3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(i,jj,i3,i4,i5,i6))=A(k,s(i,jj,i3,i4,i5,i6)) + val
          end if
        end do ! jj

        ! do 2 d^2U/dz^2
        call rderiv1(2,nz,i3,hz,czz)
        do j3=1,nz
          val = 2.0*czz(j3)
          if(j3.eq.1 .or. j3.eq.nz) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(j3),tg(i4),ug(i5),vg(i6))
          else
            A(k,s(i,ii,j3,i4,i5,i6))=A(k,s(i,ii,j3,i4,i5,i6)) + val
          end if
        end do ! j3

        ! do 2 d^2U/dt^2
        call rderiv1(2,nt,i4,ht,ctt)
        do j4=1,nt
          val = 2.0*ctt(j4)
          if(j4.eq.1 .or. j4.eq.nt) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(j4),ug(i5),vg(i6))
          else
            A(k,s(i,ii,i3,j4,i5,i6))=A(k,s(i,ii,i3,j4,i5,i6)) + val
          end if
        end do ! j4

        ! do 2 d^2U/du^2
        call rderiv1(2,nu,i5,hu,cuu)
        do j5=1,nu
          val = 2.0*cuu(j5)
          if(j5.eq.1 .or. j5.eq.nu) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(j5),vg(i6))
          else
            A(k,s(i,ii,i3,i4,j5,i6))=A(k,s(i,ii,i3,i4,j5,i6)) + val
          end if
        end do ! j5

        ! do 2 d^2U/dv^2
        call rderiv1(2,nv,i6,hv,cvv)
        do j6=1,nv
          val = 2.0*cvv(j6)
          if(j6.eq.1 .or. j6.eq.nv) then
            R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(i5),vg(j6))
          else
            A(k,s(i,ii,i3,i4,i5,j6))=A(k,s(i,ii,i3,i4,i5,j6)) + val
          end if
        end do ! j6

        ! do 6 u(x,y,z,t,u,v)
        A(k,s(i,ii,i3,i4,i5,i6))=A(k,s(i,ii,i3,i4,i5,i6)) + 6.0

        ! do f(x,y,z,t,u,v) RHS, none
        ! R(k) = R(k) + f(x,y,z,t,u,v)

        ! end terms

              end do ! i6
            end do ! i5
          end do ! i4
        end do ! i3
      end do ! ii 
    end do ! i
  end subroutine buildA

  function eval_derivative(xord, yord, zord, tord, uord, vord, &
                           i, ii, i3, i4, i5, i6, UU) result (discrete)
    use global
    use uses
    implicit none
    integer, intent(in) :: xord, yord, zord, tord, uord, vord
    integer, intent(in) :: i, ii, i3, i4, i5, i6
    double precision, dimension(*), intent(in) :: UU
    double precision :: discrete
    double precision, dimension(nx) :: cx
    double precision, dimension(ny) :: cy
    double precision, dimension(nz) :: cz
    double precision, dimension(nt) :: ct
    double precision, dimension(nu) :: cu
    double precision, dimension(nv) :: cv

    double precision :: x, y, z, t, u, v
    integer :: j, jj, j3, j4, j5, j6

    call rderiv1(xord, nx, i, hx, cx)
    x = xg(i)
    call rderiv1(yord, ny, ii, hy, cy)
    y = yg(ii)
    call rderiv1(zord, nz, i3, hz, cz)
    z = zg(i3)
    call rderiv1(tord, nt, i4, ht, ct)
    t = tg(i4)
    call rderiv1(uord, nu, i5, hu, cu)
    u = ug(i5)
    call rderiv1(vord, nv, i6, hv, cv)
    v = vg(i6)
    discrete = 0.0

    ! six cases of one partial  x, y, z, t, u, v  to any order
    if(xord.ne.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.eq.0 &
       .and. uord.eq.0 .and. vord.eq.0) then
      ! Ux  Uxx  Uxxx  Uxxxx
      do j=1,nx
        if(j.eq.1 .or. j.eq.nx) then
          discrete = discrete + cx(j)*uana(xg(j),y,z,t,u,v)
        else
          discrete = discrete + cx(j)*UU(s(j,ii,i3,i4,i5,i6))
        end if
      end do ! j

    else if(xord.eq.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.eq.0 &
            .and. uord.eq.0 .and. vord.eq.0) then
      ! Uy  Uyy  Uyyy  Uyyyy
      do jj=1,ny
        if(jj.eq.1 .or. jj.eq.ny) then
          discrete = discrete + cy(jj)*uana(x,yg(jj),z,t,u,v)
        else
          discrete = discrete + cy(jj)*UU(s(i,jj,i3,i4,i5,i6))
        end if
      end do ! jj

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.eq.0 &
            .and. uord.eq.0 .and. vord.eq.0) then
      ! Uz  Uzz  Uzzz  Uzzzz
      do j3=1,nz
        if(j3.eq.1 .or. j3.eq.nz) then
          discrete = discrete + cz(j3)*uana(x,y,zg(j3),t,u,v)
        else
          discrete = discrete + cz(j3)*UU(s(i,ii,j3,i4,i5,i6))
        end if
      end do ! j3

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.ne.0 &
            .and. uord.eq.0 .and. vord.eq.0) then
      ! Ut  Utt  Uttt  Utttt
      do j4=1,nt
        if(j4.eq.1 .or. j4.eq.nt) then
          discrete = discrete + ct(j4)*uana(x,y,z,ug(j4),u,v)
        else
           discrete = discrete + ct(j4)*UU(s(i,ii,i3,j4,i5,i6))
        end if
      end do ! j4

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.eq.0 &
            .and. uord.ne.0 .and. vord.eq.0) then
      ! Uu  Uuu  Uuuu  Uuuuu
      do j5=1,nu
        if(j5.eq.1 .or. j5.eq.nu) then
          discrete = discrete + cu(j5)*uana(x,y,z,t,ug(j5),v)
        else
           discrete = discrete + cu(j5)*UU(s(i,ii,i3,i4,j5,i6))
        end if
      end do ! j5

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.eq.0 &
            .and. uord.eq.0 .and. vord.ne.0) then
      ! Uv  Uvv  Uvvv  Uvvvv
      do j6=1,nv
        if(j6.eq.1 .or. j6.eq.nv) then
          discrete = discrete + cv(j6)*uana(x,y,z,t,u,vg(j6))
        else
           discrete = discrete + cv(j6)*UU(s(i,ii,i3,i4,i5,j6))
        end if
      end do ! j6

    end if 
    ! return discrete
  end function eval_derivative

  subroutine check_soln(UU) ! for specific PDE
    ! solve Uxxxx(x,y,z,t,u,v) + Uyyyy(x,y,z,t,u,v) + Uzzzz(x,y,z,t,u,v) +
    !       Utttt(x,y,z,t,u,v) + Uuuuu(x,y,z,t,u,v) + Uvvvv(x,y,z,t,u,v) +
    !       2 Uxx(x,y,z,t,u,v) + 2 Uyy(x,y,z,t,u,v) + 2 Uzz(x,y,z,t,u,v) +
    !       2 Utt(x,y,z,t,u,v) + 2 Uuu(x,y,z,t,u,v) + 2 Uvv(x,y,z,t,u,v) + 
    !       6 U(x,y,z,t,u,v) = 0

    use global
    implicit none
    interface
      function eval_derivative(xord, yord, zord, tord, uord, vord, &
                               i, ii, i3, i4, i5, i6, UU) result (discrete)
        integer, intent(in) :: xord, yord, zord, tord, uord, vord
        integer, intent(in) :: i, ii, i3, i4, i5, i6
        double precision, dimension(*), intent(in) :: UU
        double precision :: discrete
      end function eval_derivative
    end interface
    double precision, dimension(*), intent(in) :: UU
    double precision :: val, err, smaxerr
    integer :: i, ii, i3, i4, i5, i6

    smaxerr = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do i3=2,nz-1
          do i4=2,nt-1
            do i5=2,nu-1
              do i6=2,nv-1
              val = 0.0
              ! Uxxxx(x,y,z,t,u,v)
              val = val + eval_derivative(4,0,0,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + Uyyyy(x,y,z,t,u,v)
              val = val + eval_derivative(0,4,0,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + Uzzzz(x,y,z,t,u,v)
              val = val + eval_derivative(0,0,4,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + Utttt(x,y,z,t,u,v)
              val = val + eval_derivative(0,0,0,4,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + Uuuuu(x,y,z,t,u,v)
              val = val + eval_derivative(0,0,0,0,4,0,i,ii,i3,i4,i5,i6,UU)
              ! + Uvvvv(x,y,z,t,u,v)
              val = val + eval_derivative(0,0,0,0,0,4,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Uxx(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(2,0,0,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Uyy(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(0,2,0,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Uzz(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(0,0,2,0,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Utt(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(0,0,0,2,0,0,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Uuu(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(0,0,0,0,2,0,i,ii,i3,i4,i5,i6,UU)
              ! + 2 Uvv(x,y,z,t,u,v)
              val = val + 2.0*eval_derivative(0,0,0,0,0,2,i,ii,i3,i4,i5,i6,UU)
              ! + 5 U(x,y,z,t,u,v)
              val = val + 6.0*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(i5),vg(i6))
              ! -       no f(x,y,z,t,u,v) should be zero
              err = abs(val)
              if(err>smaxerr) then
                smaxerr = err
              end if
              end do !i6
            end do !i5
          end do !i4
        end do !i3
      end do ! ii
    end do ! i
    print *, 'check_soln maxerr=',smaxerr
  end subroutine check_soln

  subroutine write_soln(UU) ! for specific PDE
    use global
    use uses
    implicit none
    double precision, dimension(*), intent(in) :: UU
    integer :: i, ii, i3, i4, i5, i6

    print *, 'write file pde46h_eq_f90.dat'
    open (unit = 11, file = "pde46h_eq_f90.dat", action = "write")
    write(unit = 11, fmt = "(7I3)")&
          6,nx,ny,nz,nt,nu,nv
    do i=2,nx-1
      do ii=2,ny-1
        do i3=2,nz-1
          do i4=2,nt-1
            do i5=2,nu-1
              do i6=2,nv-1
            write(unit = 11, fmt = "(7F12.6)")&
                  xg(i),yg(ii),zg(i3),tg(i4),ug(i5),vg(i6),&
                  UU(s(i,ii,i3,i4,i5,i6))
              end do !i6
            end do !i5
          end do !i4
        end do !i3
      end do ! ii
    end do ! i
    print *, 'file written '
  end subroutine write_soln

program pde46h_eq
  use global
  use formain
  use uses
  implicit none
  integer :: i, j, ii, jj, i3, j3, i4, j4, i5, j5, i6, j6
  integer :: k
  double precision :: x, y, z, t, u, v
  double precision :: val, err, avgerr, maxerr
  double precision :: start, now, prev ! for time
  integer :: clock_count, clock_rate, clock_max

  print *, 'pde46h_eq.f90 running'
  call system_clock(count=clock_count, count_rate=clock_rate, &
                    count_max=clock_max)
  start = (1.0D0+clock_count)/clock_rate
  prev = start
  print *,' solve Uxxxx(x,y,z,t,u,v) + Uyyyy(x,y,z,t,u,v) + Uzzzz(x,y,z,t,u,v) +'
  print *,'       Utttt(x,y,z,t,u,v) + Uuuuu(x,y,z,t,u,v) + Uvvvv(x,y,z,t,u,v) +'
  print *,'       2 Uxx(x,y,z,t,u,v) + 2 Uyy(x,y,z,t,u,v) + 2 Uzz(x,y,z,t,u,v) +'
  print *,'       2 Utt(x,y,z,t,u,v) + 2 Uuu(x,y,z,t,u,v) + 2 Uvv(x,y,z,t,u,v) +' 
  print *,'       6 U(x,y,z,t,u,v) = 0 '
  print *,' '
  print *,' boundary conditions computed using U(x,y,z,t,u,v)'
  print *,' analytic solution is U(x,y,z,t,u,v) = sin(x+y+z+t+u+v)'
  print *, ' '
  print *, 'xmin=',xmin,', xmax =',xmax
  print *, 'ymin=',ymin,', ymax =',ymax
  print *, 'zmin=',zmin,', zmax =',zmax
  print *, 'tmin=',tmin,', tmax =',tmax
  print *, 'umin=',umin,', umax =',umax
  print *, 'vmin=',vmin,', vmax =',vmax
  print *, 'nx=',nx,', ny=',ny,', nz=',nz
  print *, 'nt=',nt,', nu=',nu,', nv=',nv 

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymin,zmin,tmin,umin,vmin'
  do i=1,nx 
    xg(i) = xmin+(i-1)*hx
    Ua(i) = uana(xg(i),ymin,zmin,tmin,umin,vmin) ! not saved
    print *, 'i=',i,', Ua(',xg(i),')=',Ua(i)
  end do
  print *, ' '

  hy = (ymax-ymin)/(ny-1) ! does not have to be uniform spacing
  print *, 'y grid and analytic solution, at xmin,zmin,tmin,umin,vmin'
  do ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    Ua(ii) = uana(xmin,yg(ii),zmin,tmin,umin,vmin)
    print *, 'ii=',ii,', Ua(',yg(ii),')=',Ua(ii)
  end do
  print *, ' '

  hz = (zmax-zmin)/(nz-1) ! does not have to be uniform spacing
  print *, 'z grid and analytic solution, at xmin,ymin,tmin,umin,vmin'
  do i3=1,nz 
    zg(i3) = zmin+(i3-1)*hz
    Ua(i3) = uana(xmin,ymin,zg(i3),tmin,umin,vmin)
    print *, 'i3=',i3,', Ua(',zg(i3),')=',Ua(i3)
  end do
  print *, ' '

  ht = (tmax-tmin)/(nt-1) ! does not have to be uniform spacing
  print *, 't grid and analytic solution, at xmin,ymin,zmin,umin,vmin'
  do i4=1,nt 
    tg(i4) = tmin+(i4-1)*ht
    Ua(i4) = uana(xmin,ymin,zmin,tg(i4),umin,vmin)
    print *, 'i4=',i4,', Ua(',tg(i4),')=',Ua(i4)
  end do
  print *, ' '

  hu = (umax-umin)/(nu-1) ! does not have to be uniform spacing
  print *, 'u grid and analytic solution, at xmin,ymin,zmin,tmin,vmin'
  do i5=1,nu
    ug(i5) = umin+(i5-1)*hu
    Ua(i5) = uana(xmin,ymin,zmin,tmin,ug(i5),vmin) ! not saved
    print *, 'i5=',i5,', Ua(',ug(i5),')=',Ua(i5)
  end do
  print *, ' '

  hv = (vmax-vmin)/(nv-1) ! does not have to be uniform spacing
  print *, 'v grid and analytic solution, at xmin,ymin,zmin,tmin,umin'
  do i6=1,nv
    vg(i6) = vmin+(i6-1)*hv
    Ua(i6) = uana(xmin,ymin,zmin,tmin,umin,vg(i6)) ! not saved
    print *, 'i6=',i6,', Ua(',ug(i6),')=',Ua(i6)
  end do
  print *, ' '

  ! put solution in: Ua vector
  do i=2,nx-1
    x = xg(i)
    do ii=2,ny-1
      y = yg(ii)
      do i3=2,nz-1
        z = zg(i3)
        do i4=2,nt-1
          t = tg(i4)
          do i5=2,nu-1
            u = ug(i5)
            do i6=2,nv-1
              v = vg(i6)
              Ua(s(i,ii,i3,i4,i5,i6)) = uana(x,y,z,t,u,v)
            end do ! i6
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i

  ! initialize k and f  inside boundary
  do i=2,nx-1
    do j=2,nx-1 
      do ii=2,ny-1
        do jj=2,ny-1 
          do i3=2,nz-1 
            do j3=2,nz-1 
              do i4=2,nt-1 
                do j4=2,nt-1 
                  do i5=2,nu-1 
                    do j5=2,nu-1 
                      do i6=2,nv-1 
                        do j6=2,nv-1 
                          A(s(i,ii,i3,i4,i5,i6),s(j,jj,j3,j4,j5,j6)) = 0.0
                        end do ! j6
                      end do ! i6
                    end do ! j5
                  end do ! i5
                end do ! j4
              end do ! i4
            end do ! j3
          end do ! i3
        end do ! jj
      end do ! ii
    end do ! j
  end do ! i

  do i=2,nx-1
    do ii=2,ny-1
      do i3=2,nz-1 
        do i4=2,nt-1 
          do i5=2,nu-1 
            do i6=2,nv-1 
              R(s(i,ii,i3,i4,i5,i6)) = 0.0
            end do ! i6
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i

  ! compute entries in stiffness matrix
  call buildA()
  call system_clock(count=clock_count)
  now = (1.0D0+clock_count)/clock_rate
  print *, 'time to build equations ',now-prev,' seconds'
  prev = now

  print *, 'solve ',nxyztuv,' equations in ',nxyztuv,' unknowns'
  call simeq(nxyztuv, nxyztuv, A, R, Uk)
  call system_clock(count=clock_count)
  now = (1.0D0+clock_count)/clock_rate
  print *, 'time to solve simultaneous equations ',now-prev,' seconds'
  prev = now

  print *, 'check check_soln against solution'
  call check_soln(Ua)
  print *, 'check computed solution'
  call check_soln(Uk)
  call write_soln(Uk)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'U computed discrete, Ua analytic, error'
  do i=2,nx-1 
    do ii=2,ny-1 
      do i3=2,nz-1 
        do i4=2,nt-1 
          do i5=2,nu-1 
            do i6=2,nv-1 
      err = abs(Uk(s(i,ii,i3,i4,i5,i6))-Ua(s(i,ii,i3,i4,i5,i6)))
      if(err>maxerr) then
        maxerr = err
      end if
      avgerr = avgerr + err
      print "('U(',I1,',',I1,',',I1,',',I1,',',I1,',',I1,')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
              i,ii,i3,i4,i5,i6,Uk(s(i,ii,i3,i4,i5,i6)),&
              Ua(s(i,ii,i3,i4,i5,i6)), &
              (Uk(s(i,ii,i3,i4,i5,i6))-Ua(s(i,ii,i3,i4,i5,i6)))
            end do ! i6
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i

  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nxyztuv)
  call system_clock(count=clock_count)
  now = (1.0D0+clock_count)/clock_rate
  print *, 'end pde46h_eq.f90  in ',now-start,' seconds'

end program pde46h_eq