-- pde3b_eq.adb  3D non equal h boundary value PDE  non iterative    
--  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  
--     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  
-- d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  
--     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  
-- d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  
--     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  
-- d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  
--                                                                 
-- solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)              
-- by second order method on cube -1,-1,-1 to 1.2,1,1              

-- second order numerical difference                              
-- d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                            
-- 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +                    
-- 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +                    
-- 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)             
--  = c(x,y,z) = 12x + 10y - 8z                                   

-- solve system of linear equations for                           
-- u(i-1,j,k)/hx^2 + u(i+1,j,k)/hx^2 +                            
-- u(i,j-1,k)/hy^2 + u(i,j+1,k)/hy^2 +                            
-- u(i,j,k-1)/hz^2 + u(i,j,k+1)/hz^2 -                            
-- (2/hx^2 + 2/hy^2 + 2/hz^2)*u(i,j,k) = c(x,y,z)                 
-- for i=1,nx-2 j=1,ny-2 k=1,nz-2 while substituting              
-- u(i,j,k) boundary values i=0 or nx-1, j=0 or ny-1, k=0 or nz-1 
-- x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz                          

-- create two dimensional array with                              
-- (nx-2)*(ny-2)*(nz-2) columns, one for each equation            
-- and (nx-2)*(ny-2)*(nz-2)+1 rows for variables plus constant    
-- Initialize the matrix, see init_matrix in the source code      
-- Then solve the system of linear equations to get the solution. 


with ada.text_io; use ada.text_io;

