[BACK]Return to standard_matrix_inversion.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Schubert
File: [local] / OpenXM_contrib / PHC / Ada / Schubert / standard_matrix_inversion.adb (download)

Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:33 2000 UTC (23 years, 6 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD
Changes since 1.1: +0 -0 lines 

Import the second public release of PHCpack.

OKed by Jan Verschelde.

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;