Return to evaluated_minors.adb CVS log

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

Annotation of OpenXM_contrib/PHC/Ada/Schubert/evaluated_minors.adb, Revision 1.1.1.1

1.1       maekawa     1: with Standard_Natural_Vectors;           use Standard_Natural_Vectors;
                      2: with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
                      3: with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
                      4:
                      5: package body Evaluated_Minors is
                      6:
                      7:   function Determinant ( m : Standard_Floating_Matrices.Matrix )
                      8:                        return double_float is
                      9:
                     10:     res : double_float;
                     11:     wrk : Standard_Floating_Matrices.Matrix(m'range(1),m'range(2));
                     12:     piv : Standard_Natural_Vectors.Vector(m'range(1));
                     13:     inf : natural;
                     14:
                     15:   begin
                     16:     for i in m'range(1) loop
                     17:       piv(i) := i;
                     18:       for j in m'range(2) loop
                     19:         wrk(i,j) := m(i,j);
                     20:       end loop;
                     21:     end loop;
                     22:     lufac(wrk,m'last(1),piv,inf);
                     23:     if inf /= 0
                     24:      then res := 0.0;
                     25:      else res := 1.0;
                     26:           for i in m'range(1) loop
                     27:             res := res*wrk(i,i);
                     28:           end loop;
                     29:           for i in piv'range loop
                     30:             if piv(i) > i
                     31:              then res := -res;
                     32:             end if;
                     33:           end loop;
                     34:     end if;
                     35:     return res;
                     36:   end Determinant;
                     37:
                     38:   function Determinant ( m : Standard_Floating_Matrices.Matrix; b : Bracket )
                     39:                        return double_float is
                     40:
                     41:     res : double_float;
                     42:     sqm : Standard_Floating_Matrices.Matrix(b'range,b'range);
                     43:     piv : Standard_Natural_Vectors.Vector(b'range);
                     44:     inf : natural;
                     45:
                     46:   begin
                     47:     for i in b'range loop
                     48:       piv(i) := i;
                     49:       for j in b'range loop
                     50:         sqm(i,j) := m(b(i),j);
                     51:       end loop;
                     52:     end loop;
                     53:     lufac(sqm,b'last,piv,inf);
                     54:     if inf /= 0
                     55:      then res := 0.0;
                     56:      else res := 1.0;
                     57:           for i in b'range loop
                     58:             res := res*sqm(i,i);
                     59:           end loop;
                     60:           for i in piv'range loop
                     61:             if piv(i) > i
                     62:              then res := -res;
                     63:             end if;
                     64:           end loop;
                     65:     end if;
                     66:     return res;
                     67:   end Determinant;
                     68:
                     69:   function Determinant ( m : Standard_Complex_Matrices.Matrix )
                     70:                        return Complex_Number is
                     71:
                     72:     res : Complex_Number;
                     73:     wrk : Standard_Complex_Matrices.Matrix(m'range(1),m'range(2));
                     74:     piv : Standard_Natural_Vectors.Vector(m'range(1));
                     75:     inf : natural;
                     76:
                     77:   begin
                     78:     for i in m'range(1) loop
                     79:       piv(i) := i;
                     80:       for j in m'range(2) loop
                     81:         wrk(i,j) := m(i,j);
                     82:       end loop;
                     83:     end loop;
                     84:     lufac(wrk,wrk'last(1),piv,inf);
                     85:     if inf /= 0
                     86:      then res := Create(0.0);
                     87:      else res := Create(1.0);
                     88:           for i in wrk'range(1) loop
                     89:             res := res*wrk(i,i);
                     90:           end loop;
                     91:           for i in piv'range loop
                     92:             if piv(i) > i
                     93:              then res := -res;
                     94:             end if;
                     95:           end loop;
                     96:     end if;
                     97:     return res;
                     98:   end Determinant;
                     99:
                    100:   function Determinant ( m : Standard_Complex_Matrices.Matrix; b : Bracket )
                    101:                        return Complex_Number is
                    102:
                    103:     res : Complex_Number;
                    104:     sqm : Standard_Complex_Matrices.Matrix(b'range,b'range);
                    105:     piv : Standard_Natural_Vectors.Vector(b'range);
                    106:     inf : natural;
                    107:
                    108:   begin
                    109:     for i in b'range loop
                    110:       piv(i) := i;
                    111:       for j in b'range loop
                    112:         sqm(i,j) := m(b(i),j);
                    113:       end loop;
                    114:     end loop;
                    115:     lufac(sqm,b'last,piv,inf);
                    116:     if inf /= 0
                    117:      then res := Create(0.0);
                    118:      else res := Create(1.0);
                    119:           for i in sqm'range(1) loop
                    120:             res := res*sqm(i,i);
                    121:           end loop;
                    122:           for i in piv'range loop
                    123:             if piv(i) > i
                    124:              then res := -res;
                    125:             end if;
                    126:           end loop;
                    127:     end if;
                    128:     return res;
                    129:   end Determinant;
                    130:
                    131: end Evaluated_Minors;