-- crout.adb  solve simultaneous equation  A x = c  by Crout method
with ada.text_io; use ada.text_io;

procedure crout is
   type real_matrix is array(integer range <>,integer range <>)of float;
   type real_vector is array(integer range <>)of float;
   n:integer:=4;
   a:real_matrix(1..n,1..n):=((1.5,2.0,3.0,3.5),
                              (2.0,2.0,3.0,4.0),
                              (3.0,3.0,3.0,4.0),
                              (4.0,4.0,4.0,4.0));
   c:real_vector(1..n):=(1.0,2.0,3.0,4.0);
   x:real_vector(1..n);
   sum:float;
   diag:float;
   i:integer;
   j:integer;

   -- for checking ...
   aa:real_matrix(1..n,1..n):=a;
   cc:real_vector(1..n):=c;
   package float_io is new ada.text_io.float_io(float);
   use float_io;

begin
   put_line("Testing Crout method of solving simultaneous equations");
   default_aft:=4;
   x(1) := 0.0;

   for m in 1..n loop
       diag := a(m,m);

       -- below diagonal computations
       j := m;
       for i in m..n loop
         sum := 0.0;
         for k in 1..j-1 loop
           sum := sum + a(i,k)*a(k,j);
         end loop; -- k
         a(i,j) := a(i,j) - sum;
       end loop; -- i

       -- above diagonal computations
       i := m;
       for j in m+1..n loop
         sum := 0.0;
         for k in 1..i-1 loop
           sum := sum + a(i,k)*a(k,j);
         end loop; -- k
         a(i,j) := (a(i,j)-sum)/a(i,i);
       end loop; -- j
       sum := 0.0;
       for k in 1..i-1 loop
         sum := sum + a(i,k)*c(k);
       end loop; -- k
       c(i) := (c(i)-sum)/a(i,i);
   end loop; -- m

   -- back substitute
   for i in reverse 1..n loop
     sum := 0.0;
     for k in i+1..n loop
       sum := sum + a(i,k)*x(k);
     end loop; -- k
     x(i) := c(i) - sum;
   end loop; -- i

   -- print results
   put_line("initial matrix, initial c vector");
   for i in 1..n loop
     for j in 1..n loop
       put(aa(i,j));put("  ");
     end loop; -- j
     put(cc(i));
     new_line;
   end loop; -- i

   put_line("modified matrix, c vector, x vector");
   for i in 1..n loop
     for j in 1..n loop
       put(a(i,j));put("  ");
     end loop; -- j
     put(c(i));put("  ");
     put(x(i));
     new_line;
   end loop; -- i

   -- check results
   put_line("check");
   for i in 1..n loop
   sum := 0.0;
     for j in 1..n loop
       sum := sum + aa(i,j)*x(j);
     end loop; -- j
   put(sum);put(" should be ");put(cc(i));new_line;
   end loop; -- i
end crout;