! pde44_sparse.f90  build discrete system of sparse equations and solve
! solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
!     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
!     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
!     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
!     F(x,y,z,t)
!     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
!                  180 z^2*t + 200 y^2z*t + 44 t^3 +
!                  110 y^2z^2t + 66 xy + 12t^4 +
!                  24 z^4 + 36 y^4 + 48 x^4 +
!                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
!
! boundary conditions computed using u(x,y,z,t)
! analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
!                                   5 y^2 z^2 t^2 + 6 x y t + 
!                                   7 x z + 8 t + 9
! solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
! A * U = R, solve simultaneous equations for U

module global
  use sparse
  implicit none
  ! problem parameters
  integer, parameter :: nx=6
  integer, parameter :: ny=6
  integer, parameter :: nz=6
  integer, parameter :: nt=6
  integer, parameter :: nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2)
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: ymin = 0.0
  double precision :: ymax = 1.0
  double precision :: zmin = 0.0
  double precision :: zmax = 1.0
  double precision :: tmin = 0.0
  double precision :: tmax = 1.0
  double precision :: hx, hy, hz, ht
  type(cell), dimension(:), pointer :: A ! sparse stiffness matrix
  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(nxyzt) :: R ! RHS
  double precision, dimension(nxyzt) :: U ! unknown solution
  double precision, dimension(nxyzt) :: Ua ! known solution for checking

  interface
    function FF(x,y,z,t) result(v) ! forcing function
      double precision, intent(in) :: x,y,z,t
      double precision :: v
    end function FF
    function uana(x,y,z,t) result(v) ! analytic solution and boundary
      double precision, intent(in) :: x,y,z,t
      double precision :: v
    end function uana
    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
  end interface
end module formain

module uses
  interface
    function s(i, ii, iii, iiii) result(v)
      integer, intent(in) :: i, ii, iii, iiii
      integer :: v
    end function s
  end interface
end module uses

module evald
  interface
    function eval_derivative(xord, yord, zord, tord, &
                             i, ii, iii, iiii, UU) result (discrete)
      integer, intent(in) :: xord, yord, zord, tord, i, ii, iii, iiii
      double precision, dimension(*), intent(in) :: UU
      double precision :: discrete
    end function eval_derivative
  end interface
end module evald

