-- pdenu_eq.adb  2d non uniform grid boundary value PDE  non iterative
--               known solution for testing method
--               u(x,y) = exp(x*y)*sin(x+y)
--
--               differential equation to solve
-- d^2u/dx^2 + 2*d^2u/dy^2 - 3*d^2u/dxdy + y*du/dx - x*du/dy =
-- exp(x*y)*sin(x+y)*(x*x + 2*y*y -3*x*y -3)
-- 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) 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 pdenu_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 exp(x*y)*sin(x+y);
  end u;

  function c(x:real; y:real) return real is
           -- RHS of differential equation to solve
  begin
    return exp(x*y)*sin(x+y)*(x*x + 2.0*y*y -3.0*x*y -3.0);
  end c;

  procedure printexact is
  begin
    new_line;
    put_line("exact solution");
    for i in 0..nx loop
      for j in 0..ny loop
        put_line("u("&integer'image(i)&","&integer'image(j)&
                 ")="&real'image(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_line("us("&integer'image(i)&","&integer'image(j)&
                 ")="&real'image(us((i-1)+(nx-1)*(j-1))));
      end loop;
    end loop;
  end printus;

  procedure printdiff is
  begin
    new_line;
    put_line("difference exact-computed solution");
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        put_line("diff("&integer'image(i)&","&integer'image(j)&
                 ")="&real'image(u(xg(i),yg(j))-us((i-1)+(nx-1)*(j-1))));
      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
        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);
              -- put_line("kj="&integer'image(kj)&" at k="&integer'image(k));
              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);
              -- put_line("ik="&integer'image(ik)&" at k="&integer'image(k));
              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 := yg(j);
          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 := -xg(i);
          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 := 0.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
    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  -- pdenu_eq

  put_line("pdenu_eq.adb running");
  put_line("xmin="&real'image(xg(0))&", xmax="&real'image(xg(nx))&
           ", nx="&integer'image(nx));
  put_line("ymin="&real'image(yg(0))&", ymax="&real'image(yg(ny))&
           ", ny="&integer'image(ny));
  init_matrix;
  put_line("matrix initialized");
  put_line("initial matrix, left side, upper");
  for i in 0..8 loop
    for j in 0..8 loop
      put(ut(i,j),4,3,0);
    end loop;
    new_line;
  end loop;
  new_line;

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

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

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

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

  printexact;
  printus;
  printdiff;
  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 pdenu_eq;