Return to global_dynamic_triangulation.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift

Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/global_dynamic_triangulation.adb, Revision 1.1

1.1     ! maekawa     1: with Standard_Random_Numbers;            use Standard_Random_Numbers;
        !             2: with Standard_Integer_VecVecs;           use Standard_Integer_VecVecs;
        !             3: with Standard_Integer_Matrices;          use Standard_Integer_Matrices;
        !             4: with Standard_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
        !             5: with Vertices;
        !             6:
        !             7: package body Global_Dynamic_Triangulation is
        !             8:
        !             9:   function Ordered_Initial_Simplex
        !            10:                ( n : in natural; pts : in List ) return List is
        !            11:
        !            12:   -- DESCRIPTION :
        !            13:   --   The given list pts has at least n+1 points.
        !            14:   --   The order of the points should be taken into account when
        !            15:   --   constructing the initial cell.
        !            16:   --   Returns the list of points which span the initial simplex.
        !            17:
        !            18:   -- ALGORITHM :
        !            19:   --   Determines the rank of the matrix, defined by the first n+1 points.
        !            20:   --   If necessary, more extremal points will be sought.
        !            21:
        !            22:     tmp,res,res_last : List;
        !            23:     mat : matrix(1..n,1..n);
        !            24:     sh,pt : Link_to_Vector;
        !            25:     rnk : natural;
        !            26:
        !            27:   begin
        !            28:     sh := Head_Of(pts);
        !            29:     tmp := Tail_Of(pts);
        !            30:     for i in 1..n loop
        !            31:       pt := Head_Of(tmp);
        !            32:       for j in pt'range loop
        !            33:         mat(i,j) := pt(j) - sh(j);
        !            34:       end loop;
        !            35:       tmp := Tail_Of(tmp);
        !            36:     end loop;
        !            37:     rnk := Rank(mat);
        !            38:     if rnk = n
        !            39:      then tmp := pts;
        !            40:           for i in 1..n+1 loop
        !            41:             pt := Head_Of(tmp);
        !            42:             Append(res,res_last,pt.all);
        !            43:             tmp := Tail_Of(tmp);
        !            44:           end loop;
        !            45:     end if;
        !            46:     return res;
        !            47:   end Ordered_Initial_Simplex;
        !            48:
        !            49:   procedure Initial_Simplex
        !            50:                ( pts : in List; order : in boolean;
        !            51:                  s : out Simplex; rest : in out List ) is
        !            52:
        !            53:     n : natural;
        !            54:     res,tmp,rest_last : List;
        !            55:     pt : Link_to_Vector;
        !            56:
        !            57:   begin
        !            58:     if Is_Null(pts)
        !            59:      then
        !            60:        s := Null_Simplex;
        !            61:      else
        !            62:        n := Head_Of(pts)'last;
        !            63:        if Length_Of(pts) < n+1
        !            64:         then
        !            65:           s := Null_Simplex;
        !            66:         else
        !            67:           if order
        !            68:            then res := Ordered_Initial_Simplex(n,pts);
        !            69:                 if Is_Null(res)
        !            70:                  then res := Vertices.Extremal_Points(n,pts);
        !            71:                 end if;
        !            72:            else res := Vertices.Extremal_Points(n,pts);
        !            73:           end if;
        !            74:           if Length_Of(res) = n+1
        !            75:            then
        !            76:              declare
        !            77:                points : VecVec(1..n+1);
        !            78:              begin
        !            79:                tmp := res;
        !            80:                for k in points'range loop
        !            81:                  pt := Head_Of(tmp);
        !            82:                  points(k) := new vector(1..n+1);
        !            83:                  points(k)(1..n) := pt.all;
        !            84:                  points(k)(n+1) := 0;
        !            85:                  tmp := Tail_Of(tmp);
        !            86:                end loop;
        !            87:                s := Create(points);
        !            88:              end;
        !            89:             -- construct the list with the rest of the points :
        !            90:              tmp := pts;
        !            91:              while not Is_Null(tmp) loop
        !            92:                pt := Head_Of(tmp);
        !            93:                if not Is_In(res,pt.all)
        !            94:                 then Append(rest,rest_last,pt.all);
        !            95:                end if;
        !            96:                tmp := Tail_Of(tmp);
        !            97:              end loop;
        !            98:            else s := Null_Simplex;
        !            99:           end if;
        !           100:           Deep_Clear(res);
        !           101:        end if;
        !           102:     end if;
        !           103:   end Initial_Simplex;
        !           104:
        !           105: -- COMPUTING AN EXTREMAL POINT :
        !           106:
        !           107:   function Max_Extreme ( l : List; k : natural ) return Link_to_Vector is
        !           108:
        !           109:     res : Link_to_Vector := Head_Of(l);
        !           110:     mx : integer := res(k);
        !           111:     tmp : List := Tail_Of(l);
        !           112:     pt : Link_to_Vector;
        !           113:
        !           114:   begin
        !           115:     while not Is_Null(tmp) loop
        !           116:       pt := Head_Of(tmp);
        !           117:       if pt(k) > mx
        !           118:        then mx := pt(k);
        !           119:             res := pt;
        !           120:       end if;
        !           121:       tmp := Tail_Of(tmp);
        !           122:     end loop;
        !           123:     return res;
        !           124:   end Max_Extreme;
        !           125:
        !           126:   function Max_Extreme ( l : List; weights : vector ) return Link_to_Vector is
        !           127:
        !           128:     res : Link_to_Vector := Head_Of(l);
        !           129:     tmp : List := Tail_Of(l);
        !           130:     mxs : integer := res.all*weights;
        !           131:     sp : integer;
        !           132:     pt : Link_to_Vector;
        !           133:
        !           134:   begin
        !           135:     while not Is_Null(tmp) loop
        !           136:       pt := Head_Of(tmp);
        !           137:       sp := pt.all*weights;
        !           138:       if sp > mxs
        !           139:        then mxs := sp;
        !           140:             res := pt;
        !           141:       end if;
        !           142:       tmp := Tail_Of(tmp);
        !           143:     end loop;
        !           144:     return res;
        !           145:   end Max_Extreme;
        !           146:
        !           147:   function Max_Extreme ( l : List; n : natural; low,upp : integer )
        !           148:                        return Link_to_Vector is
        !           149:
        !           150:     w : Vector(1..n);
        !           151:
        !           152:   begin
        !           153:     for k in w'range loop
        !           154:       w(k) := Random(low,upp);
        !           155:     end loop;
        !           156:     return Max_Extreme(l,w);
        !           157:   end Max_Extreme;
        !           158:
        !           159: end Global_Dynamic_Triangulation;