! functions called by program
  function FF(x,y,z,t) result(v) ! forcing function
    implicit none
    double precision, intent(in) :: x,y,z,t
    double precision :: v

    v= 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t + &
       200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y + &
       12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x + &
       60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t
  end function FF

  function uana(x,y,z,t) result(v) ! analytic solution and boundary
    implicit none
    double precision, intent(in) :: x,y,z,t
    double precision :: v

    v =  t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x + &
         5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0
  end function uana

  function s(i, ii, iii, iiii) result(v)
    use global
    implicit none
    integer, intent(in) :: i, ii, iii, iiii
    integer :: v

    v = (i-2)*(ny-2)*(nz-2)*(nt-2)+(ii-2)*(nz-2)*(nt-2)+(iii-2)*(nt-2)+iiii-1
  end function s

  subroutine buildA() ! dpends on specific PDE
    use global
    use uses
    implicit none
    integer :: i, j, ii, jj, iii, jjj, iiii, jjjj, k
    double precision :: val = 0
    double precision, dimension(nx) :: cx,      cxxx, cxxxx
    double precision, dimension(ny) :: cy, cyy,       cyyyy
    double precision, dimension(nz) :: cz, czz,       czzzz
    double precision, dimension(nt) :: ct, ctt,       ctttt
    double precision :: coefdxxxx, coefdyyyy, coefdzzzz, coefdtttt
    double precision :: coefdx, coefdy, coefdz, coefdt
    double precision :: coefdxxx, coefdyy, coefdzz, coefdtt
    double precision :: coefdxyt, coefdyyzz, coefdzzt, coefdztt
    double precision :: coefdyyt, coefdzt, coefu

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


    print *, 'check rderiv1(4,nx,1,hx,cxxxx)'
    call rderiv1(4,nx,1,hx,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    print *, 'check rderiv1(4,nx,nx,hx,cxxxx)'
    call rderiv1(4,nx,nx,hx,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    print *, ' '

    val = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do iii=2,nz-1
          do iiii=2,nt-1
            k = s(i,ii,iii,iiii) ! row index

            ! do each term=first equation

            ! do d^4u/dx^4
            coefdxxxx = 1.0
            call rderiv1(4,nx,i,hx,cxxxx)
            do j=1,nx
              val = coefdxxxx*cxxxx(j)
              if(j.eq.1 .or. j.eq.nx) then
                R(k) = R(k) - val*uana(xg(j),yg(ii),zg(iii),tg(iiii))
              else
                call add(A,k,s(j,ii,iii,iiii),val)
              end if
            end do ! j

            ! do d^4u/dy^4
            coefdyyyy = 2.0
            call rderiv1(4,ny,ii,hy,cyyyy)
            do jj=1,ny
              val = coefdyyyy*cyyyy(jj)
              if(jj.eq.1 .or. jj.eq.ny) then
                R(k) = R(k) - val*uana(xg(i),yg(jj),zg(iii),tg(iiii))
              else
                call add(A,k,s(i,jj,iii,iiii),val)
              end if
            end do ! jj

            ! do d^4u/dz^4
            coefdzzzz = 3.0
            call rderiv1(4,nz,iii,hz,czzzz)
            do jjj=1,nz
              val = coefdzzzz*czzzz(jjj)
              if(jjj.eq.1 .or. jjj.eq.nz) then
                R(k) = R(k) - val*uana(xg(i),yg(ii),zg(jjj),tg(iiii))
              else
                call add(A,k,s(i,ii,jjj,iiii),val)
              end if
            end do ! jjj

            ! do d^4u/dt^4
            coefdtttt = 4.0
            call rderiv1(4,nt,iiii,ht,ctttt)
            do jjjj=1,nt
              val = coefdtttt*ctttt(jjjj)
              if(jjjj.eq.1 .or. jjjj.eq.nt) then
                R(k) = R(k) - val*uana(xg(i),yg(ii),zg(iii),tg(jjjj))
              else
                call add(A,k,s(i,ii,iii,jjjj),val)
              end if
            end do ! jjjj

            ! do d^3u/dxdydt
            coefdxyt  = 5.0
            call rderiv1(1,nx,i,hx,cx)
            call rderiv1(1,ny,ii,hy,cy)
            call rderiv1(1,nt,iiii,ht,ct)
            do j=1,nx
              do jj=1,ny
                do jjjj=1,nt
                  val = coefdxyt*cx(j)*cy(jj)*ct(jjjj)
                  if(j.eq.1 .or. j.eq.nx .or. jj.eq.1 .or. jj.eq.ny &
                     .or. jjjj.eq.1 .or. jjjj.eq.nt) then
                    R(k) = R(k) - val*uana(xg(j),yg(jj),zg(iii),tg(jjjj))
                  else
                    call add(A,k,s(j,jj,iii,jjjj),val)
                  end if
                end do ! jjjj
              end do ! jj
            end do ! j

            ! do d^4u/dy^2dz^2
            coefdyyzz = 6.0
            call rderiv1(2,ny,ii,hy,cyy)
            call rderiv1(2,nz,iii,hz,czz)
            do jj=1,ny
              do jjj=1,nz
                val = coefdyyzz*cyy(jj)*czz(jjj)
                if(jj.eq.1 .or. jj.eq.ny .or. jjj.eq.1 .or. jjj.eq.nz) then
                  R(k) = R(k) - val*uana(xg(i),yg(jj),zg(jjj),tg(iiii))
                else
                  call add(A,k,s(i,jj,jjj,iiii),val)
                end if
              end do ! jjj
            end do ! jj

            ! do d^3u/dzdt^2
            coefdztt  = 7.0
            call rderiv1(1,nz,iii,hz,cz)
            call rderiv1(2,nt,iiii,ht,ctt)
            do jjj=1,nz
              do jjjj=1,nt
                val = coefdztt*cz(jjj)*ctt(jjjj)
                if(jjj.eq.1 .or. jjj.eq.nz .or. jjjj.eq.1 .or. jjjj.eq.nt) then
                  R(k) = R(k) - val*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
                else
                  call add(A,k,s(i,ii,jjj,jjjj),val)
                end if
              end do ! jjjj
            end do ! jjj

            ! do d^3u/dx^3
            coefdxxx  = 8.0
            call rderiv1(3,nx,i,hx,cxxx)
            do j=1,nx
              val = coefdxxx*cxxx(j)
              if(j.eq.1 .or. j.eq.nx) then
                R(k) = R(k) - val*uana(xg(j),yg(ii),zg(iii),tg(iiii))
              else
                call add(A,k,s(j,ii,iii,iiii),val)
              end if
            end do ! j

            ! do d^3u/dy^2dt
            coefdyyt  = 9.0
            call rderiv1(2,ny,ii,hy,cyy)
            call rderiv1(1,nt,iiii,ht,ct)
            do jj=1,ny
              do jjjj=1,nt
                val = coefdyyt*cyy(jj)*ct(jjjj)
                if(jj.eq.1 .or. jj.eq.ny .or. jjjj.eq.1 .or. jjjj.eq.nt) then
                  R(k) = R(k) - val*uana(xg(i),yg(jj),zg(iii),tg(jjjj))
                else
                  call add(A,k,s(i,jj,iii,jjjj),val)
                end if
              end do ! jjjj
            end do ! jj

            ! do d^2u/dzdt
            coefdzt   = 10.0
            call rderiv1(1,nz,iii,hz,cz)
            call rderiv1(1,nt,iiii,ht,ct)
            do jjj=1,nz
              do jjjj=1,nt
                val = coefdzt*cz(jjj)*ct(jjjj)
                if(jjj.eq.1 .or. jjj.eq.nz .or. jjjj.eq.1 .or. jjjj.eq.nt) then
                  R(k) = R(k) - val*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
                else
                  call add(A,k,s(i,ii,jjj,jjjj),val)
                end if
              end do ! jjjj
            end do ! jjj

            ! do du/dt
            coefdt    = 11.0
            call rderiv1(1,nt,iiii,ht,ct)
            do jjjj=1,nt
              val = coefdt*ct(jjjj)
              if(jjjj.eq.1 .or. jjjj.eq.nt) then
                R(k) = R(k) - val*uana(xg(i),yg(ii),zg(iii),tg(jjjj))
              else
                call add(A,k,s(i,ii,iii,jjjj),val)
              end if
            end do ! jjjj
            ! do u(x,y,z,t)
            coefu     = 12.0
            call add(A,k,s(i,ii,iii,iiii),coefu)

            ! do f(x,y,z,t) RHS constant
            R(k) = R(k) + FF(xg(i),yg(ii),zg(iii),tg(iiii))

            ! end terms

          end do ! iiii
        end do ! iii
      end do ! ii 
    end do ! i
  end subroutine buildA

  function eval_derivative(xord, yord, zord, tord, &
                           i, ii, iii, iiii, UU) result (discrete)
    use global
    use uses
    implicit none
    integer, intent(in) :: xord, yord, zord, tord, i, ii, iii, iiii
    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

    integer, parameter ::nn =  max(nx, max(ny, max(nz, nt)))
    double precision, dimension(nn) :: p1, p2, p3
    double precision :: x, y, z, t
    integer :: j, jj, jjj, jjjj

    call rderiv1(xord, nx, i, hx, cx)
    x = xg(i)
    call rderiv1(yord, ny, ii, hy, cy)
    y = yg(ii)
    call rderiv1(zord, nz, iii, hz, cz)
    z = zg(iii)
    call rderiv1(tord, nt, iiii, ht, ct)
    t = tg(iiii)
    discrete = 0.0

    ! four cases of one partial x, y, z, t to any order
    if(xord.ne.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.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)
        else
          discrete = discrete + cx(j)*UU(s(j,ii,iii,iiii))
        end if
      end do ! j
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.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)
        else
          discrete = discrete + cy(jj)*UU(s(i,jj,iii,iiii))
        end if
      end do ! jj
    else if(xord.eq.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Uz  Uzz  Uzzz  Uzzzz
      do jjj=1,nz
        if(jjj.eq.1 .or. jjj.eq.nz) then
          discrete = discrete + cz(jjj)*uana(x,y,zg(jjj),t)
        else
          discrete = discrete + cz(jjj)*UU(s(i,ii,jjj,iiii))
        end if
      end do ! jjj
    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Ut  Utt  Uttt  Utttt
      do jjjj=1,nt
        if(jjjj.eq.1 .or. jjjj.eq.nt) then
          discrete = discrete + ct(jjjj)*uana(x,y,z,tg(jjjj))
        else
           discrete = discrete + ct(jjjj)*UU(s(i,ii,iii,jjjj))
        end if
      end do ! jjjj

    ! six cases of two partials xy, xz, xt, yz, yt, zt to any order
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.eq.0) then
      ! Ux^n y^m
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          if(j.eq.1 .or. j.eq.nx .or. jj.eq.1 .or. jj.eq.ny) then
            p1(j) = p1(j) + cy(jj)*uana(xg(j),yg(jj),zg(iii),tg(iiii))
          else
             p1(j) = p1(j) + cy(jj)*UU(s(j,jj,iii,iiii))
          end if
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Ux^n z^m
      do j=1,nx
        p1(j) = 0.0
        do jjj=1,nz
          if(j.eq.1 .or. j.eq.nx .or. jjj.eq.1 .or. jjj.eq.nz) then
            p1(j) = p1(j) + cz(jjj)*uana(xg(j),yg(ii),zg(jjj),tg(iiii))
          else
            p1(j) = p1(j) + cz(jjj)*UU(s(j,ii,jjj,iiii))
          end if
        end do ! jjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Ux^n t^m
      do j=1,nx
        p1(j) = 0.0
        do jjjj=1,nt
          if(j.eq.1 .or. j.eq.nx .or. jjjj.eq.1 .or. jjjj.eq.nt) then
            p1(j) = p1(j) + ct(jjjj)*uana(xg(j),yg(ii),zg(iii),tg(jjjj))
          else
            p1(j) = p1(j) + ct(jjjj)*UU(s(j,ii,iii,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Uy^n z^m
      do jj=1,ny
        p1(jj) = 0.0
        do jjj=1,nz
          if(jj.eq.1 .or. jj.eq.ny .or. jjj.eq.1 .or. jjj.eq.nz) then
            p1(jj) = p1(jj) + cz(jjj)*uana(xg(i),yg(jj),zg(jjj),tg(iiii))
          else
            p1(jj) = p1(jj) + cz(jjj)*UU(s(i,jj,jjj,iiii))
          end if
        end do ! jjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Uy^n t^m
      do jj=1,ny
        p1(jj) = 0.0
        do jjjj=1,nt
          if(jj.eq.1 .or. jj.eq.ny .or. jjjj.eq.1 .or. jjjj.eq.nt) then
             p1(jj) = p1(jj) + ct(jjjj)*uana(xg(i),yg(jj),zg(iii),tg(jjjj))
          else
             p1(jj) = p1(jj) + ct(jjjj)*UU(s(i,jj,iii,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj
    else if(xord.eq.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Uz^n t^m
      do jjj=1,nz
        p1(jjj) = 0.0
        do jjjj=1,nt
          if(jjj.eq.1 .or. jjj.eq.nz .or. jjjj.eq.1 .or. jjjj.eq.nt) then
            p1(jjj) = p1(jjj) + ct(jjjj)*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
          else
            p1(jjj) = p1(jjj) + ct(jjjj)*UU(s(i,ii,jjj,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cz(jjj)*p1(jjj)
      end do ! jjj

    ! four cases of three partials xyz, xyt, xzt, yzt to any order
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Ux^n y^m z^p
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjj=1,nz
            if(j.eq.1 .or. j.eq.nx .or.      &
               jj.eq.1 .or. jj.eq.ny .or.    &
               jjj.eq.1 .or. jjj.eq.nz) then
               p2(jj) = p2(jj) + cz(jjj)*uana(xg(j),yg(jj),zg(jjj),tg(iiii))
            else
               p2(jj) = p2(jj) + cz(jjj)*UU(s(j,jj,jjj,iiii))
            end if
          end do ! jjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Ux^n y^m t^p
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjjj=1,nt
            if(j.eq.1 .or. j.eq.nx .or.   &
               jj.eq.1 .or. jj.eq.ny .or. &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jj) = p2(jj) + ct(jjjj)*uana(xg(j),yg(jj),zg(iii),tg(jjjj))
            else
              p2(jj) = p2(jj) + ct(jjjj)*UU(s(j,jj,iii,jjjj))
            end if
          end do ! jjjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete +cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Ux^n z^m t^p
      do j=1,nx
        p1(j) = 0.0
        do jjj=1,nz
          p2(jjj) = 0.0
          do jjjj=1,nt
            if(j.eq.1 .or. j.eq.nx .or.      &
               jjj.eq.1 .or. jjj.eq.nz .or.  &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jjj) = p2(jjj) + ct(jjjj)*uana(xg(j),yg(ii),zg(jjj),tg(jjjj))
            else
              p2(jjj) = p2(jjj) + ct(jjjj)*UU(s(j,ii,jjj,jjjj))
            end if
          end do ! jjjj
          p1(j) = p1(j) + cz(jjj)*p2(jjj)
        end do ! jjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Uy^n z^m t^p
      do jj=1,ny
        p1(jj) = 0.0
        do jjj=1,nz
          p2(jjj) = 0.0
          do jjjj=1,nt
            if(jj.eq.1 .or. jj.eq.ny .or.   &
               jjj.eq.1 .or. jjj.eq.nz .or. &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jjj) = p2(jjj) + ct(jjjj)*uana(xg(i),yg(jj),zg(jjj),tg(jjjj))
            else
              p2(jjj) = p2(jjj) + ct(jjjj)*UU(s(i,jj,jjj,jjjj))
            end if
          end do ! jjjj
          p1(jj) = p1(jj) + cz(jjj)*p2(jjj)
        end do ! jjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj

    ! one case of four partials
    else
      ! Uxyzt
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjj=1,nz
            p3(jjj) = 0.0
            do jjjj=1,nt
              if(j.eq.1 .or. j.eq.nx .or.     &
                 jj.eq.1 .or. jj.eq.ny .or.   &
                 jjj.eq.1 .or. jjj.eq.nz .or. &
                 jjjj.eq.1 .or. jjjj.eq.nt) then
                p3(jjj) = p3(jjj) + ct(jjjj)*uana(xg(j),yg(jj),zg(jjj),tg(jjjj))
              else
                p3(jjj) =  p3(jjj) + ct(jjjj)*UU(s(j,jj,jjj,jjjj))
              end if
            end do ! jjjj
            p2(jj) = p2(jj) + cz(jjj)*p3(jjj)
          end do ! jjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    end if
    ! return discrete
  end function eval_derivative

  subroutine check_soln(UU) ! for specific PDE
    ! solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
    !     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
    !     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
    !     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
    !     F(x,y,z,t)
    use global
    use uses
    use evald
    implicit none
    double precision, dimension(*), intent(in) :: UU
    double precision :: val, err, smaxerr
    integer :: i, ii, iii, iiii

    smaxerr = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do iii=2,nz-1
          do iiii=2,nt-1
            val = 0.0
            ! uxxxx(x,y,z,t)
            val = val + eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, UU)
            ! + 2 uyyyy(x,y,z,t)
            val = val + 2.0*eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, UU)
            ! + 3 uzzzz(x,y,z,t)
            val = val + 3.0*eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, UU)
            ! + 4 utttt(x,y,z,t)
            val = val + 4.0*eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, UU)
            ! + 5  uxyt(x,y,z,t)
            val = val + 5.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, UU)
            ! + 6 uyyzz(x,y,z,t)
            val = val + 6.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, UU)
            ! + 7  uztt(x,y,z,t)
            val = val + 7.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, UU)
            ! + 8  uxxx(x,y,z,t)
            val = val + 8.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, UU)
            ! + 9  uyyt(x,y,z,t)
            val = val + 9.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, UU)
            ! +10   uzt(x,y,z,t)
            val = val + 10.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, UU)
            ! +11    ut(x,y,z,t)
            val = val + 11.0*eval_derivative(0, 0, 0, 1, i, ii, iii, iiii, UU)
            ! +12     u(x,y,z,t)
            val = val + 12.0*uana(xg(i),yg(ii),zg(iii),tg(iiii))
            ! -       f(x,y,z,t) should be zero
            err = abs(val - FF(xg(i),yg(ii),zg(iii),tg(iiii)))
            if(err>smaxerr) then
              smaxerr = err
            end if
          end do !iiii
        end do !iii
      end do ! ii
    end do ! i
    print *, 'check_soln maxerr=',smaxerr
  end subroutine check_soln

