-- pde47h_eq.adb  Biharmonic 4th order PDE in seven dimensions
-- sample problem:
-- solve Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) +
--       Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) +
--       Uwwww(x,y,z,t,u,v,w) +
--       2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) +
--       2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) +
--       2 Uww(x,y,z,t,u,v,w) +
--       8 U(x,y,z,t,u,v,w) = f(x,y,z,t,u,v,w)
--
-- boundary conditions computed using Ub(x,y,z,t,u,v,w)
-- analytic solution is U(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)

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 pde47h_eq is
  nx : integer := 5; -- problem parameters
  ny : integer := 5;
  nz : integer := 5;
  nt : integer := 5;
  nu : integer := 5;
  nv : integer := 5;
  nw : integer := 5;
  nxyztuvw : integer := (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2);

  A: real_matrix(1..nxyztuvw,1..nxyztuvw); -- A * U = R
  R : real_vector(1..nxyztuvw);  -- RHS known
  UU : real_vector(1..nxyztuvw);  -- solution being computed - not boundary
  Ua : real_vector(1..nxyztuvw); -- solution for checking - not boundary
  xg : real_vector(1..nx);
  yg : real_vector(1..ny);
  zg : real_vector(1..nz);
  tg : real_vector(1..nt);
  ug : real_vector(1..nu);
  vg : real_vector(1..nv);
  wg : real_vector(1..nw);
  xmin : real := 0.0; -- problem parameters
  xmax : real := 0.5;
  ymin : real := 0.0;
  ymax : real := 0.5;
  zmin : real := 0.0;
  zmax : real := 0.5;
  tmin : real := 0.0;
  tmax : real := 0.5;
  umin : real := 0.0;
  umax : real := 0.5;
  vmin : real := 0.0;
  vmax : real := 0.5;
  wmin : real := 0.0;
  wmax : real := 0.5;
  hx, hy, hz, ht, hu, hv, hw : real;

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

  -- internal functions

  function S(i:integer; ii:integer; iii:integer;
             iiii:integer; iiiii:integer; iiiiii:integer;
             iiiiiii:integer)
            return integer is -- for internal x,y,z,t
  begin
    return (i-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2) +
           (ii-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2) +
           (iii-2)*(Nt-2)*(nu-2)*(nv-2)*(nw-2) +
           (iiii-2)*(nu-2)*(nv-2)*(nw-2) +
           (iiiii-2)*(nv-2)*(nw-2) +
           (iiiiii-2)*(nw-2) +
           (iiiiiii-1);
  end S;

  -- problem RHS,
  function f(x:real; y:real; z:real; t:real; u:real;
           v:real; w:real) return real is
  begin
     return sin(x+y+z+t+u+v+w);
  end f;

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

  function eval_derivative(xord:integer; yord:integer; zord:integer;
                           tord:integer; uord:integer; vord:integer;
                           word:integer;
                           i:integer; ii:integer; iii:integer;
                           iiii:integer; iiiii:integer; iiiiii:integer;
                           iiiiiii:integer ;
                           UU: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);
    cu : real_vector(1..nu);
    cv : real_vector(1..nv);
    cw : real_vector(1..nw);
    discrete : real;
    x, y, z, t, u, v, w : real;
  begin
    x := xg(i);
    y := yg(ii);
    z := zg(iii);
    t := tg(iiii);
    u := ug(iiiii);
    v := vg(iiiiii);
    w := wg(iiiiiii);
    discrete := 0.0;

    -- seven cases of one partial x, y, z, t, u, v, w to any order
    if xord/=0 and yord=0 and zord=0 and tord=0 and uord=0 and
       vord=0 and word=0 then
      rderiv(xord, nx, i, hx, cx);
      for j in 1..nx loop
        if j=1 or j=nx then
          discrete := discrete + cx(j)*Ub(xg(j),y,z,t,u,v,w);
        else
           discrete := discrete + cx(j)*UU(s(j,ii,iii,iiii,iiiii,
                                            iiiiii,iiiiiii));
        end if;
      end loop; -- j x

    elsif xord=0 and yord/=0 and zord=0 and tord=0 and uord=0 and
          vord=0 and word=0 then
      rderiv(yord, ny, ii, hy, cy);
      for jj in 1..ny loop
        if jj=1 or jj=ny then
          discrete := discrete + cy(jj)*Ub(x,yg(jj),z,t,u,v,w);
        else
          discrete := discrete + cy(jj)*UU(s(i,jj,iii,iiii,iiiii,
                                            iiiiii,iiiiiii));
        end if;
      end loop; -- jj y

    elsif xord=0 and yord=0 and zord/=0 and tord=0 and uord=0 and
          vord=0 and word=0 then
      rderiv(zord, nz, iii, hz, cz);
      for jjj in 1..nz loop
        if jjj=1 or jjj=nz then
          discrete := discrete + cz(jjj)*Ub(x,y,zg(jjj),t,u,v,w);
        else
          discrete := discrete + cz(jjj)*UU(s(i,ii,jjj,iiii,iiiii,
                                             iiiiii,iiiiiii));
        end if;
      end loop; -- jjj z

    elsif xord=0 and yord=0 and zord=0 and tord/=0 and uord=0 and
          vord=0 and word=0 then
      rderiv(tord, nt, iiii, ht, ct);
      for jjjj in 1..nt loop
        if jjjj=1 or jjjj=nt then
          discrete := discrete + ct(jjjj)*Ub(x,y,z,tg(jjjj),u,v,w);
        else
          discrete := discrete + ct(jjjj)*UU(s(i,ii,iii,jjjj,iiiii,
                                              iiiiii,iiiiiii));
        end if;
      end loop; -- jjjj t

    elsif xord=0 and yord=0 and zord=0 and tord=0 and uord/=0 and
          vord=0 and word=0 then
      rderiv(uord, nu, iiiii, hu, cu);
      for jjjjj in 1..nu loop
        if jjjjj=1 or jjjjj=nu then
          discrete := discrete + cu(jjjjj)*Ub(x,y,z,t,ug(jjjjj),v,w);
        else
          discrete := discrete + cu(jjjjj)*UU(s(i,ii,iii,iiii,jjjjj,
                                               iiiiii,iiiiiii));
        end if;
      end loop; -- jjjjj u

    elsif xord=0 and yord=0 and zord=0 and tord=0 and uord=0 and
          vord/=0 and word=0 then
      rderiv(vord, nv, iiiiii, hv, cv);
      for jjjjjj in 1..nv loop
        if jjjjjj=1 or jjjjjj=nv then
          discrete := discrete + cv(jjjjjj)*Ub(x,y,z,t,u,vg(jjjjjj),w);
        else
          discrete := discrete + cv(jjjjjj)*UU(s(i,ii,iii,iiii,iiiii,
                                                jjjjjj,iiiiiii));
        end if;
      end loop; -- jjjjjj v

    elsif xord=0 and yord=0 and zord=0 and tord=0 and uord=0 and
          vord=0 and word/=0 then
      rderiv(word, nw, iiiiiii, hw, cw);
      for jjjjjjj in 1..nw loop
        if jjjjjjj=1 or jjjjjjj=nw then
          discrete := discrete + cw(jjjjjjj)*Ub(x,y,z,t,u,v,wg(jjjjjjj));
        else
          discrete := discrete + cw(jjjjjjj)*UU(s(i,ii,iii,iiii,iiiii,
                                                iiiiii,jjjjjjj));
        end if;
      end loop; -- jjjjjjj w
    end if;

    -- other cases not needed for this PDE
    return discrete;
  end eval_derivative;

  procedure check_soln(UUU:real_vector) is
    val, err, smaxerr : real;
    bi, bii, biii, biiii, biiiii, biiiiii, biiiiiii : integer;
  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
            for iiiii in 2..nu-1 loop
              for iiiiii in 2..nv-1 loop
                for iiiiiii in 2..nw-1 loop
    val := 0.0;
    -- + Uxxxx(x,y,z,t,u,v,w)
    val := val + eval_derivative(4,0,0,0,0,0,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Uyyyy(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,4,0,0,0,0,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Uzzzz(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,0,4,0,0,0,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Utttt(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,0,0,4,0,0,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Uuuuu(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,0,0,0,4,0,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Uvvvv(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,0,0,0,0,4,0,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + Uwwww(x,y,z,t,u,v,w)
    val := val + eval_derivative(0,0,0,0,0,0,4,
                                 i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uxx(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(2,0,0,0,0,0,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uyy(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,2,0,0,0,0,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uzz(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,0,2,0,0,0,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Utt(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,0,0,2,0,0,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uuu(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,0,0,0,2,0,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uvv(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,0,0,0,0,2,0,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 2 Uww(x,y,z,t,u,v,w)
    val := val + 2.0*eval_derivative(0,0,0,0,0,0,2,
                                     i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU);
    -- + 8 u(x,y,z,t,u,v,w)
    val := val + 8.0*Ub(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii),
                        vg(iiiiii),wg(iiiiiii));
    --  - RHS
    val := val - f(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii),
                   vg(iiiiii),wg(iiiiiii));
    if(i=nx-1 and ii=ny-1 and iii=nz-1 and iiii=nt-1 and iiiii=nu-1 and
         iiiiii=nv-1 and iiiiiii=nw-1) then
      put_line("UUU="&real'image(UUU(s(i,ii,iii,iiii,iiiii,iiiiii,Iiiiiii))));
      put_line("w''''="&real'image(eval_derivative(0,0,0,0,0,0,4,
               i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU)));
      put_line("2w''="&real'image(2.0*eval_derivative(0,0,0,0,0,0,2,
               i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,UUU)));
      put_line("8Ub="&real'image(8.0*Ub(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii),
                   vg(iiiiii),wg(iiiiiii))));
      put_line("f="&real'image(f(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii),
                   vg(iiiiii),wg(iiiiiii))));
    end if;

    -- should be zero
    err := abs(val);
    if err>smaxerr then
      smaxerr := err;
      bi := i;
      bii := ii;
      biii := iii;
      biiii := iiii;
      biiiii := iiiii;
      biiiiii := iiiiii;
      biiiiiii := iiiiiii;
    end if;
                end loop; --iiiiiii
              end loop; --iiiiii
            end loop; --iiiii
          end loop; --iiii
        end loop; --iii
      end loop; -- ii
    end loop; -- i
    put_line("check_soln maxerr="&real'image(smaxerr));
    put_line("at="&integer'image(bi)&integer'image(bii)&
               integer'image(biiii)&integer'image(biiii)&
               integer'image(biiiii)&integer'image(biiiiii)&
               integer'image(biiiiiii));
  end Check_soln;

  procedure write_soln(Which:Integer; tag:String) is
    fout : File_Type;
  begin
    create(fout, out_file, "pde47h_eq_ada"&tag&".dat");

    put_line("writing pde47h_eq_ada"&tag&".dat");
    put_line(fout,"7 "&integer'image(nx)&integer'image(ny)&
                       integer'image(nz)&integer'image(nt)&
                       integer'image(nu)&integer'image(nv)&
                       integer'image(nw));
    for  i in 1..nx loop
      for ii in 1..ny loop
        for  iii in 1..nz loop
          for iiii in 1..nt loop
            for  iiiii in 1..nu loop
              for iiiiii in 1..nv loop
                for iiiiiii in 1..nw 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," ");
                  put(fout,ug(iiiii),2,6,0);
                  put(fout," ");
                  put(fout,vg(iiiiii),2,6,0);
                  put(fout," ");
                  put(fout,wg(iiiiiii),2,6,0);
                  put(fout," ");
                  if Which=1 then -- expected solution
                    put(fout,Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                           ug(iiiii),vg(iiiiii),wg(iiiiiii)),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 or iiiii=1 or iiiii=nu or
                      iiiiii=1 or iiiiii=nv or
                      iiiiiii=1 or iiiiiii=nw then -- boundary
                      put(fout,Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                               ug(iiiii),vg(iiiiii),wg(iiiiiii)),4,6,0);
                    else
                      put(fout,UU(S(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)),4,6,0);
                    end if;
                  end if;
                  new_line(fout);
                end loop; -- iiiiiii
              end loop; -- iiiiii
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    close(fout);
    put_line("finished writing pde47h_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 iiiii in 2..nu-1 loop
              for iiiiii in 2..nv-1 loop
                for iiiiiii in 2..nw-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
                          for jjjjj in 2..nu-1 loop
                            for jjjjjj in 2..nv-1 loop
                              for jjjjjjj in 2..nw-1 loop
                                A(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii),
                                  s(j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj)):=0.0;
                              end loop; -- jjjjjjj
                            end loop; -- jjjjjj
                          end loop; -- jjjjj
                        end loop; -- jjjj
                      end loop; -- jjj
                    end loop; -- jj
                  end loop; -- j
                  R(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) := 0.0;
                end loop; -- iiiiiii
              end loop; -- iiiiii
            end loop; -- iiiii
          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);
    ctt, ctttt : real_vector(1..nt);
    cuu, cuuuu : real_vector(1..nu);
    cvv, cvvvv : real_vector(1..nv);
    cww, cwwww : real_vector(1..nw);
  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
            for iiiii in 2..nu-1 loop
              for iiiiii in 2..nv-1 loop
                for iiiiiii in 2..nw-1 loop
                  k := s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii); -- 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*Ub(xg(j),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(j,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(j,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- j x

                -- 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*Ub(xg(i),yg(jj),zg(iii),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,jj,iii,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,jj,iii,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jj y

                -- 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*Ub(xg(i),yg(ii),zg(jjj),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,ii,jjj,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,jjj,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjj z

                -- 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*Ub(xg(i),yg(ii),zg(iii),tg(jjjj),
                                          ug(iiiii),vg(iiiiii),Wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,jjjj,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,iii,jjjj,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjjj t

                -- for d^4U/du^4
                rderiv(4,nu,iiiii,hu,cuuuu);
                for jjjjj in 1..nu loop
                  val := cuuuu(jjjjj);
                  if jjjjj=1 or jjjjj=nu then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(jjjjj),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjjjj u

                -- for d^4U/dv^4
                rderiv(4,nv,iiiiii,hv,cvvvv);
                for jjjjjj in 1..nv loop
                  val := cvvvv(jjjjjj);
                  if jjjjjj=1 or jjjjjj=nv then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(jjjjjj),Wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,iiii,iiiii,jjjjjj,iiiiiii)) :=
                    A(k,s(i,ii,iii,iiii,iiiii,jjjjjj,iiiiiii)) + val;
                  end if;
                end loop; -- jjjjjj v

                -- for d^4U/dw^4
                rderiv(4,nw,iiiiiii,hw,cwwww);
                for jjjjjjj in 1..nw loop
                  val := cwwww(jjjjjjj);
                  if jjjjjjj=1 or jjjjjjj=nw then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(jjjjjjj));
                  else
                    A(k,s(i,ii,iii,iiii,iiiii,iiiiii,jjjjjjj)) :=
                    A(k,s(i,ii,iii,iiii,iiiii,iiiiii,jjjjjjj)) + val;
                  end if;
                end loop; -- jjjjjjj w

                -- 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*Ub(xg(j),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(j,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(j,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- j x

                -- 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*Ub(xg(i),yg(jj),zg(iii),tg(iiii),
                                        ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,jj,iii,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,jj,iii,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jj y

                -- 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*Ub(xg(i),yg(ii),zg(jjj),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,ii,jjj,iiii,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,jjj,iiii,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjj z

                -- for d^2U/dt^2
                rderiv(2,nt,iiii,ht,ctt);
                for jjjj in 1..nt loop
                  val := 2.0*ctt(jjjj);
                  if jjjj=1 or jjjj=nt then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(jjjj),
                                          ug(iiiii),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,jjjj,iiiii,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,iii,jjjj,iiiii,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjjj t

                -- for d^2U/du^2
                rderiv(2,nu,iiiii,hu,cuu);
                for jjjjj in 1..nu loop
                  val := 2.0*cuu(jjjjj);
                  if jjjjj=1 or jjjjj=nu then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(jjjjj),vg(iiiiii),wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)) :=
                    A(k,s(i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)) + val;
                  end if;
                end loop; -- jjjjj u

                -- for d^2U/dv^2
                rderiv(2,nv,iiiiii,hv,cvv);
                for jjjjjj in 1..nv loop
                  val := 2.0*cvv(jjjjjj);
                  if jjjjjj=1 or jjjjjj=nv then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(jjjjjj),wg(iiiiiii));
                  else
                    A(k,s(i,ii,iii,iiii,iiiii,jjjjjj,iiiiiii)) :=
                    A(k,s(i,ii,iii,iiii,iiiii,jjjjjj,iiiiiii)) + val;
                  end if;
                end loop; -- jjjjjj v

                -- for d^2U/dw^2
                rderiv(2,nw,iiiiiii,hw,cww);
                for jjjjjjj in 1..nw loop
                  val := 2.0*cww(jjjjjjj);
                  if jjjjjjj=1 or jjjjjjj=nw then
                    R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(iii),tg(iiii),
                                          ug(iiiii),vg(iiiiii),wg(jjjjjjj));
                  else
                    A(k,s(i,ii,iii,iiii,iiiii,iiiiii,jjjjjjj)) :=
                    A(k,s(i,ii,iii,iiii,iiiii,iiiiii,jjjjjjj)) + val;
                  end if;
                end loop; -- jjjjjjj w

                -- for U(x,y,z,t,u,v,w)
                A(k,s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) :=
                A(k,s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) + 8.0;

                -- for f(x,y,z,t,u,v,w) RHS constant
                R(k) := R(k) + f(xg(i),yg(ii),zg(iii),tg(iiii),
                                 ug(iiiii),vg(iiiiii),wg(iiiiiii));

                -- end terms

                end loop; -- iiiiiii
              end loop; -- iiiiii
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
  end BuildA;

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

  begin -- main
    put_line("pde47h_eq.adb running ");
    put_line("Given Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) + ");
    put_line("      Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) + ");
    put_line("      Uwwww(x,y,z,t,u,v,w) + ");
    put_line("      2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) + ");
    put_line("      2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) + ");
    put_line("      2 Uww(x,y,z,t,u,v,w) + ");
    put_line("      8   U(x,y,z,t,u,v,w) = f(x,y,z,t,u,v,w) ");

    put_line("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    put_line("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax Boundaries ");
    put_line("wmin<w<=wmax Boundaries ");
    put_line("Analytic solution U(x,y,z,t,u,v,w) =  sin(x+y+z+t+u+v+w)");
    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("umin="&real'image(umin)&", umax="&real'image(umax));
    put_line("vmin="&real'image(vmin)&", vmax="&real'image(vmax));
    put_line("wmin="&real'image(wmin)&", wmax="&real'image(wmax));
    put_line("nx="&integer'image(nx)&", ny="&integer'image(ny)&
             ", nz="&integer'image(nz)&", nt="&integer'image(nt)&
             ", nu="&integer'image(nu)&", nv="&integer'image(nv)&
             ", nw="&integer'image(nw));
    put_line(" ");
    time_start := seconds(clock);
    time_last := time_start;

    hx := (xmax-xmin)/real(nx-1);
    put_line("x divisor, hx="&real'image(hx));
    xg(1) := xmin;
    put_line("xg( 1)="&real'image(xg(1)));
    for i in 2..nx loop -- decreasing step size
      xg(i) := xmin+real(i-1)*hx;
      put_line("xg("&integer'image(i)&")="&real'image(xg(i)));
    end loop; -- i x
    put_line(" ");

    hy := (ymax-ymin)/real(ny-1);
    put_line("y divisor, hy="&real'image(hy));
    yg(1) := ymin;
    put_line("yg( 1)="&real'image(yg(1)));
    for ii in 2..ny loop -- decreasing step size
      yg(ii) := ymin+real(ii-1)*hy;
      put_line("yg("&integer'image(ii)&")="&real'image(yg(ii)));
    end loop; -- ii y
    put_line(" ");

    hz := (zmax-zmin)/real(nz-1);
    put_line("z divisor, hz="&real'image(hz));
    zg(1) := zmin;
    put_line("zg( 1)="&real'image(zg(1)));
    for iii in 2..nz loop -- decreasing step size
      zg(iii) := zmin+real(iii-1)*hz;
      put_line("zg("&integer'image(iii)&")="&real'image(zg(iii)));
    end loop; -- iii z
    put_line(" ");

    ht := (tmax-tmin)/real(nt-1);
    put_line("t divisor, hz="&real'image(ht));
    tg(1) := tmin;
    put_line("tg( 1)="&real'image(tg(1)));
    for iiii in 2..nt loop -- decreasing step size
      tg(iiii) := tmin+real(iiii-1)*ht;
      put_line("tg("&integer'image(iiii)&")="&real'image(tg(iiii)));
    end loop; -- iiii t
    put_line(" ");

    hu := (umax-umin)/real(nu-1);
    put_line("u divisor, hu="&real'image(hu));
    ug(1) := umin;
    put_line("ug( 1)="&real'image(ug(1)));
    for iiiii in 2..nu loop -- decreasing step size
      ug(iiiii) := umin+real(iiiii-1)*hu;
      put_line("ug("&integer'image(iiiii)&")="&real'image(ug(iiiii)));
    end loop; -- iiiii u
    put_line(" ");

    hv := (vmax-vmin)/real(nv-1);
    put_line("v divisor, hv="&real'image(hv));
    vg(1) := vmin;
    put_line("vg( 1)="&real'image(vg(1)));
    for iiiiii in 2..nv loop -- decreasing step size
      vg(iiiiii) := vmin+real(iiiiii-1)*hv;
      put_line("vg("&integer'image(iiiiii)&")="&real'image(vg(iiiiii)));
    end loop; -- iiiiii v
    put_line(" ");

    hw := (wmax-wmin)/real(nw-1);
    put_line("w divisor, hw="&real'image(hw));
    wg(1) := wmin;
    put_line("wg( 1)="&real'image(wg(1)));
    for iiiiiii in 2..nw loop -- decreasing step size
      wg(iiiiiii) := wmin+real(iiiiiii-1)*hw;
      put_line("wg("&integer'image(iiiiiii)&")="&real'image(wg(iiiiiii)));
    end loop; -- iiiiiii w
    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);
            for iiiii in 2..nu-1 loop
              u := ug(iiiii);
              for iiiiii in 2..nv-1 loop
                v := vg(iiiiii);
                for iiiiiii in 2..nw-1 loop
                  w := wg(iiiiiii);
   Ua(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)) := Ub(x,y,z,t,u,v,w);
                end loop; -- iiiiiii
              end loop; -- iiiiii
            end loop; -- iiiii
          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

    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(nxyztuvw)&" equations in "&
             integer'image(nxyztuvw)&" unknowns ");

    UU := 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(UU);
    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
            for iiiii in 2..nu-1 loop
              for iiiiii in 2..nv-1 loop
                for iiiiiii in 2..nw-1 loop
                err := abs(UU(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii))-
                          Ua(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)));
                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)&
                      ","&integer'image(iiiii)&","&integer'image(iiiiii)&
                      ","&integer'image(iiiiiii)&")=");
                put(UU(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)),3,6);
                put(", Ua=");
                put(Ua(s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)),3,6);
                put(", err=");
                put(err);
                new_line;
                end loop; -- iiiiiii
              end loop; -- iiiiii
            end loop; -- iiiii
          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(nxyztuvw)));
    put_line(" ");
  end Main;
begin -- pde47h_eq
   main;
   -- write_soln(1,"U"); -- analytic solution
   -- write_soln(2,""); -- computed solution for plot7d
   put_line("end pde47h_eq.adb");
end pde47h_eq;