-- pde2sin_eq.adb  2d uniform grid boundary value PDE
--                 known solution for testing method
--                 u(x,y) = sin(x-y)
--
--                 differential equation to solve
--                 du/dx + du/dy=0
-- uniform grid on rectangle 0 to 2Pi, 0 to 2Pi
-- set up and solve system of linear equations

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 Rderiv;

procedure pde2sin_eq is
  -- do NOT make nx=ny, this creates a singular matrix
  nx : constant := 13; -- including boundary at 1 and nx, 1..nx
  ny : constant := 12; -- including boundary at 1 and ny, 1..ny
  neqn : constant := (nx-2)*(ny-2);
  cs : integer := neqn+1;
  pi : real := 3.1415926535897932384626433832795028841971;
  xmin : real := 0.0;
  xmax : real := 2.0*pi;
  ymin : real := 0.0;
  ymax : real := 2.0*pi;
  hx : real := (xmax-xmin)/real(nx-1);
  hy : real := (ymax-ymin)/real(ny-1);
  xg : real_vector(1..nx);
  yg : real_vector(1..ny);
  -- nx-2 by ny-2  internal, boundary cells
  cx : real_vector(1..nx);
  cy : real_vector(1..ny);
  coefdx, coefdy: real;
  us : real_vector(1..neqn); -- solution vector
  ut : real_matrix(1..neqn, 1..neqn+1); -- includes RHS

  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 sin(x-y);
  end u;

  function c(x:real; y:real) return real is
           -- RHS of differential equation to solve
  begin
    return 0.0;
  end c;

  function S(i:integer; j:integer) return integer is
  begin
    return (i-2)*(ny-2)+(j-2)+1; -- index into us and ut
  end S;

  procedure print is
    error : real := 0.0;
    max_error : real := 0.0;
    avg_error : real := 0.0;
    val : real;
  begin
    new_line;
    put_line("exact solution u, computed us, error");
    for i in 1..nx loop
      for j in 1..ny loop
         put("xg("&integer'image(i)&")=");
         put(xg(i),3,3,0);
         put(", yg("&integer'image(j)&")=");
         put(yg(j),3,3,0);
         put(", u=");
         put(u(xg(i),yg(j)),3,3,0);
         if i=1 or i=nx or j=1 or j=ny then -- boundary
           val := u(xg(i),yg(j));
           error := 0.0;
         else
           val := us(S(i,j));
           error := abs(val-u(xg(i),yg(j)));
           avg_error := avg_error + error;
           if error > max_error then
              max_error := error;
           end if;
         end if;
         put(", us=");
         put(val,3,3,0);
         put(", err=");
         put_line(real'image(error));
      end loop;
    end loop;
    put_line("avg_error="&real'image(avg_error/real(neqn))&
             ", max_error="&real'image(max_error));
  end print;

  procedure check_row(row : Integer) is -- only for checking
    -- row is in ut solution coordinates
    ii : Integer;
    sum : Real;
  begin
    sum := 0.0;
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        ii := S(i,j);
        sum := sum + ut(row,ii)*u(xg(i),yg(j)); -- u is check solution
      end loop;
    end loop;
    sum := sum - ut(row,cs);
    if abs(sum)>0.001 then
       Put_Line("row="&Integer'Image(row)&" is bad err="&
                Real'Image(sum));
    end if;
  end check_row;

  procedure init_matrix is
    ii, kj, ik : integer;
  begin
    -- cs is RHS, constant column subscript
    -- zero internal cells
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        for ii in 2..nx-1 loop
          for jj in 2..ny-1 loop
            ut(S(i,j),S(ii,jj)) := 0.0;
          end loop;
        end loop;
        ut(S(i,j),cs) := 0.0;
      end loop;
    end loop;

    put_line("internal cells zeroed ");
    for i in 2..nx-1 loop -- inside boundary
      for j in 2..ny-1 loop
        ii := S(i,j); -- equation
          -- for each term

          -- coefdx * d u(x,y)/dx
          coefdx := 1.0;
          rderiv(1,nx,i,hx,cx);
          for k in 1..nx loop
            cx(k) := cx(k) * coefdx;
            if k=1 or k=nx then
              ut(ii,cs) := ut(ii,cs) - cx(k)*u(xg(k),yg(j));
            else
              kj := S(k,j);
              ut(ii,kj) := ut(ii,kj) + cx(k);
            end if;
          end loop; -- k
          -- coefdy * d u(x,y)/dy
          coefdy := 1.0;
          rderiv(1,ny,j,hy,cy);
          for k in 1..ny loop
            cy(k) := cy(k) * coefdy;
            if k=1 or k=ny then
              ut(ii,cs) := ut(ii,cs) - cy(k)*u(xg(i),yg(k));
            else
              ik := S(i,k);
              ut(ii,ik) := ut(ii,ik) + cy(k);
            end if;
          end loop; --k
          -- RHS constant
          ut(ii,cs) := ut(ii,cs) + c(xg(i),yg(j));
        -- end terms
        check_row(Ii);
      end loop; -- j
    end loop; -- i
  end init_matrix;

  procedure Print_B is
  begin
    Put_Line("matrix B includes RHS");
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        for ii in 2..nx-1 loop
          for jj in 2..ny-1 loop
             Put_Line("i="&Integer'Image(i)&", j="&Integer'Image(j)&
                      ",  ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)&
                       ", B(i,j)="&
                       Real'Image(ut(S(i,j),S(ii,jj))));
          end loop; -- jj
        end loop; -- ii
        Put_Line("i="&Integer'Image(i)&", j="&Integer'Image(j)&
                 ",                RHS="&
                       Real'Image(ut(S(i,j),cs)));
      end loop; -- j
    end loop; -- i
    New_Line;
    New_Line;
  end Print_B;

  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(1..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 := 1.0;     -- 1 norm of matrix
    matrix_data_error : exception;
  begin
    for k in 1..n loop
      row(k):= k;
    end loop;
    --                               begin main reduction loop
    for k in 1..n loop
      --                             find largest element for pivot
      pivot := b(row(k), k);
      abs_pivot := abs(pivot);
      i_pivot := k;
      for i in k+1..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 1..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 1..n loop
      x(i):= b(row(i), n+1);
    end loop;
  end simeq;

  procedure check_soln is -- check when solution unknown
    smaxerr, err, x, y, ux, uy : Real;
  begin
    --  ux(x,y) + uy(x,y) = c(x,y)
    Put_line("check_soln against PDE");
    smaxerr := 0.0;

    for i in 2..nx-1 loop -- through dof equations
      rderiv(1,nx,i,hx,cx);
      x := xg(I);
      for ii in 2..ny-1 loop
        rderiv(1,ny,ii,hy,cy);
        y := yg(ii);
        err := -c(x,y); -- add other side of equation terms
        ux := 0.0;
        for j in 1..nx loop -- compute numerical derivative at (i,ii)
          if j=1 or j=nx then
            ux := ux + cx(j)*U(xg(j),yg(ii)); -- boundary
          else
            ux := ux + cx(j)*Us(S(j,ii));
          end if;
        end loop;
        uy := 0.0;
        for jj in 1..ny loop
          if jj=1 or jj=ny then
            uy := uy + cy(jj)*U(xg(i),yg(jj)); -- boundary
          else
            uy := uy + cy(jj)*Us(S(i,jj));
          end if;
        end loop;
        err := err + ux + uy;
        if abs(err) > smaxerr then
          smaxerr := abs(err);
        end if;
      end loop; -- ii Y
    end loop; -- i X

    Put_line("check_soln max error="&Real'Image(smaxerr));
  end Check_Soln;

begin  -- pde2sin_eq

  put_line("pde2sin_eq.adb running");
  put_line("differential equation to solve");
  put_line("du/dx + du/dy = c(x,y)");
  put_line("c(x,y)=0");
  put_line("uniform grid on rectangle 0,2Pi to 0,2Pi");
  put_line("known Solution, for testing method");
  put_line("u(x,y) = sin(x-y)");
  put_line("do not make nx=ny, this creates a singular matrix");

  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;
  for j in 1..ny loop
    yg(j) := ymin+real(j-1)*hy;
    put_line("yg("&integer'image(j)&")="&real'image(yg(j)));
  end loop;
  put("xmin=");
  put(xmin,4,3,0);
  put(", xmax=");
  put(xmax,4,3,0);
  put(", hx=");
  put(hx,4,3,0);
  put_line(", nx="&integer'image(nx));
  put("ymin=");
  put(ymin,4,3,0);
  put(", ymax=");
  put(ymax,4,3,0);
  put(", hy=");
  put(hy,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;

  -- debug print
  -- Put_Line("debug print");
  -- Rderiv(1,Nx,3,Hx,Cx);
  -- for I in 1..Nx loop
  --   Put_Line("rderiv 1,3,"&Integer'Image(Nx)&" cx("&Integer'Image(I)&
  --            ")="&Real'Image(Cx(I)));
  -- end loop;
  -- Rderiv(1,Ny,3,Hy,Cy);
  -- for I in 1..Ny loop
  --    Put_Line("rderiv 1,3,"&Integer'Image(Ny)&" cy("&Integer'Image(I)&
  --             ")="&Real'Image(Cy(I)));
  -- end loop;

  init_matrix;
  put_line("matrix initialized");
  -- Print_B;
  simeq(neqn, ut, us); -- tailored version, returned us
  check_soln;
  print;
end pde2sin_eq;