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     ! 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;