-- pdenu22_eq.adb  2d non uniform grid boundary value PDE
--                 known solution for testing method
--                 u(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8
--
--                 differential equation to solve
-- d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dxdy + 4*du/dx + 5*du/dy +6*u=
-- yx^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110
-- non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3
--
-- set up and solve system of linear equations
--
-- create two dimensional array with
-- (nx-1)*(ny-1)*(ny-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 pdenu22_eq is
  nx : constant := 10; -- including boundary at 0 and nx, 0..nx
  ny : constant := 9;  -- including boundary at 0 and ny, 0..ny
  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);
  -- nx-1 by ny-1  internal, non boundary cells
  cxx : real_vector(0..nx); -- nuderiv Coefficients
  cyy : real_vector(0..ny);
  cxy : real_matrix(0..nx,0..ny);
  cx : real_vector(0..nx);
  cy : real_vector(0..ny);
  coefdxx, coefdyy, coefdxdy, coefdx, coefdy, coefu : real;
  us : real_vector(0..(nx-1)*(ny-1)-1); -- solution vector
  ut : real_matrix(0..(nx-1)*(ny-1)-1, 0..(nx-1)*(ny-1));
                                        -- temporary for partial results
                            -- nx-1 by ny-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) return real is
           -- solution for boundary and optionalcheck
  begin
    return x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*x +
           6.0*y*y + 7.0*x + 8.0;
  end u;

  function c(x:real; y:real) return real is
           -- RHS of differential equation to solve
  begin
    return 6.0*x*x*x + 12.0*y*y*y + 18.0*x*x*y + 24.0*x*y*y + 64.0*x*y +
           57.0*x*x + 82.0*y*y + 122.0*x + 114.0*y + 110.0;
  end c;

  procedure printexact is
  begin
    new_line;
    put_line("exact solution u, computed us, error");
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
         put("xg("&integer'image(i)&")=");
         put(xg(i),3,3,0);
         put(", yg("&integer'image(j)&")=");
         put(xg(j),3,3,0);
         put(", u=");
         put(u(xg(i),yg(j)),3,3,0);
         put(", us=");
         put(us((i-1)+(nx-1)*(j-1)),3,3,0);
         put(", err=");
         put_line(real'image(us((i-1)+(nx-1)*(j-1))-u(xg(i),yg(j))));
      end loop;
    end loop;
  end printexact;

  procedure printus is
  begin
    new_line;
    put_line("computed solution");
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
         put("xg("&integer'image(i)&")=");
         put(xg(i),3,3,0);
         put(", yg("&integer'image(j)&")=");
         put(xg(j),3,3,0);
      end loop;
    end loop;
  end printus;

  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
        zexact := u(xg(i),yg(j));
        err := abs(us((i-1)+(nx-1)*(j-1))-zexact);
        error := error + err;
        if err>max_error then
          max_error := err;
        end if;
      end loop;
    end loop;
  end check;

  procedure init_matrix is
    ii, kj, ik, kxky, cs : integer;
  begin
    cs := (nx-1)*(ny-1); -- constant column subscript
    -- zero internal cells
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        for ii in 1..nx-1 loop
          for jj in 1..ny-1 loop
            ut((i-1)+(nx-1)*(j-1),(ii-1)+(nx-1)*(jj-1)) := 0.0;
          end loop;
        end loop;
        ut((i-1)+(nx-1)*(j-1),cs) := 0.0;
      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
        ii := (i-1)+(nx-1)*(j-1);
        -- put_line("ii="&integer'image(ii)&" at i="&integer'image(i)&
        --                                    ", j="&integer'image(j));
        -- for each term
          -- coefdxx * d^2 u(x,y)/dx^2
          coefdxx := 1.0;
          nuderiv(2,nx+1,i,xg,cxx);
          for k in 0..nx loop
            cxx(k) := cxx(k) * coefdxx;
            if k=0 or k=nx then
              ut(ii,cs) := ut(ii,cs) - cxx(k)*u(xg(k),yg(j));
            else
              kj := (k-1)+(nx-1)*(j-1);
              ut(ii,kj) := ut(ii,kj) + cxx(k);
            end if;
          end loop;
          -- coefdyy * d^2 u(x,y)/dy^2
          coefdyy := 2.0;
          nuderiv(2,ny+1,j,yg,cyy);
          for k in 0..ny loop
            cyy(k) := cyy(k) * coefdyy;
            if k=0 or k=ny then
              ut(ii,cs) := ut(ii,cs) - cyy(k)*u(xg(i),yg(k));
            else
              ik := (i-1)+(nx-1)*(k-1);
              ut(ii,ik) := ut(ii,ik) + cyy(k);
            end if;
          end loop;
          -- coefdxdy * d^2 u(x,y)/dxdy
          coefdxdy := 3.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));
              else
                kxky := (kx-1)+(nx-1)*(ky-1);
                ut(ii,kxky) := ut(ii,kxky) + cxy(kx,ky);
              end if;
            end loop;
          end loop;
          -- coefdx * d u(x,y)/dx
          coefdx := 4.0;
          nuderiv(1,nx+1,i,xg,cx);
          for k in 0..nx loop
            cx(k) := cx(k) * coefdx;
            if k=0 or k=nx then
              ut(ii,cs) := ut(ii,cs) - cx(k)*u(xg(k),yg(j));
            else
              kj := (k-1)+(nx-1)*(j-1);
              ut(ii,kj) := ut(ii,kj) + cx(k);
            end if;
          end loop;
          -- coefdy * d u(x,y)/dy
          coefdy := 5.0;
          nuderiv(1,ny+1,j,yg,cy);
          for k in 0..ny loop
            cy(k) := cy(k) * coefdy;
            if k=0 or k=ny then
              ut(ii,cs) := ut(ii,cs) - cy(k)*u(xg(i),yg(k));
            else
              ik := (i-1)+(nx-1)*(k-1);
              ut(ii,ik) := ut(ii,ik) + cy(k);
            end if;
          end loop;
          -- coefu * u(x,y)
          coefu := 6.0;
          ut(ii,ii) := ut(ii,ii) + 1.0*coefu;
          -- RHS constant
          ut(ii,cs) := ut(ii,cs) + c(xg(i),yg(j));
        -- end terms
      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  -- pdenu22_eq

  put_line("pdenu22_eq.adb running");
  put_line("differential equation to solve");
  put_line("d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dxdy + 4*du/dx + 5*du/dy +6*u=c(x,y)");
  put_line("c(x,y)=yx^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110");
  put_line("non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3");
  put_line("known Solution, U(x,y), for testing method");
  put_line("u(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8");
  put("xmin=");
  put(xg(0),4,3,0);
  put(", xmax=");
  put(xg(nx),4,3,0);
  put_line(", nx="&integer'image(nx));
  put("ymin=");
  put(yg(0),4,3,0);
  put(", ymax=");
  put(yg(ny),4,3,0);
  put_line(", ny="&integer'image(ny));
  put("u(xg(1),yg(1))=");
  put(u(xg(1),yg(1)),5,3,0); new_line;
  put("c(xg(1),yg(1))=");
  put(c(xg(1),yg(1)),5,3,0); new_line;

  init_matrix;
  put_line("matrix initialized");

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

  printexact;

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