with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Natural_Vectors;
with Standard_Floating_Vectors;
with Standard_Complex_Vectors;
with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;

package body Standard_Matrix_Inversion is

-- ALGORITHM :
--   Solves repeatedly A*x = b, with b the i-th standard basis vector.

  function Inverse ( m : Standard_Floating_Matrices.Matrix )
                   return Standard_Floating_Matrices.Matrix is

    n : constant natural := m'last(1);
    wrk : Standard_Floating_Matrices.Matrix(m'range(1),m'range(2)) := m;
    res : Standard_Floating_Matrices.Matrix(m'range(1),m'range(2));
    piv : Standard_Natural_Vectors.Vector(m'range(1));
    rhs : Standard_Floating_Vectors.Vector(m'range(1));
    inf : natural;

  begin
    for i in piv'range loop
      piv(i) := i;
    end loop;
    lufac(wrk,n,piv,inf);
    if inf = 0
     then for j in m'range(2) loop
            rhs := (rhs'range => 0.0);
            rhs(j) := 1.0;
            lusolve(wrk,n,piv,rhs);
            for i in m'range(1) loop
              res(i,j) := rhs(i);
            end loop;
          end loop;
    end if;
    return res;
  end Inverse;

  function Inverse ( m : Standard_Complex_Matrices.Matrix )
                   return Standard_Complex_Matrices.Matrix is

    n : constant natural := m'last(1);
    wrk : Standard_Complex_Matrices.Matrix(m'range(1),m'range(2)) := m;
    res : Standard_Complex_Matrices.Matrix(m'range(1),m'range(2));
    piv : Standard_Natural_Vectors.Vector(m'range(1));
    rhs : Standard_Complex_Vectors.Vector(m'range(1));
    inf : natural;

  begin
    for i in piv'range loop
      piv(i) := i;
    end loop;
    lufac(wrk,n,piv,inf);
    if inf = 0
     then for j in m'range(2) loop
            rhs := (rhs'range => Create(0.0));
            rhs(j) := Create(1.0);
            lusolve(wrk,n,piv,rhs);
            for i in m'range(1) loop
              res(i,j) := rhs(i);
            end loop;
          end loop;
    end if;
    return res;
  end Inverse;

end Standard_Matrix_Inversion;