-- pde44_sparse.adb  discretize, build linear equations, 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
--
-- 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
-- cx(j) for first derivative of x at point j
-- cxxxx(j) for fourth derivative of x at point j
-- cyy(jj)*czz(jjj) for d^4u/dy^2dz^2 at jj, jjj

-- The subscript versions of the derivatives are substituted into the
-- equation to be solved. The resulting equation is analytically solved
-- for u(i,j,k,l) = other u terms with coefficients and + f(xi,yj,zk,tl)
-- 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 is nxyzt elements adjunct to A
-- The solution vector U is nxyzt elements.
-- The building of the A matix 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 sparse1; -- 1..Nrow
with rderiv;
with calendar;

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

  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;
  coefdxxxx : real; -- computed below if not constant
  coefdyyyy : real;
  coefdzzzz : real;
  coefdtttt : real;
  coefdxyt  : real;
  coefdyyzz : real;
  coefdztt  : real;
  coefdxxx  : real;
  coefdyyt  : real;
  coefdzt   : real;
  coefdt    : real;
  coefu     : real;

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

  Package A is new Sparse1(Nxyzt);

  -- 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
  function f(x:real; y:real; z:real; t:real) return real is
  begin
    return 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 f;

  -- problem boundaries (and analytic solution for checking)
  function uana(x:real; y:real; z:real; t:real)
           return real is
  begin
    return 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 uana;


  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);
    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);
    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) + 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)
    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);
            -- + 2 uyyyy(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u);
            -- + 3 uzzzz(x,y,z,t)
            val := val + 3.0*eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u);
            -- + 4 utttt(x,y,z,t)
            val := val + 4.0*eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u);
            -- + 5  uxyt(x,y,z,t)
            val := val + 5.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, u);
            -- + 6 uyyzz(x,y,z,t)
            val := val + 6.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, u);
            -- + 7  uztt(x,y,z,t)
            val := val + 7.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, u);
            -- + 8  uxxx(x,y,z,t)
            val := val + 8.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, u);
            -- + 9  uyyt(x,y,z,t)
            val := val + 9.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, u);
            -- +10   uzt(x,y,z,t)
            val := val + 10.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, u);
            -- +11    ut(x,y,z,t)
            val := val + 11.0*eval_derivative(0, 0, 0, 1, i, ii, iii, iiii, u);
            -- +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 - f(xg(i),yg(ii),zg(iii),tg(iiii)));
            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 buildA is -- system of equations
    val : real := 0.0;
    k, cs : integer;
    cx,      cxxx, cxxxx : real_vector(1..nx);
    cy, cyy,       cyyyy : real_vector(1..ny);
    cz, czz,       czzzz : real_vector(1..nz);
    ct, ctt,       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(" ");
    cs := Nxyzt+1; -- RHS
    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
            coefdxxxx := 1.0;
            rderiv(4,nx,i,hx,cxxxx);
            for j in 1..nx loop
              val := coefdxxxx*cxxxx(j);
              if j=1 or j=nx then
                A.add(k,cs,-val*uana(xg(j),yg(ii),zg(iii),tg(iiii)));
              else
                 A.add(k,s(j,ii,iii,iiii),Val);
              end if;
            end loop; -- j

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

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

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

            -- for d^3u/dxdydt
            coefdxyt  := 5.0;
            rderiv(1,nx,i,hx,cx);
            rderiv(1,ny,ii,hy,cy);
            rderiv(1,nt,iiii,ht,ct);
            for j in 1..nx loop
              for jj in 1..ny loop
                for jjjj in 1..nt loop
                  val := coefdxyt*cx(j)*cy(jj)*ct(jjjj);
                  if j=1 or j=nx or jj=1 or jj=ny or
                     jjjj=1 or jjjj=nt then
                    A.add(k,cs,-val*uana(xg(j),yg(jj),zg(iii),tg(jjjj)));
                  else
                    A.add(k,s(j,jj,iii,jjjj),Val);
                  end if;
                end loop; -- jjjj
              end loop; -- jj
            end loop; -- j

            -- for d^4u/dy^2dz^2
            coefdyyzz := 6.0;
            rderiv(2,ny,ii,hy,cyy);
            rderiv(2,nz,iii,hz,czz);
            for jj in 1..ny loop
              for jjj in 1..nz loop
                val := coefdyyzz*cyy(jj)*czz(jjj);
                if jj=1 or jj=ny or jjj=1 or jjj=nz then
                  A.add(k,cs,-val*uana(xg(i),yg(jj),zg(jjj),tg(iiii)));
                else
                  A.add(k,s(i,jj,jjj,iiii),Val);
                end if;
              end loop; -- jjj
            end loop; -- jj

            -- for d^3u/dzdt^2
            coefdztt  := 7.0;
            rderiv(1,nz,iii,hz,cz);
            rderiv(2,nt,iiii,ht,ctt);
            for jjj in 1..nz loop
              for jjjj in 1..nt loop
                val := coefdztt*cz(jjj)*ctt(jjjj);
                if jjj=1 or jjj=nz or jjjj=1 or jjjj=nt then
                  A.add(k,cs,-val*uana(xg(i),yg(ii),zg(jjj),tg(jjjj)));
                else
                  A.add(k,s(i,ii,jjj,jjjj),Val);
                end if;
              end loop; -- jjjj
            end loop; -- jjj

            -- for d^3u/dx^3
            coefdxxx  := 8.0;
            rderiv(3,nx,i,hx,cxxx);
            for j in 1..nx loop
              val := coefdxxx*cxxx(j);
              if j=1 or j=nx then
                A.add(k,cs,-val*uana(xg(j),yg(ii),zg(iii),tg(iiii)));
              else
                A.add(k,s(j,ii,iii,iiii),Val);
              end if;
            end loop; -- j

            -- for d^3u/dy^2dt
            coefdyyt  := 9.0;
            rderiv(2,ny,ii,hy,cyy);
            rderiv(1,nt,iiii,ht,ct);
            for jj in 1..ny loop
              for jjjj in 1..nt loop
                val := coefdyyt*cyy(jj)*ct(jjjj);
                if jj=1 or jj=ny or jjjj=1 or jjjj=nt then
                  A.add(k,cs,-val*uana(xg(i),yg(jj),zg(iii),tg(jjjj)));
                else
                  A.add(k,s(i,jj,iii,jjjj),Val);
                end if;
              end loop; -- jjjj
            end loop; -- jj

            -- for d^2u/dzdt
            coefdzt   := 10.0;
            rderiv(1,nz,iii,hz,cz);
            rderiv(1,nt,iiii,ht,ct);
            for jjj in 1..nz loop
              for jjjj in 1..nt loop
                val := coefdzt*cz(jjj)*ct(jjjj);
                if jjj=1 or jjj=nz or jjjj=1 or jjjj=nt then
                  A.add(k,cs,-val*uana(xg(i),yg(ii),zg(jjj),tg(jjjj)));
                else
                  A.add(k,s(i,ii,jjj,jjjj),Val);
                end if;
              end loop; -- jjjj
            end loop; -- jjj

            -- for du/dt
            coefdt    := 11.0;
            rderiv(1,nt,iiii,ht,ct);
            for jjjj in 1..nt loop
              val := coefdt*ct(jjjj);
              if jjjj=1 or jjjj=nt then
                A.add(k,cs,-val*uana(xg(i),yg(ii),zg(iii),tg(jjjj)));
              else
                A.add(k,s(i,ii,iii,jjjj),Val);
              end if;
            end loop; -- jjjj
            -- for u(x,y,z,t)
            coefu     := 12.0;
            A.add(k,s(i,ii,iii,iiii),Coefu);

            -- for f(x,y,z,t) RHS constant
            A.add(k,cs,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("pde44_sparse.adb running ");
    put_line("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + ");
    put_line("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + ");
    put_line("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) + ");
    put_line("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = ");
    put_line("  F(x,y,z,t) ");
    put_line("  F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +   ");
    put_line("               180 z^2*t + 200 y^2z*t + 44 t^3 +    ");
    put_line("               110 y^2z^2t + 66 xy + 12t^4 +        ");
    put_line("               24 z^4 + 36 y^4 + 48 x^4 +           ");
    put_line("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t ");
    put_line("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    put_line("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 + ");
    put_line("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 ");
    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 := calendar.seconds(calendar.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 := calendar.seconds(calendar.clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    -- initializeA;
    -- A.Init ??

    time_now := calendar.seconds(calendar.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
    -- A.Write_All; -- optional


    put_line("A computed stiffness matrix ");
    time_now := calendar.seconds(calendar.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);
    A.Simeq(U);


    put_line("equations solved ");
    time_now := calendar.seconds(calendar.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 := calendar.seconds(calendar.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 -- pde44_sparse
   main;
   put_line("end pde44_sparse.adb");
end pde44_sparse;