-- pdenu23_eq.adb  3d non uniform grid boundary value PDE
--                 known solution for testing method
--                 u(x,y) = x^3+2y^3+3z^3+4yx^2+5zy^2+6xz^2+7xy+8z+9
--
--                 differential equation to solve
-- d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dz^2 + 4*d^2u/dxdy + 5*d^2u/dxdz +
-- 6*d^2u/dydz + 7*du/dx + 8*du/dy + 9du/dz + 10*u=
-- 10x^3 + 20y^3 + 30z^3 + 40yx^2 + 50zy^2 + 60xz^2 + 53x^2 + 93y^2 + 123z^2 +
-- 126xy + 108zx + 80yz + 130x + 141y + 214z + 190
-- non uniform grid on cube (-1.0,-1.1,-1.2) to (1.2,1.3,1.4)
--
-- set up and solve system of linear equations
--
-- create two dimensional array with
-- (nx-1)*(ny-1)*(nz-1) rows, One for each equation


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;
with nuderiv;

procedure pdenu23_eq is
  nx : constant := 10; -- including boundary at 0 and nx, 0..nx
  ny : constant := 9;  -- including boundary at 0 and ny, 0..ny
  nz : constant := 8;  -- including boundary at 0 and nz, 0..nz
  xg : real_vector(0..nx) := (-1.0, -0.9, -0.8, -0.6, -0.3, 0.0,
                               0.35, 0.65, 0.85, 1.0,  1.2);
  yg : real_vector(0..ny) := (-1.1, -1.0, -0.9, -0.6, -0.2,
                               0.3,  0.8,  1.0,  1.2,  1.3);
  zg : real_vector(0..nz) := (-1.2, -1.0, -0.9, -0.5,  0.0,
                               0.5,  0.9,  1.2,  1.4);
  -- nx-1 by ny-1 by nz-1 internal, non boundary cells
  cxx : real_vector(0..nx); -- nuderiv Coefficients
  cyy : real_vector(0..ny);
  czz : real_vector(0..nz);
  cxy : real_matrix(0..nx,0..ny);
  cxz : real_matrix(0..nx,0..nz);
  cyz : real_matrix(0..ny,0..nz);
  cx : real_vector(0..nx);
  cy : real_vector(0..ny);
  cz : real_vector(0..nz);
  coefdxx, coefdyy, coefdzz, coefdxdy, coefdxdz, coefdydz : real;
  coefdx, coefdy, coefdz, coefu : real;
  us : real_vector(0..(nx-1)*(ny-1)*(nz-1)-1); -- solution vector
  ut : real_matrix(0..(nx-1)*(ny-1)*(nz-1)-1, 0..(nx-1)*(ny-1)*(nz-1));
                                        -- temporary for partial results
                       -- nx-1 by ny-1 by nz-1 internal, non boundary cells
  error : real;        -- sum of absolute exact-computed
  avg_error : real;    -- average error
  max_error : real;    -- maximum cell error

  package fltio is new float_io(real);
  use fltio;
  package intio is new integer_io(integer);
  use intio;

  function u(x:real; y:real; z:real) return real is
           -- solution for boundary and optionalcheck
  begin
    return x*x*x + 2.0*y*y*y + 3.0*z*z*z + 4.0*x*x*y + 5.0*z*y*y + 6.0*x*z*z +
           7.0*x*y + 8.0*z + 9.0;
  end u;

  function c(x:real; y:real; z:real) return real is
           -- RHS of differential equation to solve
  begin
   return 10.0*x*x*x + 20.0*y*y*y + 30.0*z*z*z + 40.0*x*x*y + 50.0*z*y*y +
          60.0*z*z*x + 53.0*x*x + 93.0*y*y + 123.0*z*z + 126.0*x*y +
          108.0*x*z + 80.0*y*z + 130.0*x + 141.0*y + 214.0*z + 190.0;
  end c;

  procedure printexact is
  begin
    new_line;
    put_line("exact solution");
    for i in 0..1 loop
      for j in 0..1 loop
        for k in 0..1 loop
          put("xg("&integer'image(1+i*(nx-2))&")=");
          put(xg(1+i*(nx-2)),3,2,0);
          put(", yg("&integer'image(1+j*(ny-2))&")=");
          put(yg(1+j*(ny-2)),3,2,0);
          put(", zg("&integer'image(1+k*(nz-2))&")=");
          put(zg(1+k*(nz-2)),3,2,0);
          put(", u(x,y,z)=");
          put(u(xg(1+i*(nx-2)),yg(1+j*(ny-2)),zg(1+k*(nz-2))),5,3,0);
          new_line;
        end loop;
      end loop;
    end loop;
  end printexact;

  procedure printus is
  begin
    new_line;
    put_line("computed solution");
    for i in 0..1 loop
      for j in 0..1 loop
        for k in 0..1 loop
          put("xg("&integer'image(1+i*(nx-2))&")=");
          put(xg(1+i*(nx-2)),3,2,0);
          put(", yg("&integer'image(1+j*(ny-2))&")=");
          put(yg(1+j*(ny-2)),3,2,0);
          put(", zg("&integer'image(1+k*(nz-2))&")=");
          put(zg(1+k*(nz-2)),3,2,0);
          put(", u(x,y,z)=");
          put(us((i*(nx-2))+(nx-1)*(j*(ny-2))+(nx-1)*(ny-1)*(k*(nz-2))),5,3,0);
          new_line;
        end loop;
      end loop;
    end loop;
  end printus;

  procedure printdiff is
  begin
    new_line;
    put_line("difference exact-computed solution");
    for i in 0..1 loop
      for j in 0..1 loop
        for k in 0..1 loop
          put("xg("&integer'image(1+i*(nx-2))&")=");
          put(xg(1+i*(nx-2)),3,2,0);
          put(", yg("&integer'image(1+j*(ny-2))&")=");
          put(yg(1+j*(ny-2)),3,2,0);
          put(", zg("&integer'image(1+k*(nz-2))&")=");
          put(zg(1+k*(nz-2)),3,2,0);
          put(", diff=");
          put(u(xg(1+i*(nx-2)),yg(1+j*(ny-2)),zg(1+k*(nz-2)))-
              us((i*(nx-2))+(nx-1)*(j*(ny-2))+(nx-1)*(ny-1)*(k*(nz-2))));
          new_line;
        end loop;
      end loop;
    end loop;
  end printdiff;

  procedure check is
    zexact, err : real;
  begin
    error := 0.0;
    max_error := 0.0;
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        for k in 1..nz-1 loop
          zexact := u(xg(i),yg(j),zg(k));
          err := abs(us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1))-zexact);
          error := error + err;
          if err>max_error then
            max_error := err;
          end if;
        end loop;
      end loop;
    end loop;
  end check;

  procedure init_matrix is
    ii, mjk, imk, ijm, kxky, kxkz, kykz, cs : integer;
  begin
    cs := (nx-1)*(ny-1)*(nz-1); -- constant column subscript
    -- zero internal cells
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        for k in 1..nz-1 loop
          for ii in 1..nx-1 loop
            for jj in 1..ny-1 loop
              for kk in 1..nz-1 loop
                 ut((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1),
                    (ii-1)+(nx-1)*(jj-1)+(nx-1)*(ny-1)*(kk-1)) := 0.0;
              end loop;
            end loop;
          end loop;
          ut((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1),cs) := 0.0;
        end loop;
      end loop;
    end loop;

    put_line("internal cells zeroed ");
    for i in 1..nx-1 loop -- inside boundary
      for j in 1..ny-1 loop
        for k in 1..nz-1 loop
          ii := (i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1);
          -- for each term
            -- coefdxx * d^2 u(x,y,z)/dx^2
            coefdxx := 1.0;
            nuderiv(2,nx+1,i,xg,cxx);
            for m in 0..nx loop
              cxx(m) := cxx(m) * coefdxx;
              if m=0 or m=nx then
                ut(ii,cs) := ut(ii,cs) - cxx(m)*u(xg(m),yg(j),zg(k));
              else
                mjk := (m-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1);
                ut(ii,mjk) := ut(ii,mjk) + cxx(m);
              end if;
            end loop;

            -- coefdyy * d^2 u(x,y,z)/dy^2
            coefdyy := 2.0;
            nuderiv(2,ny+1,j,yg,cyy);
            for m in 0..ny loop
              cyy(m) := cyy(m) * coefdyy;
              if m=0 or m=ny then
                ut(ii,cs) := ut(ii,cs) - cyy(m)*u(xg(i),yg(m),zg(k));
              else
                imk := (i-1)+(nx-1)*(m-1)+(nx-1)*(ny-1)*(k-1);
                ut(ii,imk) := ut(ii,imk) + cyy(m);
              end if;
            end loop;

            -- coefdzz * d^2 u(x,y,z)/dz^2
            coefdzz := 3.0;
            nuderiv(2,nz+1,k,zg,czz);
            for m in 0..nz loop
              czz(m) := czz(m) * coefdzz;
              if m=0 or m=nz then
                ut(ii,cs) := ut(ii,cs) - czz(m)*u(xg(i),yg(j),zg(m));
              else
                ijm := (i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(m-1);
                ut(ii,ijm) := ut(ii,ijm) + czz(m);
              end if;
            end loop;

            -- coefdxdy * d^2 u(x,y,z)/dxdy
            coefdxdy := 4.0;
            nuderiv(1,nx+1,i,xg,cx);
            nuderiv(1,ny+1,j,yg,cy);
            for kx in 0..nx loop
              for ky in 0..ny loop
                cxy(kx,ky) := coefdxdy * cx(kx) * cy(ky);
                if kx=0 or kx=nx or ky=0 or ky=ny then
                  ut(ii,cs) := ut(ii,cs) - cxy(kx,ky)*u(xg(kx),yg(ky),zg(k));
                else
                  kxky := (kx-1)+(nx-1)*(ky-1)+(nx-1)*(ny-1)*(k-1);
                  ut(ii,kxky) := ut(ii,kxky) + cxy(kx,ky);
                end if;
              end loop;
            end loop;

            -- coefdxdz * d^2 u(x,y,z)/dxdz
            coefdxdz := 5.0;
            nuderiv(1,nx+1,i,xg,cx);
            nuderiv(1,nz+1,k,zg,cz);
            for kx in 0..nx loop
              for kz in 0..nz loop
                cxz(kx,kz) := coefdxdz * cx(kx) * cz(kz);
                if kx=0 or kx=nx or kz=0 or kz=nz then
                  ut(ii,cs) := ut(ii,cs) - cxz(kx,kz)*u(xg(kx),yg(j),zg(kz));
                else
                  kxkz := (kx-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(kz-1);
                  ut(ii,kxkz) := ut(ii,kxkz) + cxz(kx,kz);
                end if;
              end loop;
            end loop;

            -- coefdydz * d^2 u(x,y,z)/dydz
            coefdydz := 6.0;
            nuderiv(1,ny+1,j,yg,cy);
            nuderiv(1,nz+1,k,zg,cz);
            for ky in 0..ny loop
              for kz in 0..nz loop
                cyz(ky,kz) := coefdydz * cy(ky) * cz(kz);
                if ky=0 or ky=ny or kz=0 or kz=nz then
                  ut(ii,cs) := ut(ii,cs) - cyz(ky,kz)*u(xg(i),yg(ky),zg(kz));
                else
                  kykz := (i-1)+(nx-1)*(ky-1)+(nx-1)*(ny-1)*(kz-1);
                  ut(ii,kykz) := ut(ii,kykz) + cyz(ky,kz);
                end if;
              end loop;
            end loop;

            -- coefdx * d u(x,y)/dx
            coefdx := 7.0;
            nuderiv(1,nx+1,i,xg,cx);
            for m in 0..nx loop
              cx(m) := cx(m) * coefdx;
              if m=0 or m=nx then
                ut(ii,cs) := ut(ii,cs) - cx(m)*u(xg(m),yg(j),zg(k));
              else
                mjk := (m-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1);
                ut(ii,mjk) := ut(ii,mjk) + cx(m);
              end if;
            end loop;

            -- coefdy * d u(x,y)/dy
            coefdy := 8.0;
            nuderiv(1,ny+1,j,yg,cy);
            for m in 0..ny loop
              cy(m) := cy(m) * coefdy;
              if m=0 or m=ny then
                ut(ii,cs) := ut(ii,cs) - cy(m)*u(xg(i),yg(m),zg(k));
              else
                imk := (i-1)+(nx-1)*(m-1)+(nx-1)*(ny-1)*(k-1);
                ut(ii,imk) := ut(ii,imk) + cy(m);
              end if;
            end loop;

            -- coefdy * d u(x,y)/dy
            coefdz := 9.0;
            nuderiv(1,nz+1,k,zg,cz);
            for m in 0..nz loop
              cz(m) := cz(m) * coefdz;
              if m=0 or m=nz then
                ut(ii,cs) := ut(ii,cs) - cz(m)*u(xg(i),yg(j),zg(m));
              else
                ijm := (i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(m-1);
                ut(ii,ijm) := ut(ii,ijm) + cz(m);
              end if;
            end loop;

            -- coefu * u(x,y)
            coefu := 10.0;
            ut(ii,ii) := ut(ii,ii) + 1.0*coefu;

            -- RHS constant
            ut(ii,cs) := ut(ii,cs) + c(xg(i),yg(j),zg(k));
          -- end terms
        end loop;
      end loop;
    end loop;
  end init_matrix;

  procedure simeq(n : integer; b : in out real_matrix; x : out real_vector) is
    --  tailored version for this application, zero based, in place
    row : array(0..n)of integer;   -- row interchange indices
    hold , i_pivot : integer;      -- pivot indices
    pivot : real;            -- pivot element value
    abs_pivot : real;        -- abs of pivot element
    norm1 : real := 0.0;     -- 1 norm of matrix
    matrix_data_error : exception;
  begin
    for k in 0..n loop
      row(k):= k;
    end loop;
    --                               begin main reduction loop
    for k in 0..n loop
      --                             find largest element for pivot
      pivot := b(row(k), k);
      abs_pivot := abs(pivot);
      i_pivot := k;
      for i in k..n loop
        if abs(b(row(i), k)) > abs_pivot then
          i_pivot := i;
          pivot := b(row(i), k);
          abs_pivot := abs(pivot);
        end if;
      end loop;
      --                             check for near singular
      if abs_pivot < real'epsilon * norm1 then
        raise matrix_data_error;
      end if;
      --                             have pivot, interchange row pointers
      hold := row(k);
      row(k):= row(i_pivot);
      row(i_pivot):= hold;
      --                           reduce about pivot
      for j in k+1..n+1 loop
        b(row(k), j) := b(row(k), j) / b(row(k), k);
      end loop;
      --                           inner reduction loop
      for i in 0..n loop
        if i /= k then
          for j in k+1..n+1 loop
            b(row(i), j) := b(row(i), j) - b(row(i), k) * b(row(k), j);
          end loop;
        end if;
      end loop;
      --                             finished inner reduction
    end loop;
    --                               build  x  for return, unscrambling rows
    for i in 0..n loop
      x(i):= b(row(i), n+1);
    end loop;
  end simeq;

begin  -- pdenu23_eq

  put_line("pdenu23_eq.adb running");
  put("xmin=");
  put(xg(0),3,3,0);
  put(", xmax=");
  put(xg(nx),3,3,0);
  put_Line(", nx="&integer'image(nx));
  put("ymin=");
  put(yg(0),3,3,0);
  put(", ymax=");
  put(yg(ny),3,3,0);
  put_line(", ny="&integer'image(ny));
  put("zmin=");
  put(zg(0),3,3,0);
  put(", zmax=");
  put(zg(nz),3,3,0);
  put_line(", nz="&integer'image(nz));
  put("u(xg(1),yg(1),zg(1))=");
  put(u(xg(1),yg(1),zg(1)),4,3,0); new_line;
  put("c(xg(1),yg(1),zg(1))=");
  put(c(xg(1),yg(1),zg(1)),4,3,0); new_line;

  init_matrix;
  put_line("matrix initialized");
  put_line("initial matrix, left side, upper");
  for i in 0..6 loop
    for j in 0..6 loop
      put(ut(i,j),6,3,0);
    end loop;
    new_line;
  end loop;
  new_line;

  put_line("initial matrix, right side, upper");
  for i in 0..6 loop
    for j in (nx-1)*(ny-1)*(nz-1)-6..(nx-1)*(ny-1)*(nz-1) loop
      put(ut(i,j),6,3,0);
    end loop;
    new_line;
  end loop;
  new_line;

  put_line("initial matrix, left side, lower");
  for i in (nx-1)*(ny-1)*(nz-1)-7..(nx-1)*(ny-1)*(nz-1)-1 loop
    for j in 0..6 loop
      put(ut(i,j),6,3,0);
    end loop;
    new_line;
  end loop;
  new_line;

  put_line("initial matrix, right side, lower");
  for i in (nx-1)*(ny-1)*(nz-1)-7..(nx-1)*(ny-1)*(nz-1)-1 loop
    for j in (nx-1)*(ny-1)*(nz-1)-6..(nx-1)*(ny-1)*(nz-1) loop
      put(ut(i,j),6,3,0);
    end loop;
    new_line;
  end loop;
  new_line;

  simeq((nx-1)*(ny-1)*(nz-1)-1, ut, us); -- tailored version, returned us

  printexact;
  printus;
  printdiff;
  check;
  avg_error := error/(real(nx-1)*real(ny-1)*real(nz-1));
  put_line("error="&real'image(error)&" avg="&real'image(avg_error)&
           ", max="&real'image(max_error));
end pdenu23_eq;