-- pde44h_eq.adb  homogenous harmonic equation
-- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
--       utttt(x,y,z,t) + 2  uxx(x,y,z,t) + 2 uyy(x,y,z,t) +
--       2 uzz(x,y,z,t) + 2  u(x,y,z,t) = 0
--
-- boundary conditions computed using u(x,y,z,t)
-- analytic solution is u(x,y,z,t) = sin(x+y+z+t)
--
-- replace continuous derivatives with discrete derivatives
-- solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
-- A * U = R, solve simultaneous equations for U


-- fourth order numerical difference rderiv order=4, npoints=5, term=1
-- d^4u/dx^4 = cxxxx(1)*u(x,y,z,t) + cxxxx(2)*u(x+hx,y,z,t) +
-- cxxxx(3)*u(x+2hx,y,z,t) + cxxxx(4)*u(x+3hx,y,z,t) + cxxxx(5)*u(x+4hx,y,z,t))
--
-- d^4u/dx^4 using subscripts i,j,k,l for i a minimum subscript of x
-- d^4u/dx^4 = cxxxx(1)*u(i,j,k,l) + cxxxx(2)*u(i+1,j,k,l) +
-- cxxxx(3)*u(i+2,j,k,l) + cxxxx(4)*u(i+3,j,k,l) + cxxxx(5)*u(i+4,j,k,l))
--
-- d^4u/dy^4 using subscripts i,j,k,l for j a minimum subscript of y, etc
-- d^4u/dy^4 = cyyyy(1)*u(i,j,k,l) + cyyyy(2)*u(i,j+1,k,l) +
-- cyyyy(3)*u(i,j+2,k,l) + cyyyy(4)*u(i,j+3,k,l) + cyyyy(5)*u(i,j+4,k,l))
--
-- subscripts in code use xg(i), yg(ii), zg(iii), tg(iiii)
-- derivatives are computed with hx, order, point included as in
-- cxx(j) for second derivative of x at point j
-- cxxxx(j) for fourth derivative of x at point j

-- The subscript versions of the derivatives are substituted into the
-- equation to be solved. The resulting equation is numerically solved
-- for u(i,j,k,l)
-- The boundaries xmin..xmax, ymin..ymax, zmin..zmax, tmin..tmax are known
-- The number of grid points in each dimension is known nx, ny, nz, nt
-- hx=(xmax-xmin)/(nx-1), hy=(ymax-ymin)/(ny-1),
-- hz=(zmax-zmin)/(nz-1), ht=(tmax-tmin)/(tx-1),
-- thus xi=xmin+i*hx, yj=ymin+j*hy, zk=zmin+k*hz, tl=tmin+l*ht
-- then the linear system of equations is built,
-- nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2) equations because the boundary value
--                                   is known on each end.
-- The matrix A is nxyzt rows by nxyzt columns (nxyzt+1) columns including F
-- The forcing function vector F, zero here, is nxyzt elements adjunct to A
-- The solution vector U is nxyzt elements.
-- The building of the A matrix requires 4 nested loops for rows
-- i  in 2..nx-1, ii in 2..ny-1, iii in 2..nz-1, iiii in 2..nt-1 for rows of A
-- then each term in the differential equation fills in that row
-- and cs = (nx-2)*(ny-2)*(nz-2)*(nt-2)+1 the RHS, fills the last column of A
-- subscripting functions are used to compute linear subscripts of A, U
-- row=(i-1)*(ny-2)*(nz-2)*(nt-2)+(ii-1)*(nz-2)*(nt-2)+(iii-1)*(nt-2)+(iiii-1)
-- row*(nxyzt+1)+col for A
--
-- care must be taken to move the negative of boundary elements to
-- the right hand size of the equation to be solved. These are subtracted
-- from the computed value of the forcing function, cs, for that row.
-- tests are j=1 or j=nx   ...  jjjj=1 or jjjj=nt for boundaries

with ada.text_io;
use  ada.text_io;
with ada.numerics;
use  ada.numerics;
with ada.numerics.Long_elementary_Functions;
use  ada.numerics.Long_elementary_Functions;
with real_arrays;
use  real_arrays; -- defines type real
with simeq;
with rderiv;
with ada.calendar;
use  ada.Calendar;

