Return to solve_pieri_leaves.adb CVS log

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

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

1.1       maekawa     1: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
                      2: with Standard_Complex_Vectors_io;        use Standard_Complex_Vectors_io;
                      3: with Standard_Natural_Vectors;
                      4: with Standard_Complex_Vectors;
                      5: with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
                      6:
                      7: function Solve_Pieri_Leaves
                      8:            ( file : file_type; b1,b2 : Bracket; m : Matrix ) return Matrix is
                      9:
                     10:   res : Matrix(m'range(1),b1'range);
                     11:   n : constant natural := m'last(1);
                     12:   dim : constant natural := m'last(2)-1;  -- dimension of linear system
                     13:   equ : Bracket(1..dim);                  -- equations : res(equ(i),*) = 0
                     14:   sys : Matrix(1..dim,1..dim);
                     15:   rhs : Standard_Complex_Vectors.Vector(1..dim);
                     16:   gen : Standard_Complex_Vectors.Vector(1..n);
                     17:
                     18:   procedure Initialize_Solution_Plane is
                     19:
                     20:   -- DESCRIPTION :
                     21:   --   Fills in the patterns of the resulting p-plane and sets up the
                     22:   --   equations that the resulting plane has to satisfy.
                     23:   --   These equations define the space C.
                     24:
                     25:     cnt : natural := 0;
                     26:     zero_row : boolean;
                     27:
                     28:   begin
                     29:     for i in 1..n loop
                     30:       zero_row := true;
                     31:       for j in b1'range loop
                     32:         if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
                     33:          then null;
                     34:          else zero_row := false;
                     35:         end if;
                     36:       end loop;
                     37:       if zero_row
                     38:        then cnt := cnt+1;
                     39:             equ(cnt) := i;
                     40:       end if;
                     41:     end loop;
                     42:   end Initialize_Solution_Plane;
                     43:
                     44:   procedure Solve_Linear_System is
                     45:
                     46:   -- DESCRIPTION :
                     47:   --   Solves the linear system that determines the solution plane.
                     48:   --   The solution to the linear system gives the coefficients in the
                     49:   --   linear combination of the columns of the given m-plane that
                     50:   --   determines the generator for the line of intersection with the
                     51:   --   space C and the given m-plane.
                     52:
                     53:     ipvt : Standard_Natural_Vectors.Vector(1..dim);
                     54:     info : natural;
                     55:    -- orgsys : Standard_Complex_Matrices.Matrix(1..dim,1..dim);
                     56:    -- orgrhs,resi : Standard_Complex_Vectors.Vector(1..dim);
                     57:     use Standard_Complex_Vectors,Standard_Complex_Matrices;
                     58:
                     59:   begin
                     60:     for j in 1..dim loop
                     61:       for i in equ'range loop
                     62:         sys(i,j) := m(equ(i),j);
                     63:       end loop;
                     64:     end loop;
                     65:    -- orgsys := sys;
                     66:     for i in 1..dim loop
                     67:       rhs(i) := m(equ(i),m'last(2));
                     68:     end loop;
                     69:    -- orgrhs := rhs;
                     70:     lufac(sys,dim,ipvt,info);
                     71:     lusolve(sys,dim,ipvt,rhs);
                     72:     for i in 1..n loop       -- linear combination of the columns of m-plane
                     73:       gen(i) := -m(i,m'last(2));
                     74:       for j in rhs'range loop
                     75:         gen(i) := gen(i) + rhs(j)*m(i,j);
                     76:       end loop;
                     77:     end loop;
                     78:    -- put_line(file,"The solution to the linear system : ");
                     79:    -- put_line(file,rhs); new_line(file);
                     80:    -- resi := orgrhs - orgsys*rhs;
                     81:    -- put_line(file,"The residual vector of the linear system : ");
                     82:    -- put(file,resi,2); new_line(file);
                     83:    -- put_line(file,"The generator of the intersection : ");
                     84:    -- put_line(file,gen); new_line(file);
                     85:   end Solve_Linear_System;
                     86:
                     87:   procedure Determine_Solution_Plane is
                     88:
                     89:   -- DESCRIPTION :
                     90:   --   Determines the solution plane from the generator of the intersection
                     91:   --   of the space C with the given m-plane.
                     92:
                     93:   begin
                     94:     for i in 1..n loop
                     95:       for j in b1'range loop
                     96:         if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
                     97:          then res(i,j) := Create(0.0);
                     98:          else res(i,j) := gen(i);
                     99:         end if;
                    100:       end loop;
                    101:     end loop;
                    102:   end Determine_Solution_Plane;
                    103:
                    104: begin
                    105:   Initialize_Solution_Plane;
                    106:   Solve_Linear_System;
                    107:   Determine_Solution_Plane;
                    108:   return res;
                    109: end Solve_Pieri_Leaves;