Return to maximal_minors.adb CVS log

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

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

1.1       maekawa     1: with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
                      2: with Standard_Natural_Vectors;           use Standard_Natural_Vectors;
                      3: with Standard_Natural_Vectors_io;        use Standard_Natural_Vectors_io;
                      4: with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
                      5:
                      6: procedure Maximal_Minors ( file : in file_type;
                      7:                            n,d : in natural; mat : in Matrix;
                      8:                            min,max : out double_float ) is
                      9:
                     10:   function Determinant
                     11:               ( mat : Matrix; rows : Standard_Natural_Vectors.Vector )
                     12:               return double_float is
                     13:
                     14:   -- DESCRIPTION :
                     15:   --   Computes the determinant of the matrix obtained by selecting rows.
                     16:
                     17:     res : double_float := 1.0;
                     18:     sqm : Matrix(rows'range,rows'range);
                     19:     piv : Standard_Natural_Vectors.Vector(rows'range);
                     20:     inf : natural;
                     21:
                     22:   begin
                     23:     for i in rows'range loop
                     24:       piv(i) := i;
                     25:       for j in rows'range loop
                     26:         sqm(i,j) := mat(rows(i),j);
                     27:       end loop;
                     28:     end loop;
                     29:     lufac(sqm,rows'last,piv,inf);
                     30:     for i in rows'range loop
                     31:       res := res*sqm(i,i);
                     32:     end loop;
                     33:     for i in piv'range loop
                     34:       if piv(i) > i
                     35:        then res := -res;
                     36:       end if;
                     37:     end loop;
                     38:     return res;
                     39:   end Determinant;
                     40:
                     41:   procedure Main is
                     42:
                     43:     rows : Standard_Natural_Vectors.Vector(1..d);
                     44:     first : boolean := true;
                     45:     mindet,maxdet : double_float;
                     46:
                     47:     procedure Select_Rows ( k,start : in natural ) is
                     48:
                     49:       det : double_float;
                     50:
                     51:     begin
                     52:       if k > d
                     53:        then det := Determinant(mat,rows);
                     54:             put(file,"Minor "); put(file,rows); put(file," equals ");
                     55:             put(file,det); new_line(file);
                     56:             det := abs(det);
                     57:             if first
                     58:              then mindet := det; maxdet := det; first := false;
                     59:              else if det > maxdet
                     60:                    then maxdet := det;
                     61:                    elsif det < mindet
                     62:                        then mindet := det;
                     63:                   end if;
                     64:             end if;
                     65:        else for j in start..n loop
                     66:               rows(k) := j;
                     67:               Select_Rows(k+1,j+1);
                     68:             end loop;
                     69:       end if;
                     70:     end Select_Rows;
                     71:
                     72:   begin
                     73:     Select_Rows(1,1);
                     74:     put(file,"Min : ");       put(file,mindet,3,3,3);
                     75:     put(file,"  Max : ");     put(file,maxdet,3,3,3);
                     76:     put(file,"  Max/Min : "); put(file,maxdet/mindet,3,3,3); new_line(file);
                     77:     min := mindet; max := maxdet;
                     78:   end;
                     79:
                     80: begin
                     81:   Main;
                     82: end Maximal_Minors;