procedure pde44h_eq is
  nx : integer := 7; -- problem parameters
  ny : integer := 7;
  nz : integer := 7;
  nt : integer := 7;
  nxyzt : integer := (nx-2)*(ny-2)*(nz-2)*(nt-2);

  A: real_matrix(1..nxyzt,1..nxyzt); -- A * U = R
  R : real_vector(1..nxyzt);  -- RHS known
  U : real_vector(1..nxyzt);  -- solution being computed
  Ua : real_vector(1..nxyzt); -- solution for checking
  xg : real_vector(1..nx);
  yg : real_vector(1..ny);
  zg : real_vector(1..nz);
  tg : real_vector(1..nt);
  xmin : real := 0.0; -- problem parameters
  xmax : real := 1.0;
  ymin : real := 0.0;
  ymax : real := 1.0;
  zmin : real := 0.0;
  zmax : real := 1.0;
  tmin : real := 0.0;
  tmax : real := 1.0;
  hx, hy, hz, ht : real;

  package real_io is new float_io(real);
  use real_io;

  -- internal functions

  function max(x:integer;y:integer) return integer is
  begin
     if x>y then return x; end if;
     return y;
  end max;

  function s(i:integer; ii:integer; iii:integer; iiii:integer)
            return integer is -- for internal x,y,z,t
  begin
    return (i-2)*(ny-2)*(nz-2)*(nt-2) + (ii-2)*(nz-2)*(nt-2) +
           (iii-2)*(nt-2) + (iiii-1);
  end S;

  -- problem RHS, 0

  -- problem boundaries (and analytic solution for checking)
  function uana(x:real; y:real; z:real; t:real)
           return real is
  begin
    return sin(x+y+z+t);
  end uana;

  procedure printA is
     k : integer; -- no cs, R(k)
  begin
    for i in 2..nx-1 loop
      for ii in 2..ny-1 loop
        for iii in 2..nz-1 loop
          for iiii in 2..nt-1 loop
            k := s(i,ii,iii,iiii); -- row index
            for j in 2..nx-1 loop
              for jj in 2..ny-1 loop
                for jjj in 2..nz-1 loop
                  for jjjj in 2..nt-1 loop
                     put_line("A(row("&integer'image(i)&","&integer'image(ii)&
                              ","&integer'image(iii)&","&integer'image(iiii)&
                              "),col("&integer'image(j)&","&integer'image(jj)&
                              ","&integer'image(jjj)&","&integer'image(jjjj)&
                              ")) ="&
                              real'image(A(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj))));
                  end loop; -- jjjj
                end loop; -- jjj
              end loop; -- jj
            end loop; -- j
            put_line("RHS R("&integer'image(k)&")="&real'image(R(k)));
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    put_line(" ");
  end printA;


  function eval_derivative(xord:integer; yord:integer; zord:integer; tord:integer;
                           i:integer; ii:integer; iii:integer; iiii:integer;
                           u:real_vector) return real is
    cx : real_vector(1..nx); -- can be used for any order
    cy : real_vector(1..ny);
    cz : real_vector(1..nz);
    ct : real_vector(1..nt);
    nn : integer := max(nx, max(ny, max(nz, nt)));
    p1 : real_vector(1..nn); -- another dimension of any order
    p2 : real_vector(1..nn);
    p3 : real_vector(1..nn);
    discrete : real;
    x, y, z, t : real;
  begin
    rderiv(xord, nx, i, hx, cx);
    x := xg(i);
    rderiv(yord, ny, ii, hy, cy);
    y := yg(ii);
    rderiv(zord, nz, iii, hz, cz);
    z := zg(iii);
    rderiv(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/=0 and yord=0 and zord=0 and tord=0 then
      for j in 1..nx loop
        if j=1 or j=nx then
          discrete := discrete + cx(j)*uana(xg(j),y,z,t);
        else
          discrete := discrete + cx(j)*u(s(j,ii,iii,iiii));
        end if;
      end loop; -- j
    elsif xord=0 and yord/=0 and zord=0 and tord=0 then
      for jj in 1..ny loop
        if jj=1 or jj=ny then
          discrete := discrete + cy(jj)*uana(x,yg(jj),z,t);
        else
          discrete := discrete + cy(jj)*u(s(i,jj,iii,iiii));
        end if;
      end loop; -- jj
    elsif xord=0 and yord=0 and zord/=0 and tord=0 then
      for jjj in 1..nz loop
        if jjj=1 or jjj=nz then
          discrete := discrete + cz(jjj)*uana(x,y,zg(jjj),t);
        else
          discrete := discrete + cz(jjj)*u(s(i,ii,jjj,iiii));
        end if;
      end loop; -- jjj
    elsif xord=0 and yord=0 and zord=0 and tord/=0 then
      for jjjj in 1..nt loop
        if jjjj=1 or jjjj=nt then
          discrete := discrete + ct(jjjj)*uana(x,y,z,tg(jjjj));
        else
           discrete := discrete + ct(jjjj)*u(s(i,ii,iii,jjjj));
        end if;
      end loop; -- jjjj

    -- six cases of two partials xy, xz, xt, yz, yt, zt to any order
    elsif xord/=0 and yord/=0 and zord=0 and tord=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jj in 1..ny loop
          if j=1 or j=nx or jj=1 or jj=ny then
            p1(j) := p1(j) + cy(jj)*uana(xg(j),yg(jj),zg(iii),tg(iiii));
          else
             p1(j) := p1(j) + cy(jj)*u(s(j,jj,iii,iiii));
          end if;
        end loop; -- jj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    elsif xord/=0 and yord=0 and zord/=0 and tord=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jjj in 1..nz loop
          if j=1 or j=nx or jjj=1 or jjj=nz then
            p1(j) := p1(j) + ct(jjj)*uana(xg(j),yg(ii),zg(jjj),tg(iiii));
          else
            p1(j) := p1(j) + cz(jjj)*u(s(j,ii,jjj,iiii));
          end if;
        end loop; -- jjj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    elsif xord/=0 and yord=0 and zord=0 and tord/=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jjjj in 1..nt loop
          if j=1 or j=nx or jjjj=1 or jjjj=nt then
            p1(j) := p1(j) + ct(jjjj)*uana(xg(j),yg(ii),zg(iii),tg(jjjj));
          else
             p1(j) := p1(j) + ct(jjjj)*u(s(j,ii,iii,jjjj));
          end if;
        end loop; -- jjjj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    elsif xord=0 and yord/=0 and zord/=0 and tord=0 then
      for jj in 1..ny loop
        p1(jj) := 0.0;
        for jjj in 1..nz loop
          if jj=1 or jj=ny or jjj=1 or jjj=nz then
            p1(jj) := p1(jj) + cz(jjj)*uana(xg(i),yg(jj),zg(jjj),tg(iiii));
          else
             p1(jj) := p1(jj) + cz(jjj)*u(s(i,jj,jjj,iiii));
          end if;
        end loop; -- jjj
        discrete := discrete + cy(jj)*p1(jj);
      end loop; -- jj
    elsif xord=0 and yord/=0 and zord=0 and tord/=0 then
      for jj in 1..ny loop
        p1(jj) := 0.0;
        for jjjj in 1..nt loop
          if jj=1 or jj=ny or jjjj=1 or jjjj=nt then
            p1(jj) := p1(jj) + ct(jjjj)*uana(xg(i),yg(jj),zg(iii),tg(jjjj));
          else
             p1(jj) := p1(jj) + ct(jjjj)*u(s(i,jj,iii,jjjj));
          end if;
        end loop; -- jjjj
        discrete := discrete + cy(jj)*p1(jj);
      end loop; -- jj
    elsif xord=0 and yord=0 and zord/=0 and tord/=0 then
      for jjj in 1..nz loop
        p1(jjj) := 0.0;
        for jjjj in 1..nt loop
          if jjj=1 or jjj=nz or jjjj=1 or jjjj=nt then
            p1(jjj) := p1(jjj) + ct(jjjj)*uana(xg(i),yg(ii),zg(jjj),tg(jjjj));
          else
             p1(jjj) := p1(jjj) + ct(jjjj)*u(s(i,ii,jjj,jjjj));
          end if;
        end loop; -- jjjj
        discrete := discrete + cz(jjj)*p1(jjj);
      end loop; -- jjj

    -- four cases of three partials xyz, xyt, xzt, yzt to any order
    elsif xord/=0 and yord/=0 and zord/=0 and tord=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jj in 1..ny loop
          p2(jj) := 0.0;
          for jjj in 1..nz loop
            if j=1 or j=nx or
               jj=1 or jj=ny or
               jjj=1 or jjj=nz then
              p2(jj) := p2(jj) + cz(jjj)*uana(xg(j),yg(jj),zg(jjj),tg(iiii));
            else
               p2(jj) := p2(jj) + cz(jjj)*u(s(j,jj,jjj,iiii));
            end if;
          end loop; -- jjj
          p1(j) := p1(j) + cy(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    elsif xord/=0 and yord/=0 and zord=0 and tord/=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jj in 1..ny loop
          p2(jj) := 0.0;
          for jjjj in 1..nt loop
            if j=1 or j=nx or
               jj=1 or jj=ny or
               jjjj=1 or jjjj=nt then
              p2(jj) := p2(jj) + ct(jjjj)*uana(xg(j),yg(jj),zg(iii),tg(jjjj));
            else
              p2(jj) := p2(jj) + ct(jjjj)*u(s(j,jj,iii,jjjj));
            end if;
          end loop; -- jjjj
          p1(j) := p1(j) + cy(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete +cx(j)*p1(j);
      end loop; -- j
    elsif xord/=0 and yord=0 and zord/=0 and tord/=0 then
      for j in 1..nx loop
        p1(j) := 0.0;
        for jjj in 1..nz loop
          p2(jjj) := 0.0;
          for jjjj in 1..nt loop
            if j=1 or j=nx or
               jjj=1 or jjj=nz or
               jjjj=1 or jjjj=nt then
              p2(jjj) := p2(jjj) + ct(jjjj)*uana(xg(j),yg(ii),zg(jjj),tg(jjjj));
            else
              p2(jjj) := p2(jjj) + ct(jjjj)*u(s(j,ii,jjj,jjjj));
            end if;
          end loop; -- jjjj
          p1(j) := p1(j) + cz(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    elsif xord=0 and yord/=0 and zord/=0 and tord/=0 then
      for jj in 1..ny loop
        p1(jj) := 0.0;
        for jjj in 1..nz loop
          p2(jjj) := 0.0;
          for jjjj in 1..nt loop
            if jj=1 or jj=ny or
               jjj=1 or jjj=nz or
               jjjj=1 or jjjj=nt then
              p2(jjj) := p2(jjj) + ct(jjjj)*uana(xg(i),yg(jj),zg(jjj),tg(jjjj));
            else
               p2(jjj) := p2(jjj) + ct(jjjj)*u(s(i,jj,jjj,jjjj));
            end if;
          end loop; -- jjjj
          p1(jj) := p1(jj) + cz(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cy(jj)*p1(jj);
      end loop; -- jj

    -- one case of four partials
    else
      for j in 1..nx loop
        p1(j) := 0.0;
        for jj in 1..ny loop
          p2(jj) := 0.0;
          for jjj in 1..nz loop
            p3(jjj) := 0.0;
            for jjjj in 1..nt loop
              if   j=1 or j=nx or
                   jj=1 or jj=ny or
                   jjj=1 or jjj=nz or
                   jjjj=1 or jjjj=nt then
                p3(jjj) := p3(jjj) + ct(jjjj)*uana(xg(j),yg(jj),zg(jjj),tg(jjjj));
              else
                p3(jjj) :=  p3(jjj) + ct(jjjj)*u(s(j,jj,jjj,jjjj));
              end if;
            end loop; -- jjjj
            p2(jj) := p2(jj) + cz(jjj)*p3(jjj);
          end loop; -- jjj
          p1(j) := p1(j) + cy(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cx(j)*p1(j);
      end loop; -- j
    end if;
    return discrete;
  end eval_derivative;

  procedure check_soln(U:real_vector) is
    -- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    --       utttt(x,y,z,t) + 2  uxx(x,y,z,t) + 2 uyy(x,y,z,t) +
    --       uzz(x,y,z,t) + 2 u(x,y,z,t) = 0.0
    val, err, smaxerr : real;
  begin
    smaxerr := 0.0;
    for i in 2..nx-1 loop
      for ii in 2..ny-1 loop
        for iii in 2..nz-1 loop
          for iiii in 2..nt-1 loop
            val := 0.0;
            -- uxxxx(x,y,z,t)
            val := val + eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u);
            -- + uyyyy(x,y,z,t)
            val := val + eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u);
            -- + uzzzz(x,y,z,t)
            val := val + eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u);
            -- + utttt(x,y,z,t)
            val := val + eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u);
            -- + 2 uxx(x,y,z,t)
            val := val + 2.0*eval_derivative(2, 0, 0, 0, i, ii, iii, iiii, u);
            -- + 2 uyy(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 2, 0, 0, i, ii, iii, iiii, u);
            -- + 2 uzz(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 0, 2, 0, i, ii, iii, iiii, u);
            -- + 2 u(x,y,z,t)
            val := val + 2.0*uana(xg(i),yg(ii),zg(iii),tg(iiii));
            --  should be zero
            err := abs(val);
            if err>smaxerr then
              smaxerr := err;
            end if;
          end loop; --iiii
        end loop; --iii
      end loop; -- ii
    end loop; -- i
    put_line("check_soln maxerr="&real'image(smaxerr));
  end Check_Soln;

  procedure Write_Soln(Which:Integer; Tag:String) is
    Fout : File_Type;
  begin
    Create(Fout, Out_File, "pde44h_eq_ada"&Tag&".dat");

    Put_Line("writing pde44h_eq_ada"&Tag&".dat");
    Put_Line(Fout,"4 "&integer'image(nx)&integer'image(ny)&
                       integer'image(nz)&integer'image(nt));
    for  i in 1..nx loop
      for ii in 1..ny loop
        for  iii in 1..nz loop
          for iiii in 1..nt loop
            Put(Fout," ");
            Put(Fout,xg(i),2,6,0);
            Put(Fout," ");
            Put(Fout,yg(ii),2,6,0);
            Put(Fout," ");
            Put(Fout,zg(iii),2,6,0);
            Put(Fout," ");
            Put(Fout,tg(iiii),2,6,0);
            Put(Fout," ");
            if Which=1 then -- expected solution
               Put(Fout,uana(xg(i),yg(ii),zg(iii),tg(iiii)),4,6,0);
            else -- computed solution
              if i=1 or i=nx or ii=1 or ii=ny or iii=1 or
                 iii=nz or iiii=1 or iiii=nt then -- boundary
                Put(Fout,uana(xg(i),yg(ii),zg(iii),tg(iiii)),4,6,0);
              else
                Put(Fout,U(S(i,ii,iii,iiii)),4,6,0);
              end if;
            end if;
            New_Line(Fout);
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    Close(Fout);
    Put_Line("finished writing pde44h_eq_ada"&Tag&".dat");
  end Write_Soln;

  procedure initializeA is
  begin
    put_line("initialize A ");
    -- initialize A
    for i in 2..nx-1 loop
      for ii in 2..ny-1 loop
        for iii in 2..nz-1 loop
          for iiii in 2..nt-1 loop
            for j in 2..nx-1 loop
              for jj in 2..ny-1 loop
                for jjj in 2..nz-1 loop
                  for jjjj in 2..nt-1 loop
                    A(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) := 0.0;
                  end loop; -- jjjj
                end loop; -- jjj
              end loop; -- jj
            end loop; -- j
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
  end initializeA;


  procedure buildA is -- system of equations
    val : real := 0.0;
    k: integer;
    cxx, cxxxx : real_vector(1..nx);
    cyy, cyyyy : real_vector(1..ny);
    czz, czzzz : real_vector(1..nz);
    ctttt : real_vector(1..nt);
  begin
    -- compute entries in A matrix, inside boundary
    put_line("compute A matrix and Y RHS ");


    put_line("check rderiv(4,nx,1,hx,cxxxx");
    rderiv(4,nx,1,hx,cxxxx);
    for i in 1..nx loop put_line("cxxxx("&integer'image(i)&
                                 ")="&real'image(cxxxx(i))); end loop;
    put_line("check rderiv(4,nx,nx,hx,cxxxx");
    rderiv(4,nx,nx,hx,cxxxx);
    for i in 1..nx loop put_line("cxxxx("&integer'image(i)&
                                 ")="&real'image(cxxxx(i))); end loop;
    put_line(" ");

    val := 0.0;
    for i in 2..nx-1 loop
      for ii in 2..ny-1 loop
        for iii in 2..nz-1 loop
          for iiii in 2..nt-1 loop
            k := s(i,ii,iii,iiii); -- row index

            -- for each term in first equation

            -- for d^4u/dx^4
            rderiv(4,nx,i,hx,cxxxx);
            for j in 1..nx loop
              val := cxxxx(j);
              if j=1 or j=nx then
                R(k) := R(k) - val*uana(xg(j),yg(ii),zg(iii),tg(iiii));
              else
                 A(k,s(j,ii,iii,iiii)) := A(k,s(j,ii,iii,iiii)) + val;
              end if;
            end loop; -- j

            -- for d^4u/dy^4
            rderiv(4,ny,ii,hy,cyyyy);
            for jj in 1..ny loop
              val := cyyyy(jj);
              if jj=1 or jj=ny then
                R(k) := R(k) - val*uana(xg(i),yg(jj),zg(iii),tg(iiii));
              else
                A(k,s(i,jj,iii,iiii)) := A(k,s(i,jj,iii,iiii)) + val;
              end if;
            end loop; -- jj

            -- for d^4u/dz^4
            rderiv(4,nz,iii,hz,czzzz);
            for jjj in 1..nz loop
              val := czzzz(jjj);
              if jjj=1 or jjj=nz then
                R(k) := R(k) - val*uana(xg(i),yg(ii),zg(jjj),tg(iiii));
              else
                A(k,s(i,ii,jjj,iiii)) := A(k,s(i,ii,jjj,iiii)) + val;
              end if;
            end loop; -- jjj

            -- for d^4u/dt^4
            rderiv(4,nt,iiii,ht,ctttt);
            for jjjj in 1..nt loop
              val := ctttt(jjjj);
              if jjjj=1 or jjjj=nt then
                R(k) := R(k) - val*uana(xg(i),yg(ii),zg(iii),tg(jjjj));
              else
                A(k,s(i,ii,iii,jjjj)) := A(k,s(i,ii,iii,jjjj)) + val;
              end if;
            end loop; -- jjjj

            -- for d^2u/dx^2
            rderiv(2,nx,i,hx,cxx);
            for j in 1..nx loop
              val := 2.0*cxx(j);
              if j=1 or j=nx then
                R(k) := R(k) - val*uana(xg(j),yg(ii),zg(iii),tg(iiii));
              else
                 A(k,s(j,ii,iii,iiii)) := A(k,s(j,ii,iii,iiii)) + val;
              end if;
            end loop; -- j

            -- for d^2u/dy^2
            rderiv(2,ny,ii,hy,cyy);
            for jj in 1..ny loop
              val := 2.0*cyy(jj);
              if jj=1 or jj=ny then
                R(k) := R(k) - val*uana(xg(i),yg(jj),zg(iii),tg(iiii));
              else
                A(k,s(i,jj,iii,iiii)) := A(k,s(i,jj,iii,iiii)) + val;
              end if;
            end loop; -- jj

            -- for d^2u/dz^2
            rderiv(2,nz,iii,hz,czz);
            for jjj in 1..nz loop
              val := 2.0*czz(jjj);
              if jjj=1 or jjj=nz then
                R(k) := R(k) - val*uana(xg(i),yg(ii),zg(jjj),tg(iiii));
              else
                A(k,s(i,ii,jjj,iiii)) := A(k,s(i,ii,jjj,iiii)) + val;
              end if;
            end loop; -- jjj

            -- for u(x,y,z,t)
            A(k,s(i,ii,iii,iiii)) := A(k,s(i,ii,iii,iiii)) + 2.0;

            -- for f(x,y,z,t) RHS constant, zero here
            -- R(k) := R(k) + f(xg(i),yg(ii),zg(iii),tg(iiii));

            -- end terms

          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii loop
    end loop; -- i
  end BuildA;

  procedure main is
    x, y, z, t : real;
    err, avgerr, maxerr : real;
    time_start : duration;
    time_now : duration;
    time_last : duration;

  begin -- main
    put_line("pde44h_eq.adb running ");
    put_line("Given uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + ");
    put_line("      utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) + ");
    put_line("      2 uzz(x,y,z,t) + 2 u(x,y,z,t) = 0");
    put_line("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    put_line("Analytic solution u(x,y,z,t) =  sin(x+y+z+t)");
    put_line(" ");

    put_line("xmin="&real'image(xmin)&", xmax="&real'image(xmax));
    put_line("ymin="&real'image(ymin)&", ymax="&real'image(ymax));
    put_line("zmin="&real'image(zmin)&", zmax="&real'image(zmax));
    put_line("tmin="&real'image(tmin)&", tmax="&real'image(tmax));
    put_line("nx="&integer'image(nx)&", ny="&integer'image(ny)&
             ", nz="&integer'image(nz)&", nt="&integer'image(nt));
    time_start := seconds(clock);
    time_last := time_start;

    hx := (xmax-xmin)/real(nx-1);
    put_line("x grid, hx="&real'image(hx));
    for i in 1..nx loop
      xg(i) := xmin+real(i-1)*hx;
      put_line("xg("&integer'image(i)&")="&real'image(xg(i)));
    end loop; -- i
    put_line(" ");
    hy := (ymax-ymin)/real(ny-1);
    put_line("y grid, hy="&real'image(hy));
    for ii in 1..ny loop
      yg(ii) := ymin+real(Ii-1)*hy;
      put_line("yg("&integer'image(ii)&")="&real'image(yg(ii)));
    end loop; -- ii
    put_line(" ");
    hz := (zmax-zmin)/real(nz-1);
    put_line("z grid, hz="&real'image(hz));
    for iii in 1..nz loop
      zg(iii) := zmin+real(Iii-1)*hz;
      put_line("zg("&integer'image(iii)&")="&real'image(zg(iii)));
    end loop; -- iii
    put_line(" ");
    ht := (tmax-tmin)/real(nt-1);
    put_line("t grid, ht="&real'image(ht));
    for iiii in 1..nt loop
      tg(iiii) := tmin+real(Iiii-1)*ht;
      put_line("tg("&integer'image(iiii)&")="&real'image(tg(iiii)));
    end loop; -- iiii
    put_line(" ");

    put_line("put solution in Ua vector ");
    -- put solution in Ua vector
    for i in 2..nx-1 loop
      x := xg(i);
      for ii in 2..ny-1 loop
        y := yg(ii);
        for iii in 2..nz-1 loop
          z := zg(iii);
          for iiii in 2..nt-1 loop
            t := tg(iiii);
            Ua(s(i,ii,iii,iiii)) := uana(x,y,z,t);
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    time_now := seconds(clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    initializeA;

    time_now := seconds(clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    buildA; -- system of equations

    -- printA; -- optional output

    put_line("A computed stiffness matrix ");
    time_now := seconds(clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");


    put_line("solve "&integer'image(nxyzt)&" equations in "&
             integer'image(nxyzt)&" unknowns ");

    U := simeq(A, R);

    put_line("equations solved ");
    time_now := seconds(clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    put_line("check pde against known solution");
    check_soln(Ua);
    put_line("check pde against computed solution");
    check_soln(U);
    put_line(" ");

    avgerr := 0.0;
    maxerr := 0.0;
    put_line("          U computed, Ua analytic, error ");
    for i in 2..nx-1 loop
      for ii in 2..ny-1 loop
        for iii in 2..nz-1 loop
          for iiii in 2..nt-1 loop
            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;
            put("ug("&integer'image(i)&","&integer'image(ii)&
                ","&integer'image(iii)&","&integer'image(iiii)&")=");
            put(U(s(i,ii,iii,iiii)),3,6); put(", Ua=");
            put(Ua(s(i,ii,iii,iiii)),3,6); put(", err=");
            put(err);
            new_line;
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i

    time_now := seconds(clock);
    put_line("total CPU time  = "&Duration'image(time_now-time_start)&
             " seconds");
    put_line(" ");


    put_line("                     maxerr="&real'image(maxerr)&
             ", avgerr="&real'image(avgerr/real(nx*ny*nz*nt)));
    put_line(" ");
  end Main;
begin -- pde44h_eq
   main;
   Write_Soln(1,"");
   Write_Soln(2,"1");
   put_line("end first test");
   new_line;
   xmax := 1.57;
   ymax := 1.57;
   zmax := 1.57;
   tmax := 1.57;
   main;
   Write_Soln(2,"2");
   put_line("end second test");
   new_line;
   xmin := -1.0;
   ymin := -1.0;
   zmin := -1.0;
   tmin := -1.0;
   xmax := 1.0;
   ymax := 1.0;
   zmax := 1.0;
   tmax := 1.0;
   main;
   Write_Soln(2,"3");
   put_line("end pde44h_eq.adb");
end pde44h_eq;