procedure pde3b_eq is
  type real_matrix is array(integer range <>, integer range <>)
                      of long_float;
  type real_vector is array(integer range <>) of long_float;
  us : real_vector(0..566); -- solution being iterated
  ut : real_matrix(0..566, 0..567); -- temporary for partial results
                -- nx-1 by ny-1 by nz-1 internal, non boundary cells
  nx : integer := 10; -- one less than 'C' when Ada subscript starts at zero
  ny : integer := 10;  
  nz : integer := 8;
  xmin : long_float := -1.0;
  ymin : long_float := -1.0;
  zmin : long_float := -1.0;
  xmax : long_float :=  1.2;
  ymax : long_float :=  1.0;
  zmax : long_float :=  1.0;
  hx   : long_float;         -- x step = (xmax-xmin)/nx 
  hy   : long_float;         -- y step = (ymax-ymin)/ny 
  hz   : long_float;         -- z step = (zmax-zmin)/nz 
  error : long_float;        -- sum of absolute exact-computed         
  avg_error : long_float;    -- average error                          
  max_error : long_float;    -- maximum cell error
  
  function u(x:long_float; y:long_float; z:long_float) return long_float is
           -- solution for boundary and check 
  begin
    return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
           4.0*z*z*x + 5.0*y +6.0;
  end u;

  function c(x:long_float; y:long_float; z:long_float) return long_float is
           -- from solve
  begin
    return 20.0*x + 10.0*y + 6.0*z;
  end c;

  procedure printexact is 
  begin
    new_line;
    put_line("exact solution");
    for i in 0..nx loop
      for j in 0..ny loop
        for k in 0..nz loop
          put_line("u("&integer'image(i)&","&integer'image(j)&
                   ","&integer'image(k)&")="&long_float'image(
                   u(xmin+long_float(i)*hx,ymin+long_float(j)*hy,zmin+long_float(k)*hz)));
        end loop;
      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
        for k in 1..nz-1 loop
          put_line("us("&integer'image(i)&","&integer'image(j)&
                   ","&integer'image(k)&")="&
                   long_float'image(us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1))));
        end loop;
      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
        for k in 1..nz-1 loop
          put_line("diff("&integer'image(i)&","&integer'image(j)&
                        ","&integer'image(k)&")="&long_float'image(
          u(xmin+long_float(i)*hx,ymin+long_float(j)*hy,zmin+long_float(k)*hz)-
          us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1))));
        end loop;
      end loop;
    end loop; 
  end printdiff;

  procedure check is
    zexact, err : long_float;
  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(xmin+long_float(i)*hx,
                      ymin+long_float(j)*hy,
                      zmin+long_float(k)*hz);
          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, cs : integer;
    hx2, hy2, hz2 : long_float;
  begin
    hx2 := hx*hx;
    hy2 := hy*hy;
    hz2 := hz*hz;
    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;
        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); -- column of matrix 
          ut(ii,ii) := -(2.0/hx2 + 2.0/hy2 + 2.0/hz2);
          ut(ii,cs) := c(xmin+long_float(i)*hx, ymin+long_float(j)*hy, zmin+long_float(k)*hz);
          if (i+1)<(nx-1) then
            ut(ii,((i-1)+1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)) := 1.0/hx2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i+1)*hx,
                                       ymin+long_float(j)*hy,
                                       zmin+long_float(k)*hz)/hx2;
          end if;
          if (i-1)/=0 then
            ut(ii,((i-1)-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)) := 1.0/hx2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i-1)*hx,
                                       ymin+long_float(j)*hy,
                                       zmin+long_float(k)*hz)/hx2;
          end if;
          if (j+1)<(ny-1) then
            ut(ii,(i-1)+(nx-1)*((j-1)+1)+(nx-1)*(ny-1)*(k-1)) := 1.0/hy2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i)*hx,
                                       ymin+long_float(j+1)*hy,
                                       zmin+long_float(k)*hz)/hy2;
          end if;
          if (j-1)/=0 then
            ut(ii,(i-1)+(nx-1)*((j-1)-1)+(nx-1)*(ny-1)*(k-1)) := 1.0/hy2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i)*hx,
                                       ymin+long_float(j-1)*hy,
                                       zmin+long_float(k)*hz)/hy2;
          end if;
          if (k+1)<(nz-1) then
            ut(ii,(i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*((k-1)+1)) := 1.0/hz2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i)*hx,
                                       ymin+long_float(j)*hy,
                                       zmin+long_float(k+1)*hz)/hz2;
          end if;
          if (k-1)/=0 then 
            ut(ii,(i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*((k-1)-1)) := 1.0/hz2;
          else
            ut(ii,cs) := ut(ii,cs) - u(xmin+long_float(i)*hx,
                                       ymin+long_float(j)*hy,
                                       zmin+long_float(k-1)*hz)/hz2;
          end if;
        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
    row : array(0..n)of integer;   -- row interchange indices
    hold , i_pivot : integer;      -- pivot indices
    pivot : long_float;            -- pivot element value
    abs_pivot : long_float;        -- abs of pivot element
    norm1 : long_float := 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 < long_float'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  -- pde3b_eq
  
  put_line("pde3b_eq.adb running");
  hx := (xmax-xmin)/long_float(nx); -- step size in x 
  hy := (ymax-ymin)/long_float(ny); -- step size in y 
  hz := (zmax-zmin)/long_float(nz); -- step size in z 
  put_line("xmin="&long_float'image(xmin)&", xmax="&long_float'image(xmax)&
           ", nx="&integer'image(nx)&", hx="&long_float'image(hx));
  put_line("ymin="&long_float'image(ymin)&", ymax="&long_float'image(ymax)&
           ", ny="&integer'image(ny)&", hy="&long_float'image(hy));
  put_line("zmin="&long_float'image(zmin)&", zmax="&long_float'image(zmax)&
           ", nz="&integer'image(nz)&", hz="&long_float'image(hz));
  
  init_matrix;
  put_line("matrix initialized");
  put_line("initial matrix, left side, upper");
  for i in 0..3 loop
    for j in 0..3 loop
      put(" "&long_float'image(ut(i,j)));
    end loop;
    new_line;
  end loop;
  new_line;
  
  put_line("initial matrix, right side, upper");
  for i in 0..3 loop
    for j in 564..567 loop
      put(" "&long_float'image(ut(i,j)));
    end loop;
    new_line;
  end loop;
  new_line;
  
  put_line("initial matrix, right side, lower");
  for i in 563..566 loop
    for j in 564..567 loop
      put(" "&long_float'image(ut(i,j)));
    end loop;
    new_line;
  end loop;
  new_line;
  
  simeq(566, ut, us); -- tailored version, returned us
  
  printexact;
  printus;
  printdiff;
  check;
  avg_error := error/(long_float(nx-1)*long_float(ny-1)*long_float(nz-1));
  put_line("error="&long_float'image(error)&" avg="&long_float'image(avg_error)&
           ", max="&long_float'image(max_error));
end pde3b_eq;