program pde44_sparse
  use global
  use formain
  use uses
  implicit none
  integer :: i, j, ii, jj, iii, jjj, iiii, jjjj
  integer :: k
  double precision :: x, y, z, t
  double precision :: val, err, avgerr, maxerr

  print *, 'pde44_sparse.f90 running'
  print *,' solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +'
  print *,'     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +'
  print *,'     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +'
  print *,'     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) ='
  print *,'     F(x,y,z,t)'
  print *,'     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +'
  print *,'                  180 z^2*t + 200 y^2z*t + 44 t^3 +'
  print *,'                  110 y^2z^2t + 66 xy + 12t^4 +'
  print *,'                  24 z^4 + 36 y^4 + 48 x^4 +'
  print *,'                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t' 
  print *,' '
  print *,' boundary conditions computed using u(x,y,z,t)'
  print *,' analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +'
  print *,'                                   5 y^2 z^2 t^2 + 6 x y t + '
  print *,'                                   7 x z + 8 t + 9'
  print *, ' '
  print *, 'xmin=',xmin,', xmax =',xmax
  print *, 'ymin=',ymin,', ymax =',ymax
  print *, 'zmin=',zmin,', zmax =',zmax
  print *, 'tmin=',tmin,', tmax =',tmax
  print *, 'nx=',nx,', ny =',ny,', nz =',nz,', nt=',nt 

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymin,zmin,tmin'
  do i=1,nx 
    xg(i) = xmin+(i-1)*hx
    Ua(i) = uana(xg(i),ymin,zmin,tmin) ! 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'
  do ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    Ua(ii) = uana(xmin,yg(ii),zmin,tmin)
    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'
  do iii=1,nz 
    zg(iii) = zmin+(iii-1)*hz
    Ua(iii) = uana(xmin,ymin,zg(iii),tmin)
    print *, 'iii=',iii,', Ua(',zg(iii),')=',Ua(iii)
  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'
  do iiii=1,nt 
    tg(iiii) = tmin+(iiii-1)*ht
    Ua(iiii) = uana(xmin,ymin,zmin,tg(iiii))
    print *, 'iiii=',iiii,', Ua(',tg(iiii),')=',Ua(iiii)
  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 iii=2,nz-1
        z = zg(iii)
        do iiii=2,nt-1
          t = tg(iiii)
          Ua(s(i,ii,iii,iiii)) = uana(x,y,z,t)
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i

  ! initialize A and R  inside boundary
  call init(A,nxyzt)

  do i=2,nx-1 ! could go directly into A
    do ii=2,ny-1 
      do iii=2,nz-1 
        do iiii=2,nt-1 
          R(s(i,ii,iii,iiii)) = 0.0
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i

  ! compute entries in stiffness matrix
  call buildA()

  print *, 'A and R computed'
  call putrhs(A,R)
  call write_all(A) ! optional
  print *, ' '

  call simeq(A,U)

  call check_soln(Ua)
  call check_soln(U)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'U computed, Ua analytic, error'
  do i=2,nx-1 
    do ii=2,ny-1 
      do iii=2,nz-1 
        do iiii=2,nt-1 
          err = abs(U(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
          if(err>maxerr) then
            maxerr = err
          end if
          avgerr = avgerr + err
          print "('U(',I1,',',I1,','I1,',',I1')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
              i, ii, iii, iiii, U(s(i,ii,iii,iiii)), Ua(s(i,ii,iii,iiii)), &
              (U(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nxyzt)
  print *, 'end pde44_sparse.f90 '

end program pde44_sparse