Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/initial_mixed_cell.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
2: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
4: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
5: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
6: with Vertices; use Vertices;
7: with Frequency_Graph; use Frequency_Graph;
8:
9: procedure Initial_Mixed_Cell
10: ( n : in natural; mix : in Vector; pts : in Array_of_Lists;
11: mic : out Mixed_Cell; rest : in out Array_of_Lists ) is
12:
13: res : Mixed_Cell;
14: acc,tmp,rest_last : List;
15: pt : Link_to_Vector;
16: expts : Array_of_Lists(pts'range);
17: perm : Vector(pts'range);
18: mind,minl,dims,lent,index : integer;
19: nullvec : constant Vector := (1..n => 0);
20: shiftvecs : VecVec(pts'range);
21: grap : Matrix(1..n,pts'range);
22: fail : boolean := false;
23: rowcnt,cnt : natural := 0;
24: mat : Matrix(1..n,1..n+1);
25: ipvt : Vector(1..n+1);
26:
27: -- AUXILIAIRIES
28:
29: procedure Add ( pt,shift : in Vector; l : in out List ) is
30:
31: newpt : Link_to_Vector := new Vector'(pt);
32:
33: begin
34: Add(newpt.all,shift);
35: Construct(newpt,l);
36: end Add;
37:
38: procedure Add ( pt : in Vector; l : in out List ) is
39:
40: newpt : Link_to_Vector := new Vector'(pt);
41:
42: begin
43: Construct(newpt,l);
44: end Add;
45:
46: procedure Fill ( l : in List; rowcnt : in out natural; m : in out Matrix ) is
47: -- DESCRIPTION :
48: -- Fills the matrix with the elements in l.
49:
50: -- ON ENTRY :
51: -- l list of point vectors;
52: -- rowcnt m(m'first(1)..rowcnt) is already defined;
53: -- m matrix with m'range(2) = range of points in l.
54:
55: -- rowcnt updated row counter;
56: -- m matrix with the nonzero points of l
57:
58: tmp : List := l;
59: pt : Link_to_Vector;
60:
61: begin
62: while not Is_Null(tmp) loop
63: pt := Head_Of(tmp);
64: declare
65: nullvec : constant Vector(pt'range) := (pt'range => 0);
66: begin
67: if pt.all /= nullvec
68: then rowcnt := rowcnt + 1;
69: for k in pt'range loop
70: m(rowcnt,k) := pt(k);
71: end loop;
72: end if;
73: tmp := Tail_Of(tmp);
74: end;
75: end loop;
76: end Fill;
77:
78: procedure Initialize ( l : in List; rowcnt : out natural;
79: m : in out Matrix; ipvt : in out Vector ) is
80:
81: -- DESCRIPTION :
82: -- Given a list of linearly independent points, eventually with
83: -- 0 in l, a matrix will be constructed which contains these points
84: -- in upper triangular form.
85:
86: cnt : natural := 0;
87: piv : natural;
88:
89: begin
90: Fill(l,cnt,m);
91: for i in 1..cnt loop
92: m(i,m'last(2)) := i;
93: end loop;
94: for k in 1..cnt loop
95: Triangulate(k,m,ipvt,piv);
96: end loop;
97: rowcnt := cnt;
98: end Initialize;
99:
100: procedure Linearly_Independent
101: ( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
102: l : in List; len : out natural;
103: indep,indep_last : out List ) is
104:
105: -- DESCRIPTION :
106: -- Computes the list of points which are linearly independent w.r.t.
107: -- the matrix m.
108:
109: -- ON ENTRY :
110: -- m m(1..rowcnt,m'range(2)) is upper triangular,
111: -- can be used as work space;
112: -- rowcnt counter for the number of rows in m:
113: -- ipvt vector with the pivoting information;
114: -- l list of points to consider.
115:
116: -- ON RETURN :
117: -- len length of the list indep;
118: -- indep the list of linearly independent points;
119: -- indep_last pointer to the last element of indep.
120:
121: tmp,res,res_last : List;
122: pt : Link_to_Vector;
123: wipvt : Vector(ipvt'range) := ipvt;
124: piv,cnt : natural;
125:
126: begin
127: if rowcnt < m'last(2)
128: then tmp := l; cnt := 0;
129: while not Is_Null(tmp) loop
130: pt := Head_Of(tmp);
131: for i in pt'range loop
132: m(rowcnt+1,i) := pt(i);
133: end loop;
134: m(rowcnt+1,m'last(2)) := rowcnt+1;
135: Triangulate(rowcnt+1,m,wipvt,piv);
136: if m(rowcnt+1,rowcnt+1) /= 0
137: then cnt := cnt + 1;
138: Append(res,res_last,pt.all);
139: end if;
140: wipvt := ipvt;
141: tmp := Tail_Of(tmp);
142: end loop;
143: end if;
144: len := cnt;
145: indep := res; indep_last := res_last;
146: end Linearly_Independent;
147:
148: procedure Construct_Basis ( l : in List; rowcnt : in out natural;
149: m : in out Matrix; ipvt : in out Vector ) is
150:
151: -- DESCRIPTION :
152: -- Constructs a triangular basis for the vectors in l.
153:
154: -- REQUIRED : dimensions of the vectors must match.
155:
156: -- ON ENTRY :
157: -- l list of vector;
158: -- m triangular basis, stored in the rows of the matrix;
159: -- rowcnt index to the last significant row in m, could be 0,
160: -- when there is no initial basis to take into account;
161: -- ipvt initial pivoting vector, must be equal to the identity
162: -- permutation vector, when rowcnt = 0.
163:
164: -- ON RETURN :
165: -- m triangular basis, stored in the rows of the matrix;
166: -- rowcnt index to the last significant row in m;
167: -- ipvt contains pivoting information.
168:
169: tmp : List := l;
170: wipvt : Vector(ipvt'range);
171: pt : Link_to_Vector;
172: piv : natural;
173:
174: begin
175: while not Is_Null(tmp) loop
176: pt := Head_Of(tmp);
177: for i in pt'range loop
178: m(rowcnt+1,i) := pt(i);
179: end loop;
180: wipvt := ipvt;
181: m(rowcnt+1,m'last(2)) := rowcnt+1;
182: Triangulate(rowcnt+1,m,wipvt,piv);
183: if m(rowcnt+1,rowcnt+1) /= 0
184: then rowcnt := rowcnt + 1;
185: ipvt := wipvt;
186: end if;
187: exit when (rowcnt >= m'last(2));
188: tmp := Tail_Of(tmp);
189: end loop;
190: end Construct_Basis;
191:
192: function Linearly_Independent ( l : List; x : Vector ) return boolean is
193:
194: -- DESCRIPTION :
195: -- Returns true if the given vector x is linearly independent w.r.t.
196: -- the vectors in l.
197:
198: basis : Matrix(1..Length_Of(l)+1,x'first..x'last+1);
199: ipvt : Vector(basis'range(2));
200: rowcnt : natural := 0;
201: piv : natural;
202:
203: begin
204: for i in ipvt'range loop
205: ipvt(i) := i;
206: for j in basis'range(1) loop
207: basis(j,i) := 0;
208: end loop;
209: end loop;
210: Construct_Basis(l,rowcnt,basis,ipvt);
211: for i in x'range loop
212: basis(rowcnt+1,i) := x(i);
213: end loop;
214: basis(rowcnt+1,basis'last(2)) := rowcnt+1;
215: Triangulate(rowcnt+1,basis,ipvt,piv);
216: return (basis(rowcnt+1,rowcnt+1) /= 0);
217: end Linearly_Independent;
218:
219: procedure Linearly_Independent
220: ( l,xl : in List; indep,indep_last : in out List ) is
221:
222: -- DESCRIPTION :
223: -- Returns those points in xl that are linearly independent from the
224: -- points in l.
225:
226: begin
227: if not Is_Null(xl)
228: then
229: declare
230: x : Link_to_Vector := Head_Of(xl);
231: basis : Matrix(1..Length_Of(l)+1,x'first..x'last+1);
232: ipvt : Vector(basis'range(2));
233: rowcnt,len : natural := 0;
234: begin
235: for i in ipvt'range loop
236: ipvt(i) := i;
237: for j in basis'range(1) loop
238: basis(j,i) := 0;
239: end loop;
240: end loop;
241: Construct_Basis(l,rowcnt,basis,ipvt);
242: Linearly_Independent(basis,rowcnt,ipvt,xl,len,indep,indep_last);
243: end;
244: end if;
245: end Linearly_Independent;
246:
247: function Construct_Rest ( l : Array_of_Lists; start : natural )
248: return List is
249:
250: -- DESCRIPTION :
251: -- Collects all remaining points of l(i) into one single list,
252: -- for i in start..l'last.
253:
254: tmp,res,res_last : List;
255: pt : Link_to_Vector;
256:
257: begin
258: for i in start..l'last loop
259: tmp := l(i);
260: while not Is_Null(tmp) loop
261: pt := Head_Of(tmp);
262: if not Is_In(res,pt)
263: then Append(res,res_last,pt.all);
264: end if;
265: tmp := Tail_Of(tmp);
266: end loop;
267: end loop;
268: return res;
269: end Construct_Rest;
270:
271: procedure Sort_Co_Linearly_Independent
272: ( l : in out List; al : in Array_of_Lists;
273: start : in natural ) is
274:
275: -- DESCRIPTION :
276: -- Sorts the points in l, by putting those points that are linearly
277: -- independent w.r.t. the rest of the poins in al, in front of the list.
278:
279: rest : List := Construct_Rest(al,start);
280: tmp,indep,indep_last : List;
281: pt : Link_to_Vector;
282:
283: begin
284: -- put_line("The set of points in al just after Construct_Rest :"); put(al);
285: -- put_line("The rest of the points : "); put(rest);
286: Linearly_Independent(rest,l,indep,indep_last);
287: if Is_Null(indep)
288: then null; -- nothing to put in front of l
289: else tmp := l; -- append the other points of l to indep
290: -- put_line("Linearly co-independent points : "); put(indep);
291: -- put_line(" w.r.t. the rest : "); put(rest);
292: while not Is_Null(tmp) loop
293: pt := Head_Of(tmp);
294: if not Is_In(indep,pt)
295: then Append(indep,indep_last,pt.all);
296: end if;
297: tmp := Tail_Of(tmp);
298: end loop;
299: Copy(indep,l); Deep_Clear(indep);
300: end if;
301: Deep_Clear(rest);
302: -- put_line("The list of points in al after Sort_Co : "); put(al);
303: end Sort_Co_Linearly_Independent;
304:
305: procedure Incremental_Dimension
306: ( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
307: l : in List; dim : out natural ) is
308:
309: -- DESCRIPTION :
310: -- Computes the number of points which are linearly independent w.r.t.
311: -- the matrix m.
312:
313: -- ON ENTRY :
314: -- m m(1..rowcnt,m'range(2)) is upper triangular,
315: -- can be used as work space;
316: -- rowcnt counter for the number of rows in m:
317: -- ipvt vector with the pivoting information;
318: -- l list of points to consider.
319:
320: -- ON RETURN :
321: -- dim the number of linearly independent points.
322:
323: tmp : List;
324: pt : Link_to_Vector;
325: cnt,piv : natural;
326: newipvt,wipvt : Vector(ipvt'range);
327:
328: begin
329: wipvt := ipvt; newipvt := ipvt;
330: tmp := l; cnt := 0;
331: while not Is_Null(tmp) loop
332: pt := Head_Of(tmp);
333: cnt := cnt + 1;
334: for i in pt'range loop
335: m(rowcnt+cnt,i) := pt(i);
336: end loop;
337: m(rowcnt+cnt,m'last(2)) := rowcnt+cnt;
338: Triangulate(rowcnt+cnt,m,newipvt,piv);
339: if m(rowcnt+cnt,rowcnt+cnt) /= 0
340: then wipvt := newipvt;
341: else newipvt := wipvt;
342: cnt := cnt - 1;
343: end if;
344: tmp := Tail_Of(tmp);
345: exit when ((rowcnt + cnt) >= m'last(1));
346: end loop;
347: dim := cnt;
348: end Incremental_Dimension;
349:
350: procedure Incremental_Dimension
351: ( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
352: l : in out List; dim,len : out natural ) is
353:
354: -- DESCRIPTION :
355: -- Computes the increase in dimension by considering the points
356: -- in the list l, w.r.t. the matrix m.
357:
358: -- ON ENTRY :
359: -- m m(1..rowcnt,m'range(2)) is upper triangular,
360: -- rowcnt counter for the number of rows in m:
361: -- ipvt vector with the pivoting information;
362: -- l list of points to consider.
363:
364: -- ON RETURN :
365: -- l list of linearly independent points w.r.t. m;
366: -- dim increase in dimension;
367: -- len length of the list of linearly independent points
368: -- w.r.t. the rows in the matrix m.
369:
370: -- ALGORITHM :
371: -- works in two stages:
372: -- 1. determination of all linearly independent points;
373: -- 2. determination of the dimension.
374:
375: workm : Matrix(m'range(1),m'range(2));
376: indep,indep_last : List;
377:
378: begin
379: -- Initialize workm :
380: -- put_line("The list to investigate : "); put(l);
381: for i in 1..rowcnt loop
382: for j in m'range(2) loop
383: workm(i,j) := m(i,j);
384: end loop;
385: end loop;
386: -- Determine linearly independent points :
387: Linearly_Independent(workm,rowcnt,ipvt,l,len,indep,indep_last);
388: -- Determine the incremental dimension :
389: Incremental_Dimension(workm,rowcnt,ipvt,indep,dim);
390: Copy(indep,l); Deep_Clear(indep);
391: end Incremental_Dimension;
392:
393: procedure Next_Point ( acc,l : in List; pt : out Link_to_Vector;
394: fail : out boolean; rowcnt : in out natural;
395: m : in out Matrix; ipvt : in out Vector ) is
396:
397: -- DESCRIPTION :
398: -- A new point out of l, not in the list acc will be chosen.
399: -- The point pt has to be linearly independent from the other,
400: -- already chosen points. Therefore, the upper triangular matrix m
401: -- will be used and properly updated.
402:
403: res : Link_to_Vector;
404: tmp : List := l;
405: newipvt : Vector(ipvt'range) := ipvt;
406: done : boolean := false;
407: piv : natural;
408:
409: begin
410: while not Is_Null(tmp) loop
411: res := Head_Of(tmp);
412: if not Is_In(acc,res)
413: then --put("Checking point "); put(res); new_line;
414: rowcnt := rowcnt + 1;
415: for i in res'range loop
416: m(rowcnt,i) := res(i);
417: end loop;
418: m(rowcnt,m'last(2)) := rowcnt;
419: Triangulate(rowcnt,m,newipvt,piv);
420: if m(rowcnt,rowcnt) = 0
421: then rowcnt := rowcnt - 1;
422: newipvt := ipvt;
423: else ipvt := newipvt;
424: pt := res;
425: done := true;
426: end if;
427: end if;
428: exit when done;
429: tmp := Tail_Of(tmp);
430: end loop;
431: fail := not done;
432: end Next_Point;
433:
434: begin
435: -- INITIALIZE : compute for each list the basis points
436: res.nor := new Vector'(1..n+1 => 0); res.nor(n+1) := 1;
437: res.pts := new Array_of_Lists(mix'range);
438: for i in pts'range loop
439: expts(i) := Max_Extremal_Points(n,pts(i)); perm(i) := i;
440: -- put("Extremal points for component "); put(i,1); put_line(" :");
441: -- put(expts(i));
442: end loop;
443: -- INITIALIZE : order the lists according to occurencies
444: grap := Graph(n,expts);
445: -- put_line("The graph matrix of the extremal points : "); put(grap);
446: for i in expts'range loop
447: Sort(expts(i),grap);
448: -- put("Ordered extremal points for component "); put(i,1); put_line(" :");
449: -- put(expts(i));
450: end loop;
451: -- INITIALIZE : choose one anchor point for each component,
452: -- shift when necessary :
453: for i in expts'range loop
454: if Is_In(pts(i),nullvec)
455: then Add(nullvec,res.pts(i));
456: shiftvecs(i) := null;
457: else -- ADD LAST POINT = LEAST IMPORTANT
458: tmp := expts(i);
459: if not Is_Null(tmp)
460: then
461: while not Is_Null(Tail_Of(tmp)) loop
462: tmp := Tail_Of(tmp);
463: end loop;
464: pt := Head_Of(tmp);
465: Add(pt.all,res.pts(i));
466: shiftvecs(i) := new Vector'(pt.all);
467: -- put("The shift vector : "); put(shiftvecs(i)); new_line;
468: Shift(expts(i),shiftvecs(i));
469: -- put_line("The shifted list of points : "); put(expts(i));
470: end if;
471: end if;
472: end loop;
473: -- INITIALIZE : intermediate data structures mat and ipvt :
474: for i in ipvt'range loop
475: ipvt(i) := i;
476: end loop;
477: for i in mat'range(1) loop
478: for j in mat'range(2) loop
479: mat(i,j) := 0;
480: end loop;
481: end loop;
482: rowcnt := 0;
483: -- COMPUTE CELL : based on lists with extremal points
484: for i in expts'range loop -- MAIN LOOP
485: -- LOOK FOR SMALLEST SET :
486: index := i;
487: if i = 1
488: then mind := Length_Of(expts(i))-1;
489: for j in i+1..expts'last loop
490: dims := Length_Of(expts(j))-1;
491: if dims < mind
492: then index := j; mind := dims;
493: end if;
494: end loop;
495: else --put("Calling Incremental_Dimension for i = "); put(i,1); new_line;
496: Incremental_Dimension(mat,rowcnt,ipvt,expts(i),mind,minl);
497: --put(" dimension : "); put(mind,1);
498: --put(" length : " ); put(minl,1); new_line;
499: for j in i+1..expts'last loop
500: --put("Calling Incremental_Dimension for j = "); put(j,1); new_line;
501: Incremental_Dimension(mat,rowcnt,ipvt,expts(j),dims,lent);
502: --put(" dimension : "); put(dims,1);
503: --put(" length : " ); put(lent,1); new_line;
504: if (dims < mind) or else ((dims = mind) and then (lent < minl))
505: then index := j; mind := dims; minl := lent;
506: end if;
507: end loop;
508: end if;
509: -- put("The index : "); put(index,1); new_line;
510: -- put(" with minimal dimension : "); put(mind,1); new_line;
511: -- put("perm before permute : "); put(perm); new_line;
512: if index /= i
513: then tmp := expts(index); expts(index) := expts(i); expts(i) := tmp;
514: mind := perm(i); perm(i) := perm(index); perm(index) := mind;
515: end if;
516: -- SORT ACCORDING TO OCCURRENCIES AND TO CO-LINEARLY :
517: Sort(expts(i),grap,i-1,perm);
518: -- put_line("After first ordering : "); put(expts(i));
519: Sort_Co_Linearly_Independent(expts(i),expts,i+1);
520: -- put_line("The ordered list : "); put(expts(i));
521: -- put("perm after permute : "); put(perm); new_line;
522: -- CHOOSE THE POINTS :
523: if i = 1
524: then Add(nullvec,acc);
525: tmp := expts(i);
526: for j in 1..mix(perm(i)) loop
527: pt := Head_Of(tmp);
528: if pt.all = nullvec -- skip the null vector
529: then tmp := Tail_Of(tmp);
530: if not Is_Null(tmp)
531: then pt := Head_Of(tmp);
532: end if;
533: end if;
534: exit when Is_Null(tmp);
535: Add(pt.all,acc);
536: -- put_line("frequency graph before ignore :"); put(grap);
537: Ignore(grap,pt.all);
538: -- put_line("frequency graph after ignore :"); put(grap);
539: if shiftvecs(perm(i)) /= null
540: then Add(pt.all,shiftvecs(perm(i)).all,res.pts(perm(i)));
541: else Add(pt.all,res.pts(perm(i)));
542: end if;
543: tmp := Tail_Of(tmp);
544: exit when Is_Null(tmp);
545: end loop;
546: cnt := Length_Of(acc);
547: Initialize(acc,rowcnt,mat,ipvt);
548: -- put_line("The list acc : "); put(acc);
549: -- put_line("The matrix mat : "); put(mat,1,rowcnt);
550: -- put_line("The vector ipvt : "); put(ipvt); new_line;
551: fail := (cnt <= mix(perm(i)));
552: else for j in 1..mix(perm(i)) loop
553: Next_Point(acc,expts(i),pt,fail,rowcnt,mat,ipvt);
554: if not fail
555: then if shiftvecs(perm(i)) /= null
556: then Add(pt.all,shiftvecs(perm(i)).all,res.pts(perm(i)));
557: else Add(pt.all,res.pts(perm(i)));
558: end if;
559: Add(pt.all,acc); cnt := cnt + 1;
560: -- put_line("frequency graph before ignore :"); put(grap);
561: Ignore(grap,pt.all);
562: -- put_line("frequency graph after ignore :"); put(grap);
563: end if;
564: exit when fail;
565: end loop;
566: -- put_line("The list acc : "); put(acc);
567: -- put_line("The matrix mat : "); put(mat,1,rowcnt);
568: -- put_line("The vector ipvt : "); put(ipvt); new_line;
569: fail := (Length_Of(res.pts(perm(i))) <= mix(perm(i)));
570: end if;
571: exit when fail;
572: end loop; -- END OF MAIN LOOP
573: Deep_Clear(acc);
574: -- COMPUTE THE REST OF THE POINT LISTS :
575: if not fail
576: then for i in pts'range loop
577: tmp := pts(i);
578: rest_last := rest(i);
579: while not Is_Null(tmp) loop
580: pt := Head_Of(tmp);
581: if not Is_In(res.pts(i),pt)
582: then Append(rest(i),rest_last,pt);
583: end if;
584: tmp := Tail_Of(tmp);
585: end loop;
586: end loop;
587: end if;
588: -- GIVE THE POINTS IN THE INITIAL SIMPLEX LIFTING VALUE 0 :
589: for i in res.pts'range loop
590: tmp := res.pts(i);
591: while not Is_Null(tmp) loop
592: pt := Head_Of(tmp);
593: declare
594: lpt : Link_to_Vector;
595: begin
596: lpt := new Vector(1..n+1);
597: lpt(pt'range) := pt.all; lpt(n+1) := 0;
598: Set_Head(tmp,lpt);
599: end;
600: tmp := Tail_Of(tmp);
601: end loop;
602: end loop;
603: mic := res;
604: end Initial_Mixed_Cell;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to initial_mixed_cell.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift