-- pde_bihar44tl_eq.adb  using Ada tasks for parallel simeq
--
-- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + utttt(x,y,z,t) +
--       2 uxxyy(x,y,z,t) + 2 uxxzz(x,y,z,t) + 2 uxxtt(x,y,z,t) +
--       2 uyyzz(x,y,z,t) + 2 uyytt(x,y,z,t) + 2 uzztt(x,y,z,t) -
--      16 u(x,y,z,t) = F(x,y,z,t)
-- F(x,y,z,t) = 0.0
--
-- boundary conditions computed using ub(x,y,z,t) = sin(x,y,z,t)
-- analytic solution is u(x,y,z,t) = sin(x+y+z+t)
--

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 psimeq;
with rderiv;
with calendar;

procedure pde_bihar44tl_eq is
  nx : integer := 12; -- problem parameters
  ny : integer := 12;
  nz : integer := 12;
  nt : integer := 12;
  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
  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.57;
  ymin : real := 0.0;
  ymax : real := 1.57;
  zmin : real := 0.0;
  zmax : real := 1.57;
  tmin : real := 0.0;
  tmax : real := 1.57;
  hx, hy, hz, ht : real;
  cxx, cxxxx : real_matrix(1..nx,1..nx);
  cyy, cyyyy : real_matrix(1..ny,1..ny);
  czz, czzzz : real_matrix(1..nz,1..nz);
  ctt, ctttt : real_matrix(1..nt,1..nt);

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

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

  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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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)*ub(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 uxxyy(x,y,z,t) + 2 uxxzz(x,y,z,t) + 2 uxxtt(x,y,z,t) +
   --       2 uyyzz(x,y,z,t) + 2 uyytt(x,y,z,t) + 2 uzztt(x,y,z,t) -
   --      16 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);
            -- + 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 uxxyy(x,y,z,t)
            val := val + 2.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, u);
            -- + 2 uxxzz(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, u);
            -- + 2 uxxtt(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, u);
            -- + 2 uyyzz(x,y,z,t)
            val := val + 2.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, u);
            -- + 2 uyytt(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, u);
            -- + 2 uzztt(x,y,z,t)
            val := val + 2.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, u);
            -- - 16 u(x,y,z,t)
            val := val - 16.0*ub(xg(i),yg(ii),zg(iii),tg(iiii));
            -- - F(x,y,z,t)  zero this case
            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 against PDE maxerr="&real'image(smaxerr));
  end Check_Soln;

  procedure Write_Soln is
    Fout : File_Type;
  begin
    Create(Fout, Out_File, "pde_bihar44tl_eq_ada.dat");

    Put_Line("writing pde_bihar44tl_eq_ada.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 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,ub(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;
            New_Line(Fout);
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    Close(Fout);
    Put_Line("finished writing pde_bihar44tl_eq_ada.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 dsetup is
    ccx : real_vector(1..nx);
    ccy : real_vector(1..ny);
    ccz : real_vector(1..nz);
    cct : real_vector(1..nt);
  begin
    put_line("check rderiv(4,nx,1,hx,cxxxx");
    for I in 1..nx loop
      rderiv(4,nx,i,hx,ccx);
      for j in 1..nx loop
        cxxxx(i,j) := ccx(j);
      end loop;
    end loop;
    for i in 1..nx loop
      for j in 1..nx loop
         put_line("cxxxx("&integer'image(i)&","&integer'image(j)&
                    ")="&real'image(cxxxx(i,j)));
      end loop;
    end loop;
    put_line(" ");
    for i in 1..nx loop
      rderiv(2,nx,i,hx,ccx);
      for j in 1..nx loop
        cxx(i,j) := ccx(j);
      end loop;
    end loop;

    for ii in 1..ny loop
      rderiv(4,ny,ii,hy,ccy);
      for jj in 1..ny loop
        cyyyy(ii,jj) := ccy(jj);
      end loop;
    end loop;
    for ii in 1..ny loop
      rderiv(2,ny,ii,hy,ccy);
      for jj in 1..ny loop
        cyy(ii,jj) := ccy(jj);
      end loop;
    end loop;

    for iii in 1..nz loop
      rderiv(4,nz,iii,hz,ccz);
      for jjj in 1..nz loop
        czzzz(iii,jjj) := ccz(jjj);
      end loop;
    end loop;
    for iii in 1..nz loop
      rderiv(2,nz,iii,hz,ccz);
      for jjj in 1..nz loop
        czz(iii,jjj) := ccz(jjj);
      end loop;
    end loop;

    for iiii in 1..nt loop
      rderiv(4,nt,iiii,ht,cct);
      for jjjj in 1..nt loop
        ctttt(iiii,jjjj) := cct(jjjj);
      end loop;
    end loop;
    for iiii in 1..nt loop
      rderiv(2,nt,iiii,ht,cct);
      for jjjj in 1..nt loop
        ctt(iiii,jjjj) := cct(jjjj);
      end loop;
    end loop;

  end dsetup;

  procedure buildA is -- system of equations
    val : real := 0.0;
    k: integer;
  begin
    -- compute entries in A matrix, inside boundary
    put_line("compute A matrix and Y 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
            for j in 1..nx loop
              val := cxxxx(i,j);
              if j=1 or j=nx then
                R(k) := R(k) - val*ub(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
            for jj in 1..ny loop
              val := cyyyy(ii,jj);
              if jj=1 or jj=ny then
                R(k) := R(k) - val*ub(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
            for jjj in 1..nz loop
              val := czzzz(iii,jjj);
              if jjj=1 or jjj=nz then
                R(k) := R(k) - val*ub(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
            for jjjj in 1..nt loop
              val := ctttt(iiii,jjjj);
              if jjjj=1 or jjjj=nt then
                R(k) := R(k) - val*ub(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 2 d^4u/dx^2dy^2
            for j in 1..nx loop
              for jj in 1..ny loop
                  val := 2.0*cxx(i,j)*cyy(ii,jj);
                  if j=1 or j=nx or jj=1 or jj=ny then
                    R(k) := R(k) - val*ub(xg(j),yg(jj),zg(iii),tg(iiii));
                  else
                    A(k,s(j,jj,iii,iiii)) := A(k,s(j,jj,iii,iiii)) + val;
                  end if;
              end loop; -- jj
            end loop; -- j

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

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

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

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

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

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

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

            -- end terms
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    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("pde_bihar44tl_eq.adb running ");
    put_line("solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +");
    put_line("      utttt(x,y,z,t) + 2 uxxyy(x,y,z,t) + 2 uxxzz(x,y,z,t) +");
    put_line("    2 uxxtt(x,y,z,t) + 2 uyyzz(x,y,z,t) + 2 uyytt(x,y,z,t) +");
    put_line("    2 uzztt(x,y,z,t) - 16 u(x,y,z,t) = F(x,y,z,t)");
    put_line(" F(x,y,z,t) = 0.0 ");
    put_line(" ");
    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 := 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(" ");

    dsetup; -- precompute discretization coefficients

    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;

    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

    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 := psimeq(Nx-2, A, R);

    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 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
      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);
            err := abs(U(s(i,ii,iii,iiii))-ub(x,y,z,t));
            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(", ub=");
            put(ub(x,y,z,t),3,6); put(", err=");
            put(err);
            new_line;
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i

    write_soln;

    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 -- pde_bihar44tl_eq
   main;
   put_line("end pde_bihar44tl_eq.adb");
end pde_bihar44tl_eq;