-- sparse_abc.adb  see WEB page for development CMSC 455 supplemental
--                 The PDE to be solved is:
--        a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
--        d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y) + f1(x,y)*u(x,y) = c(x,y)
--
-- The functions and other user information is in the files:
--                      abc.ads   abc.adb
--   also compile with  sparse.ads  sparse.adb   deriv.adb
--                      real_arrays.ads   real_arrays.adb

with Ada.Text_IO;
use  Ada.Text_IO;
with Real_Arrays;
use  Real_Arrays;
with abc;
use  abc;
with sparse;
with deriv;

procedure sparse_abc is
  -- nx and ny  are defined in abc.ads
  us:Real_Vector(1..(nx-2)*(ny-2));              -- solution being computed, X
  uu:Real_Vector(1..(nx-2)*(ny-2));              -- solution if checking
  -- S ut:Real_Matrix(1..(nx-2)*(ny-2),1..(nx-2)*(ny-2)); --solution matrix, A
  -- S uc:Real_Vector(1..(nx-2)*(ny-2));               -- equation constant, Y
  --   ut uc now inside A, zero based subscript uc is nrow
                      -- nx-2 by ny-2  internal, non boundary cells
  -- all elements 1..nx by 1..ny include boundary elements
  -- xmin, xmax, ymin, ymax are defined in  abc.ads
  hx:Real := (xmax-xmin)/Real(nx-1); -- x step
  hy:Real := (ymax-ymin)/Real(ny-1); -- y step
  hx2:Real := hx*hx;
  hy2:Real := hy*hy;
  hxy:Real := hx*hy;
  -- cs;    -- index of constant vector now in uc(ij)  A.put(ij-1,nrow)
  pdxx:Real_Matrix(1..nd,1..nd);
  pdxy:Real_Matrix(1..nd,1..nd);
  pdyy:Real_Matrix(1..nd,1..nd);
  pdx :Real_Matrix(1..nd,1..nd);
  pdy :Real_Matrix(1..nd,1..nd);
  pd  :Real_Matrix(1..nd,1..nd);
  md:integer := (nd+1)/2 ;  -- mid point of 1..nd, nd is odd
  error:Real;           -- sum of absolute exact-computed
  avgerror:Real;        -- average error
  maxerror:Real;        -- maximum cell error
  nrow : Integer := (nx-2)*(ny-2);
  package A is new Sparse(nrow);
  use A;

  package Real_IO is new Float_IO(Real);
  use Real_IO;

  procedure clear is -- just zero arrays
     ii : Integer;
  begin
    for i in 1..nx-2 loop -- solution, ut, coordinates
      -- solution coordinates OK for xmin+i*hx, ymin+j*hy
      for j in 1..ny-2 loop
        ii := i+(nx-2)*(j-1);
        us(ii) := 0.0;
        -- S uc(ii) := 0.0;
        if ifcheck>0 then -- all element coordinates
          uu(ii) := u(xmin+Real(i)*hx,ymin+Real(j)*hy);
        end if;
      end loop;
    end loop;
  end clear;

  procedure printus is -- print the computed solution
   ii:Integer;
  begin
    Put_Line("computed solution");
    Put_Line("us(x,y) ");
    for i in 1..nx-2 loop
      for j in 1..ny-2 loop
        ii := i+(nx-2)*(j-1);
        Put("i="&Integer'Image(i)&
                 ", j="&Integer'Image(j)&
            ", at x=");
        Put(xmin+Real(i)*hx, 5,3,0);
        Put(", y=");
        Put(ymin+Real(j)*hy, 5,3,0);
        Put(", us(i,j)=");
        Put(us(ii), 6,3,0);
        New_Line;
      end loop;
    end loop;
    New_Line;
  end printus;

  procedure check is  -- typically not known, instrumentation
    ii:Integer;
    uexact, err:Real;
  begin
    error := 0.0;
    maxerror := 0.0;
    for i in 1..nx-2 loop -- solution coordinates
      for j in 1..ny-2 loop
        ii := i+(nx-2)*(j-1);
        uexact := u(xmin+Real(i)*hx, ymin+Real(j)*hy);
        err := abs(us(ii)-uexact);
        error := error + err;
        if err>maxerror then  maxerror:=err; end if;
      end loop;
    end loop;
  end check;

  procedure printexact is -- typically not known
  begin
  Put_Line("exact solution");
    Put_Line("u(x,y) ");
    for i in 1..nx-2 loop
      for j in 1..ny-2 loop
        Put("i="&Integer'Image(i)&
                 ", j="&Integer'Image(j)&
            ", at x=");
        Put(xmin+Real(i)*hx, 5,3,0);
        Put(", y=");
        Put(ymin+Real(j)*hy, 5,3,0);
        Put(", u(i,j)=");
        Put(u(xmin+Real(i)*hx,ymin+Real(j)*hy), 6,3,0);
        New_Line;
      end loop;
    end loop;
    New_Line;
  end printexact;

  procedure printdiff is -- only when checking
    ii:Integer;
  begin
    Put_Line("difference exact-computed solution");
    Put_Line("u(x,y)-us(x,y) ");
    for i in 1..nx-2 loop
      for j in 1..ny-2 loop
        ii := i+(nx-2)*(j-1);
        Put("i="&Integer'Image(i)&
                 ", j="&Integer'Image(j)&
            ", at x=");
        Put(xmin+Real(i)*hx, 5,3,0);
        Put(", y=");
        Put(ymin+Real(j)*hy, 5,3,0);
        Put(", diff(i,j)=");
        Put(u(xmin+Real(i)*hx,ymin+Real(j)*hy)-us(ii), 6,3,0);
        New_Line;
      end loop;
    end loop;
    New_Line;
  end printdiff;

  procedure print_coef is -- for debugging
  begin
    Put_Line("pdxx     j=1      j=2      j=3      j=4      j=5  y");
    for i in 1..nd loop
      Put("i="&Integer'Image(i));
      for j in 1..nd loop
        Put(pdxx(i,j), 6,2,0);
      end loop;
      New_Line;
    end loop;
    Put_Line("x");
    New_Line;

    Put_Line("pdxy     j=1      j=2      j=3      j=4      j=5  y");
    for i in 1..nd loop
      Put("i="&Integer'Image(i));
      for j in 1..nd loop
        Put(pdxy(i,j), 6,2,0);
      end loop;
      New_Line;
    end loop;
    New_Line;

    Put_Line("pdyy     j=1      j=2      j=3      j=4      j=5  y");
    for i in 1..nd loop
      Put("i="&Integer'Image(i));
      for j in 1..nd loop
        Put(pdyy(i,j), 6,2,0);
      end loop;
      New_Line;
    end loop;
    New_Line;

    Put_Line("pdx      j=1      j=2      j=3      j=4      j=5  y");
    for i in 1..nd loop
      Put("i="&Integer'Image(i));
      for j in 1..nd loop
        Put(pdx(i,j), 6,2,0);
      end loop;
      New_Line;
    end loop;
    New_Line;

    Put_Line("pdy      j=1      j=2      j=3      j=4      j=5  y");
    for i in 1..nd loop
      Put("i="&Integer'Image(i));
      for j in 1..nd loop
        Put(pdy(i,j), 6,2,0);
      end loop;
      New_Line;
    end loop;
    New_Line;
  end print_coef;

  procedure build_coef(ix:Integer; iy:Integer) is
    -- ix and iy in 1..nd are 'point' to "deriv" and index
    a:Real_Vector(1..nd*nd);
  begin
    -- build discretization matrices
    for i in 1..nd loop
      for j in 1..nd loop
        pdxx(i,j):=0.0; -- coefficient of a1(x,y)
        pdxy(i,j):=0.0; -- coefficient of b1(x,y)
        pdyy(i,j):=0.0; -- coefficient of c1(x,y)
        pdx (i,j):=0.0; -- coefficient of d1(x,y)
        pdy (i,j):=0.0; -- coefficient of e1(x,y)
        pd  (i,j):=0.0; -- coefficient of f1(x,y)
      end loop;
    end loop;
    pd(ix,iy) := 1.0; -- u(i,j) term
    deriv(2, nd, ix, a); -- for pdxx
    for i in 1..nd loop
      pdxx(i,iy) := a(i)/hx2;
    end loop;

    deriv(2, nd, iy, a); -- for pdyy
    for j in 1..nd loop
      pdyy(ix,j) := a(j)/hy2;
    end loop;

    deriv(1, nd, ix, a); -- for pdx
    for i in 1..nd loop
      pdx(i,iy) := a(i)/hx;
    end loop;

    deriv(1, nd, iy, a); -- for pdy
    for j in 1..nd loop
      pdy(ix,j) := a(j)/hy;
    end loop;

    for i in 1..nd loop
      for j in 1..nd loop
        pdxy(i,j) := pdx(i,iy)*pdy(ix,j);
      end loop;
    end loop;
    if ifcheck>2 then print_coef; end if;
  end build_coef;

  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 1..nx-2 loop
      for j in 1..ny-2 loop
        ii := i+(nx-2)*(j-1);
        -- S sum := sum + ut(row,ii)*uu(ii); -- uu is check solution
        sum := sum + A.get(row-1,ii-1)*uu(ii); -- uu is check solution
      end loop;
    end loop;
    -- S sum := sum - uc(row);
    sum := sum - A.get(row-1,nrow);
    if abs(sum)>0.001 then
       Put_Line("row="&Integer'Image(row)&" is bad err="&
                Real'Image(sum));
    end if;
  end check_row;

  procedure build_row(i:Integer; j:Integer; bi:Integer; bj:integer) is
    -- i and j are solution coordinates, ir and jr are all element
    ii, ij, ir, jr, ia, ja: integer;
    p, a1c, b1c, c1c, d1c, e1c, f1c, cc : real;
  begin
    ir := i + 1; -- all element coordinates
    jr := j + 1;
    ii := i+(nx-2)*(j-1); -- row of solution, ut, matrix
    if ifcheck>2 then
      Put_Line("working row="&Integer'Image(ii)&
               ", i="&Integer'Image(i)&
               ", j="&Integer'Image(j)&
               ", bi="&Integer'Image(bi)&
               ", bj="&Integer'Image(bj)&
               ", ir="&Integer'Image(ir)&
               ", jr="&Integer'Image(jr));
    end if; -- check
    a1c := a1(xmin+Real(i)*hx,ymin+Real(j)*hy); -- yes, ir-1, jr-1
    b1c := b1(xmin+Real(i)*hx,ymin+Real(j)*hy);
    c1c := c1(xmin+Real(i)*hx,ymin+Real(j)*hy);
    d1c := d1(xmin+Real(i)*hx,ymin+Real(j)*hy);
    e1c := e1(xmin+Real(i)*hx,ymin+Real(j)*hy);
    f1c := f1(xmin+Real(i)*hx,ymin+Real(j)*hy);
    cc  :=  c(xmin+Real(i)*hx,ymin+Real(j)*hy);
    if ifcheck>1 then
      Put_Line("a1="&Real'Image(a1c));
      Put_Line("b1="&Real'Image(b1c));
      Put_Line("c1="&Real'Image(c1c));
      Put_Line("d1="&Real'Image(d1c));
      Put_Line("e1="&Real'Image(e1c));
      Put_Line("f1="&Real'Image(f1c));
      Put_Line("c ="&Real'Image(cc));
    end if; -- check
    for ki in 1..nd loop -- discretization
      for kj in 1..nd loop
        ia := ir-bi+ki-md; -- all element coordinates 1 and nx are boundary
        ja := jr-bj+kj-md; -- may be outside solution, thus boundary
        ij := (ia-1)+(nx-2)*(ja-2); -- solution, ut, coordinates
        p := a1c*pdxx(ki,kj) + b1c*pdxy(ki,kj) + c1c*pdyy(ki,kj)+
             d1c* pdx(ki,kj) + e1c* pdy(ki,kj) + f1c* pd(ki,kj);
        if ifcheck>1 then
          Put_Line("ia="&Integer'Image(ia)&", ja="&Integer'Image(ja)&
                   ", ij="&Integer'Image(ij)&", p="&Real'Image(p));
        end if;
        if ia=1 then  -- boundary value tests (correct for starting with 1)
          -- S uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy);
          A.add(ii-1, nrow, -p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy));
        elsif ia=nx then
          -- S uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy);
          A.add(Ii-1, nrow, -p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy));
        elsif ja=1 then
          -- S uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy);
          A.add(ii-1, nrow, -p * u(xmin+real(ia-1)*hx, ymin+real(ja-1)*hy));
        elsif ja=ny then
          -- S uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy);
          A.add(ii-1, nrow, -p * u(xmin+real(ia-1)*hx,ymin+real(ja-1)*hy));
        else -- non boundary case
          -- S ut(ii,ij) := ut(ii,ij) + p;
          A.add(ii-1,ij-1, P);
          Put_Line("ut(ii,"&Integer'Image(ij)&")="&
                   Real'Image(A.get(ii-1,ij-1)));  -- S ut(ii,ij)));
        end if;
      end loop;
    end loop;
    -- S uc(ii) := uc(ii) + cc;
    A.add(ii-1, nrow, cc);
    Put_Line("row ii="&Integer'Image(ii)&", uc(ii)="&Real'Image(
              A.get(ii-1, nrow))); -- S uc(ii)));
    check_row(ii);
  end build_row;

  procedure init_matrix is -- build the ut matrix
    i, j : integer;
  begin
    -- zero internal cells, the matrix does not include boundary cells
    -- S for i in 1..nx-2 loop    by  "is new"
    -- S  for j in 1..ny-2 loop
    -- S    for ii in 1..nx-2 loop
    -- S      for jj in 1..ny-2 loop
    -- S        ut(i+(nx-2)*(j-1),ii+(nx-2)*(jj-1)) := 0.0;
    -- S      end loop;
    -- S   end loop;
    -- S  end loop;
    -- S end loop;
    Put_Line("internal cells zeroed, do four corners");

    -- do four corners (special case, solution, ut, coordinates)
    i:=1;
    j:=1;
    build_coef(md-1, md-1);
    build_row(i, j, -1, -1);
    i:=1;
    j:=ny-2;
    build_coef(md-1, md+1);
    build_row(i, j, -1, 1);
    i:=nx-2;
    j:=1;
    build_coef(md+1, md-1);
    build_row(i, j, 1, -1);
    i:=nx-2;
    j:=ny-2;
    build_coef(md+1, md+1);
    build_row(i, j, 1, 1);
    Put_Line("four corners generated. do central cells");

    -- central cells in solution, ut, coordinates
    build_coef(md, md);
    for i in 2..nx-3 loop -- compute standard central cells
      for j in 2..ny-3 loop
        build_row(i, j, 0, 0);
      end loop;
    end loop;
    Put_Line("standard central cells generated, do sides");

    -- do four sides, one element in. special cases
    build_coef(md-1, md);
    for j in 2..ny-3 loop -- do i=1
      i := 1;
      build_row(i, j, -1, 0);
    end loop; -- on edge

    build_coef(md+1, md);
    for j in 2..ny-3 loop -- do i=nx-2
      i := nx-2;
      build_row(i, j, 1, 0);
    end loop; -- on edge

    build_coef(md, md-1);
    for i in 2..nx-3 loop -- do j=1
      j := 1;
      build_row(i, j, 0 , -1);
    end loop; -- on edge

    build_coef(md, Md+1);
    for i in 2..nx-3 loop -- do  j=ny-2
      j := ny-2;
      build_row(i, j, 0 , 1);
    end loop; -- on edge
    Put_Line("four sides generated");
  end init_matrix;

  procedure print_mat is
    sum : Real;
  begin
    Put_Line("initial matrix");
    for i in 1..(nx-2)*(ny-2) loop
      for J in 1..(nx-2)*(ny-2) loop
        Put(A.get(i-1,j-1),5,2,0); -- S ut(i,j),5,2,0);
        if (nx-2)*(j/(nx-2))=J and j/=(nx-2)*(ny-2) then New_Line; end if;
      end loop;
      Put(A.get(i-1, nrow),5,2,0); -- S uc(i),5,2,0);
      sum := 0.0;
      for ii in 1..(nx-2)*(ny-2) loop
        sum := sum + A.get(i-1,ii-1)*uu(ii); -- uu is check solution S ut(i,ii)
      end loop;
      Put("=");
      Put(sum,5,2,0);
      New_Line;
    end loop;
    New_Line;
  end print_mat;

begin -- sparse_abc
  Put_Line("sparse_abc.adb running");
  Put_Line("nx="&Integer'Image(nx)&
           ", ny="&Integer'Image(ny));
  Put_Line("xmin="&Real'Image(xmin)&
           ", xmax="&Real'Image(xmax));
  Put_Line("ymin="&Real'Image(ymin)&
           ", ymax="&Real'Image(ymax));
  Put_Line("hx="&Real'Image(hx)&
           ", hy="&Real'Image(hy));
  Put_Line("hx2="&Real'Image(hx2)&
           ", hy2="&Real'Image(hy2));
  clear;

  if Ifcheck>0 then Put_Line("md="&Integer'Image(md)); end if;

  init_matrix;
  print_mat;
  -- S us := simeq(ut, uc);  uc now inside ut
  A.Simeq(us); -- destroys, A

  printus;
  printexact;
  printdiff;
  if ifcheck>0 then
    check;
    avgerror := error/Real((nx-2)*(ny-2));
    Put_Line("average error = "&Real'Image(avgerror));
    Put_Line("maximum error = "&Real'Image(maxerror));
  end if;

end sparse_abc;