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

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;