Return to standard_matrix_inversion.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Schubert

Annotation of OpenXM_contrib/PHC/Ada/Schubert/standard_matrix_inversion.adb, Revision 1.1

1.1     ! maekawa     1: with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
        !             2: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
        !             3: with Standard_Natural_Vectors;
        !             4: with Standard_Floating_Vectors;
        !             5: with Standard_Complex_Vectors;
        !             6: with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
        !             7: with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
        !             8:
        !             9: package body Standard_Matrix_Inversion is
        !            10:
        !            11: -- ALGORITHM :
        !            12: --   Solves repeatedly A*x = b, with b the i-th standard basis vector.
        !            13:
        !            14:   function Inverse ( m : Standard_Floating_Matrices.Matrix )
        !            15:                    return Standard_Floating_Matrices.Matrix is
        !            16:
        !            17:     n : constant natural := m'last(1);
        !            18:     wrk : Standard_Floating_Matrices.Matrix(m'range(1),m'range(2)) := m;
        !            19:     res : Standard_Floating_Matrices.Matrix(m'range(1),m'range(2));
        !            20:     piv : Standard_Natural_Vectors.Vector(m'range(1));
        !            21:     rhs : Standard_Floating_Vectors.Vector(m'range(1));
        !            22:     inf : natural;
        !            23:
        !            24:   begin
        !            25:     for i in piv'range loop
        !            26:       piv(i) := i;
        !            27:     end loop;
        !            28:     lufac(wrk,n,piv,inf);
        !            29:     if inf = 0
        !            30:      then for j in m'range(2) loop
        !            31:             rhs := (rhs'range => 0.0);
        !            32:             rhs(j) := 1.0;
        !            33:             lusolve(wrk,n,piv,rhs);
        !            34:             for i in m'range(1) loop
        !            35:               res(i,j) := rhs(i);
        !            36:             end loop;
        !            37:           end loop;
        !            38:     end if;
        !            39:     return res;
        !            40:   end Inverse;
        !            41:
        !            42:   function Inverse ( m : Standard_Complex_Matrices.Matrix )
        !            43:                    return Standard_Complex_Matrices.Matrix is
        !            44:
        !            45:     n : constant natural := m'last(1);
        !            46:     wrk : Standard_Complex_Matrices.Matrix(m'range(1),m'range(2)) := m;
        !            47:     res : Standard_Complex_Matrices.Matrix(m'range(1),m'range(2));
        !            48:     piv : Standard_Natural_Vectors.Vector(m'range(1));
        !            49:     rhs : Standard_Complex_Vectors.Vector(m'range(1));
        !            50:     inf : natural;
        !            51:
        !            52:   begin
        !            53:     for i in piv'range loop
        !            54:       piv(i) := i;
        !            55:     end loop;
        !            56:     lufac(wrk,n,piv,inf);
        !            57:     if inf = 0
        !            58:      then for j in m'range(2) loop
        !            59:             rhs := (rhs'range => Create(0.0));
        !            60:             rhs(j) := Create(1.0);
        !            61:             lusolve(wrk,n,piv,rhs);
        !            62:             for i in m'range(1) loop
        !            63:               res(i,j) := rhs(i);
        !            64:             end loop;
        !            65:           end loop;
        !            66:     end if;
        !            67:     return res;
        !            68:   end Inverse;
        !            69:
        !            70: end Standard_Matrix_Inversion;