Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/floating_linear_inequalities.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
2: with Givens_Rotations; use Givens_Rotations;
3:
4: package body Floating_Linear_Inequalities is
5:
6: -- USAGE OF LP :
7:
8: procedure Is_In_Cone ( tab : in out Matrix; lastcol : in integer;
9: rhs : in out Standard_Floating_Vectors.Vector;
10: tol : in double_float;
11: solution : out Standard_Floating_Vectors.Vector;
12: columns : out Standard_Integer_Vectors.Vector;
13: feasible : out boolean ) is
14:
15: -- ALGORITHM:
16: -- The following linear program will be solved:
17: --
18: -- min u_1 + .. + u_n
19: --
20: -- l_1*p_1i + l_2*p_2i + .. + l_m*p_mi + u_i*q_i = q_i i=1,2,..,n
21: --
22: -- to determine whether q belongs to the cone spanned by the
23: -- vectors p_1,p_2,..,p_m.
24: -- If all u_i are zero and all constraints are satisfied,
25: -- then q belongs to the cone spanned by the vectors.
26:
27: -- CONSTANTS :
28:
29: n : constant natural := rhs'last;
30: m : constant natural := 2*n; -- number of constraints
31: nb : constant natural := lastcol+n; -- number of unknowns
32:
33: -- VARIABLES :
34:
35: mat : matrix(0..m,0..nb);
36: sol : Standard_Floating_Vectors.Vector(1..nb) := (1..nb => 0.0);
37: bas : Standard_Integer_Vectors.Vector(1..m) := (1..m => 0);
38: slk : Standard_Integer_Vectors.Vector(1..nb) := (1..nb => 0);
39:
40: cnt : natural;
41: s : double_float;
42:
43: begin
44:
45: -- INITIALIZATION OF THE TARGET :
46:
47: for i in 0..(nb-n) loop
48: mat(0,i) := 0.0; -- sum of the lambda's
49: end loop;
50: for i in (nb-n+1)..nb loop
51: mat(0,i) := -1.0; -- sum of the mu's
52: end loop;
53:
54: -- INITIALIZATION OF THE CONSTRAINTS :
55:
56: for i in 1..n loop
57: for j in (nb-n+1)..nb loop
58: if i /= j-nb+n
59: then mat(i,j) := 0.0;
60: mat(i+n,j) := 0.0;
61: end if;
62: end loop;
63: end loop;
64:
65: for j in tab'first(2)..lastcol loop
66: for i in tab'range(1) loop
67: mat(i,j) := -tab(i,j);
68: mat(i+n,j) := tab(i,j);
69: end loop;
70: end loop;
71:
72: for i in rhs'range loop
73: mat(i,0) := -rhs(i);
74: mat(i+n,0) := rhs(i);
75: mat(i,i+nb-n) := -rhs(i);
76: mat(i+n,i+nb-n) := rhs(i);
77: end loop;
78:
79: -- SOLVE THE LINEAR PROGRAM :
80:
81: -- New_Tableau.Init_Basis(nb,bas,slk);
82: -- Dual_Simplex(mat,sol,bas,slk,tol,false);
83:
84: -- RETURN AND CHECK THE SOLUTION :
85:
86: cnt := 0;
87: for i in tab'first(2)..lastcol loop
88: if abs(sol(i)) > tol
89: then cnt := cnt + 1;
90: solution(cnt) := sol(i);
91: columns(cnt) := i;
92: end if;
93: end loop;
94:
95: s := 0.0;
96: for i in (nb-n+1)..nb loop
97: s := s + sol(i);
98: end loop;
99: if abs(s) > tol
100: then feasible := false; return;
101: end if;
102: feasible := true;
103:
104: end Is_In_Cone;
105:
106: -- AUXILIARIES :
107:
108: function Index_of_Maximum ( v : Standard_Floating_Vectors.Vector;
109: lst : integer ) return integer is
110:
111: -- DESCRIPTION :
112: -- Returns the index in the vector v(v'first..lst) for which the maximum
113: -- is obtained.
114:
115: res : integer := v'first;
116: max : double_float := v(res);
117:
118: begin
119: for i in v'first+1..lst loop
120: if v(i) > max
121: then max := v(i); res := i;
122: end if;
123: end loop;
124: return res;
125: end Index_of_Maximum;
126:
127: function Norm ( v : Standard_Floating_Vectors.Vector ) return double_float is
128:
129: -- DESCRIPTION :
130: -- Returns the 2-norm of the vector.
131:
132: res : double_float := 0.0;
133:
134: begin
135: for j in v'range loop
136: res := res + v(j)*v(j);
137: end loop;
138: if res + 1.0 = 1.0
139: then return 0.0;
140: else return SQRT(res);
141: end if;
142: end Norm;
143:
144: function Norm ( v : Standard_Floating_Vectors.vector;
145: frst : integer ) return double_float is
146:
147: -- DESCRIPTION :
148: -- Returns the 2-norm of the vector v(frst..v'last).
149:
150: res : double_float := 0.0;
151:
152: begin
153: for j in frst..v'last loop
154: res := res + v(j)*v(j);
155: end loop;
156: if res + 1.0 = 1.0
157: then return 0.0;
158: else return SQRT(res);
159: end if;
160: end Norm;
161:
162: procedure Scale ( v : in out Standard_Floating_Vectors.vector ) is
163:
164: -- DESCRIPTION :
165: -- Returns the normed vector, i.e. v/Norm(v).
166:
167: nrm : double_float := Norm(v);
168:
169: begin
170: if nrm + 1.0 /= 1.0
171: then for i in v'range loop
172: v(i) := v(i)/nrm;
173: end loop;
174: end if;
175: end Scale;
176:
177: function Norm ( mat : matrix; i : integer ) return double_float is
178:
179: -- DESCRIPTION :
180: -- Returns the 2-norm of the ith column in the matrix.
181:
182: res : double_float := 0.0;
183:
184: begin
185: for j in mat'range(1) loop
186: res := res + mat(j,i)*mat(j,i);
187: end loop;
188: if res + 1.0 = 1.0
189: then return 0.0;
190: else return SQRT(res);
191: end if;
192: end Norm;
193:
194: procedure Scale ( mat : in out matrix; lastcolumn : in integer ) is
195:
196: -- DESCRIPTION :
197: -- Scales the columns in the matrix, by dividing by their norm.
198:
199: nrm : double_float;
200:
201: begin
202: for i in mat'first(2)..lastcolumn loop
203: nrm := Norm(mat,i);
204: if nrm + 1.0 /= 1.0
205: then for j in mat'range(1) loop
206: mat(j,i) := mat(j,i)/nrm;
207: end loop;
208: end if;
209: end loop;
210: end Scale;
211:
212: function Inner_Product ( mat : Matrix; i : integer;
213: v : Standard_Floating_Vectors.Vector )
214: return double_float is
215: -- DESCRIPTION :
216: -- Computes the inner product of the given vector v and the vector
217: -- in the ith column of the matrix.
218:
219: res : double_float := 0.0;
220:
221: begin
222: for k in mat'range(1) loop
223: res := res + mat(k,i)*v(k);
224: end loop;
225: return res;
226: end Inner_Product;
227:
228: function Maximal_Product ( mat : matrix; i : integer;
229: v : Standard_Floating_Vectors.Vector;
230: frst : integer ) return integer is
231: -- DESCRIPTION :
232: -- Returns the index in v(frst..v'last) for which mat(index,i)*v(index)
233: -- becomes maximal.
234:
235: res : integer := frst;
236: max : double_float := mat(res,i)*v(res);
237: tmp : double_float;
238:
239: begin
240: for j in frst+1..v'last loop
241: tmp := mat(j,i)*v(j);
242: if tmp > max
243: then max := tmp; res := j;
244: end if;
245: end loop;
246: return res;
247: end Maximal_Product;
248:
249: function Inner_Product ( mat : Matrix; i,j : integer ) return double_float is
250:
251: -- DESCRIPTION :
252: -- Computes the inner product of the vectors in the columns i and j
253: -- of the matrix.
254:
255: res : double_float := 0.0;
256:
257: begin
258: for k in mat'range(1) loop
259: res := res + mat(k,i)*mat(k,j);
260: end loop;
261: return res;
262: end Inner_Product;
263:
264: function Inner_Products ( mat : matrix; lastcol : integer;
265: v : Standard_Floating_Vectors.Vector )
266: return Standard_Floating_Vectors.Vector is
267:
268: -- DESCRIPTION :
269: -- Returns a vector with all inner products of the given vector
270: -- and the column vectors of the matrix.
271:
272: res : Standard_Floating_Vectors.Vector(mat'first(2)..lastcol);
273:
274: begin
275: for i in res'range loop
276: res(i) := Inner_Product(mat,i,v);
277: end loop;
278: return res;
279: end Inner_Products;
280:
281: function Positive ( v : Standard_Floating_Vectors.Vector;
282: lst : integer; tol : double_float )
283: return boolean is
284: -- DESCRIPTION :
285: -- Returns true if all elements in v(v'first..lst) are >= 0.
286:
287: begin
288: for i in v'first..lst loop
289: if abs(v(i)) > tol and then (v(i) < 0.0)
290: then return false;
291: end if;
292: end loop;
293: return true;
294: end Positive;
295:
296: function Positive
297: ( mat : matrix; rhs : Standard_Floating_Vectors.Vector;
298: ind,col : integer;
299: solind : double_float; tol : double_float ) return boolean is
300:
301: -- DESCRIPTION :
302: -- The given matrix is upper triangular up to the column indicated
303: -- by ind. This function returns true when by replacing column ind
304: -- by column col, a positive solution would be obtained.
305: -- The number solind is sol(ind) in the solution vector.
306:
307: sol : Standard_Floating_Vectors.Vector(1..ind) := (1..ind => 0.0);
308:
309: begin
310: sol(ind) := solind;
311: for i in reverse mat'first(1)..(ind-1) loop
312: for j in i+1..(ind-1) loop
313: sol(i) := sol(i) + mat(i,j)*sol(j);
314: end loop;
315: sol(i) := sol(i) + mat(i,col)*sol(ind);
316: sol(i) := (rhs(i) - sol(i))/mat(i,i);
317: if (abs(sol(i)) > tol) and then (sol(i) < 0.0)
318: then --put_line("The solution before returning false : ");
319: --put(sol,3,3,3); new_line;
320: return false;
321: end if;
322: end loop;
323: --put_line("The solution before returning true : ");
324: --put(sol,3,3,3); new_line;
325: return true;
326: end Positive;
327:
328: function Positive_Diagonal
329: ( mat : Matrix; rhs : Standard_Floating_Vectors.Vector;
330: ind,col : integer;
331: solind : double_float; tol : double_float ) return boolean is
332:
333: -- DESCRIPTION :
334: -- The given matrix is diagonal up to the column indicated
335: -- by ind. This function returns true when by replacing column ind
336: -- by column col, a positive solution would be obtained.
337: -- The number solind is sol(ind) in the solution vector.
338:
339: sol : Standard_Floating_Vectors.Vector(1..ind) := (1..ind => 0.0);
340:
341: begin
342: sol(ind) := solind;
343: for i in reverse mat'first(1)..(ind-1) loop
344: sol(i) := (rhs(i) - mat(i,col)*sol(ind))/mat(i,i);
345: if (abs(sol(i)) > tol) and then (sol(i) < 0.0)
346: then --put_line("The solution before returning false : ");
347: --put(sol,3,3,3); new_line;
348: return false;
349: end if;
350: end loop;
351: --put_line("The solution before returning true : ");
352: --put(sol,3,3,3); new_line;
353: return true;
354: end Positive_Diagonal;
355:
356: -- PIVOTING PROCEDURES :
357:
358: procedure Interchange ( v : in out Standard_Floating_Vectors.Vector;
359: i,j : in integer ) is
360:
361: -- DESCRIPTION :
362: -- The entries v(i) and v(j) are interchanged.
363:
364: temp : double_float;
365:
366: begin
367: temp := v(i); v(i) := v(j); v(j) := temp;
368: end Interchange;
369:
370: procedure Interchange_Columns ( mat : in out matrix; i,j : in integer ) is
371:
372: -- DESCRIPTION :
373: -- The columns i and j of the given matrix are interchanged.
374:
375: temp : double_float;
376:
377: begin
378: for k in mat'range(1) loop
379: temp := mat(k,i); mat(k,i) := mat(k,j); mat(k,j) := temp;
380: end loop;
381: end Interchange_Columns;
382:
383: procedure Interchange_Rows
384: ( mat : in out Matrix; lastcol : in integer;
385: v : in out Standard_Floating_Vectors.Vector;
386: i,j : in integer ) is
387:
388: -- DESCRIPTION :
389: -- The rows i and j of the given matrix and vector are interchanged.
390:
391: temp : double_float;
392:
393: begin
394: for k in mat'first(2)..lastcol loop
395: temp := mat(i,k); mat(i,k) := mat(j,k); mat(j,k) := temp;
396: end loop;
397: temp := v(i); v(i) := v(j); v(j) := temp;
398: end Interchange_Rows;
399:
400: procedure Rotate ( mat : in out Matrix; lastcol,i : in integer;
401: v : in out Standard_Floating_Vectors.Vector;
402: tol : in double_float ) is
403:
404: -- DESCRIPTION :
405: -- Performs Givens rotations on the matrix and vector such that
406: -- mat(j,i) = 0, for all j > i.
407:
408: -- NOTE :
409: -- A diagonalization procedure constructs an orthonormal basis
410: -- and does not preserve the angles between the vectors.
411:
412: cosi,sine : double_float;
413:
414: begin
415: for j in i+1..mat'last(1) loop
416: if abs(mat(j,i)) > tol
417: then Givens_Factors(mat,i,j,cosi,sine);
418: Givens_Rotation(mat,lastcol,i,j,cosi,sine);
419: Givens_Rotation(v,i,j,cosi,sine);
420: end if;
421: end loop;
422: end Rotate;
423:
424: procedure Upper_Triangulate ( mat : in out Matrix; lastcol,i : in integer;
425: v : in out Standard_Floating_Vectors.Vector;
426: tol : in double_float ) is
427:
428: -- DESCRIPTION :
429: -- Makes the upper diagonal elements of the ith colume zero:
430: -- mat(j,i) = 0, for all j /= i.
431:
432: fac : double_float;
433:
434: begin
435: for j in mat'first(1)..(i-1) loop
436: if (abs(mat(j,i)) > tol)
437: then fac := mat(j,i)/mat(i,i);
438: for k in i..lastcol loop
439: mat(j,k) := mat(j,k) - fac*mat(i,k);
440: end loop;
441: end if;
442: end loop;
443: end Upper_Triangulate;
444:
445: -- INITIALIZE : COMPUTE CLOSEST VECTOR TO RHS
446:
447: procedure Initialize
448: ( tableau : in out Matrix; lastcol : in integer;
449: rhs,sol : in out Standard_Floating_Vectors.Vector;
450: tol : in double_float;
451: columns : in out Standard_Integer_Vectors.Vector;
452: cosines : out Standard_Floating_Vectors.Vector;
453: feasible : out boolean ) is
454:
455: -- DESCRIPTION :
456: -- Scales the tableau and right hand side by dividing by its norm.
457: -- Then computes the vector of cosines.
458: -- If no vector with positive cosine with the right hand side vector
459: -- can be found, then the problem is not feasible, otherwise,
460: -- the vector with largest positive cosine will be the first column.
461: -- At last, this first column will be transformed to the unit vector,
462: -- by means of Givens rotations.
463: -- Then the first solution component equals the first component of
464: -- the transformed right hand side vector.
465:
466: ips : Standard_Floating_Vectors.Vector(tableau'first(2)..lastcol);
467: index : integer;
468:
469: begin
470: Scale(tableau,lastcol); Scale(rhs);
471: -- put_line("The tableau after scaling : "); put(tableau,3,3,3);
472: -- put_line(" with scaled right hand side : "); put(rhs,3,3,3); new_line;
473: ips := Inner_Products(tableau,lastcol,rhs);
474: index := Index_of_Maximum(ips,lastcol);
475: if (abs(ips(index)) < tol) or else (ips(index) < 0.0)
476: then if index <= rhs'last
477: then feasible := (abs(rhs(index)) < tol);
478: else feasible := false;
479: end if;
480: else feasible := true;
481: columns(columns'first) := index;
482: columns(index) := columns'first;
483: Interchange_Columns(tableau,tableau'first(2),index);
484: Interchange(ips,tableau'first(2),index);
485: index := Maximal_Product(tableau,tableau'first(2),rhs,rhs'first);
486: if index /= tableau'first(2)
487: then Interchange_Rows(tableau,lastcol,rhs,tableau'first(2),index);
488: end if;
489: Rotate(tableau,lastcol,tableau'first(2),rhs,tol);
490: -- Upper_Triangulate(tableau,lastcol,tableau'first(2),rhs,tol);
491: sol(sol'first) := rhs(rhs'first);
492: end if;
493: cosines := ips;
494: -- put_line("At the end of Initialize : ");
495: -- put(" cosines : "); put(ips,3,3,3); new_line;
496: -- put_line("The tableau : "); put(tableau,3,3,3);
497: -- put_line(" with right hand side : "); put(rhs,3,3,3); new_line;
498: -- put(" and last column "); put(lastcol,1); new_line;
499: end Initialize;
500:
501: procedure Initialize
502: ( tableau : in out Matrix; lastcol : in integer;
503: rhs,sol : in out Standard_Floating_Vectors.Vector;
504: tol : in double_float;
505: columns : in out Standard_Integer_Vectors.Vector;
506: feasible : out boolean ) is
507:
508: -- DESCRIPTION :
509: -- Scales the tableau and right hand side by dividing by its norm.
510: -- Then computes the vector of cosines.
511: -- If no vector with positive cosine with the right hand side vector
512: -- can be found, then the problem is not feasible, otherwise,
513: -- the vector with largest positive cosine will be the first column.
514: -- At last, this first column will be transformed to the unit vector,
515: -- by means of Givens rotations.
516: -- Then the first solution component equals the first component of
517: -- the transformed right hand side vector.
518:
519: ips : Standard_Floating_Vectors.Vector(tableau'first(2)..lastcol);
520: index : integer;
521:
522: begin
523: Scale(tableau,lastcol); Scale(rhs);
524: -- put_line("The tableau after scaling : "); put(tableau,3,3,3);
525: -- put_line(" with scaled right hand side : "); put(rhs,3,3,3); new_line;
526: ips := Inner_Products(tableau,lastcol,rhs);
527: index := Index_of_Maximum(ips,lastcol);
528: if (abs(ips(index)) < tol) or else (ips(index) < 0.0)
529: then if index <= rhs'last
530: then feasible := (abs(rhs(index)) < tol);
531: else feasible := false;
532: end if;
533: else feasible := true;
534: columns(columns'first) := index;
535: columns(index) := columns'first;
536: Interchange_Columns(tableau,tableau'first(2),index);
537: index := Maximal_Product(tableau,tableau'first(2),rhs,rhs'first);
538: if index /= tableau'first(2)
539: then Interchange_Rows(tableau,lastcol,rhs,tableau'first(2),index);
540: end if;
541: Rotate(tableau,lastcol,tableau'first(2),rhs,tol);
542: -- Upper_Triangulate(tableau,lastcol,tableau'first(2),rhs,tol);
543: sol(sol'first) := rhs(rhs'first);
544: end if;
545: -- put_line("At the end of Initialize : ");
546: -- put_line("The tableau : "); put(tableau,3,3,3);
547: -- put_line(" with right hand side : "); put(rhs,3,3,3); new_line;
548: -- put(" and last column "); put(lastcol,1); new_line;
549: end Initialize;
550:
551: -- PERFORM THE NEXT STEPS :
552:
553: procedure Select_Column
554: ( mat : in Matrix; lastcol : in integer;
555: rhs,cosines : in Standard_Floating_Vectors.Vector;
556: index : in integer;
557: tol : in double_float; row,col : out integer ) is
558: -- DESCRIPTION :
559: -- Selects a column col with in matrix for which
560: -- 1) the cosine is maximal and
561: -- 2) row = Maximal_Product(mat,col,rhs,index), with mat(row,col) > 0.0
562: -- If the resulting column is smaller than index, then no such column
563: -- has been found.
564:
565: cl : integer := index;
566: maxcos,prod : double_float;
567: rw : integer := Maximal_Product(mat,cl,rhs,index);
568: mp : integer;
569:
570: begin
571: -- put_line("The tableau when Select_Column :"); put(mat,3,3,3);
572: -- put_line(" with right hand side : "); put(rhs,3,3,3); new_line;
573: -- put_line("and vector of cosines : "); put(cosines,3,3,3); new_line;
574: -- put("The index : "); put(index,1); new_line;
575: prod := mat(rw,cl)*rhs(rw);
576: if prod > tol
577: and then Positive(mat,rhs,index,cl,rhs(rw)/mat(rw,cl),tol)
578: -- Positive_Diagonal(mat,rhs,index,cl,rhs(rw)/mat(rw,cl),tol)
579: then maxcos := cosines(cl);
580: else maxcos := -2.0;
581: end if;
582: for i in index+1..lastcol loop
583: mp := Maximal_Product(mat,i,rhs,index);
584: -- put("mp : "); put(mp,1); put(" when testing "); put(i,1);
585: -- put_line("th column");
586: prod := mat(mp,i)*rhs(mp);
587: if prod > tol and then cosines(i) > maxcos
588: then if Positive(mat,rhs,index,i,rhs(mp)/mat(mp,i),tol)
589: -- Positive_Diagonal(mat,rhs,index,i,rhs(mp)/mat(mp,i),tol)
590: then cl := i; rw := mp; maxcos := cosines(i);
591: end if;
592: end if;
593: end loop;
594: -- put("maximal cosine : "); put(maxcos,3,3,3); new_line;
595: if maxcos >= -1.0
596: then col := cl; row := rw;
597: -- put("column : "); put(cl,1);
598: -- put(" and row : "); put(rw,1); new_line;
599: else col := index-1; row := index-1;
600: -- put("column : "); put(index-1,1);
601: -- put(" and row : "); put(index-1,1); new_line;
602: end if;
603: end Select_Column;
604:
605: procedure Select_Column
606: ( mat : in Matrix; lastcol : in integer;
607: rhs : in Standard_Floating_Vectors.Vector;
608: index,column : in integer;
609: tol : in double_float; row,col : out integer ) is
610: -- DESCRIPTION :
611: -- Selects first column col with in matrix for which
612: -- row = Maximal_Product(mat,col,rhs,index), with mat(row,col) > 0.0.
613: -- If the resulting column is smaller than index, then no such column
614: -- has been found.
615:
616: cl : integer := index-1;
617: prod : double_float;
618: rw,mp : integer;
619:
620: begin
621: -- put_line("The tableau when Select_Column :"); put(mat,3,3,3);
622: -- put_line(" with right hand side : "); put(rhs,3,3,3); new_line;
623: -- put("The index : "); put(index,1); new_line;
624: for i in column..lastcol loop
625: mp := Maximal_Product(mat,i,rhs,index);
626: -- put("mp : "); put(mp,1); put(" when testing "); put(i,1);
627: -- put_line("th column");
628: prod := mat(mp,i)*rhs(mp);
629: if prod > tol
630: and then Positive(mat,rhs,index,i,rhs(mp)/mat(mp,i),tol)
631: -- Positive_Diagonal(mat,rhs,index,i,rhs(mp)/mat(mp,i),tol)
632: then cl := i; rw := mp;
633: end if;
634: exit when (cl >= index);
635: end loop;
636: if cl >= index
637: then col := cl; row := rw;
638: -- put("column : "); put(cl,1);
639: -- put(" and row : "); put(rw,1); new_line;
640: else col := index-1; row := index-1;
641: -- put("column : "); put(index-1,1);
642: -- put(" and row : "); put(index-1,1); new_line;
643: end if;
644: end Select_Column;
645:
646: procedure Next ( tableau : in out Matrix; lastcol : in integer;
647: rhs,sol,cosines : in out Standard_Floating_Vectors.Vector;
648: pivots : in out Standard_Integer_Vectors.Vector;
649: tol : in double_float;
650: index : in integer; feasi : in out boolean ) is
651:
652: -- DESCRIPTION :
653: -- Selects one new column out of the tableau.
654:
655: -- ON ENTRY :
656: -- tableau upper triangular up to the row and column index-1;
657: -- lastcol index of the last column of interest in the tableau;
658: -- rhs right hand side vector;
659: -- sol solution vector for the components 1..index;
660: -- cosines vector of cosines;
661: -- pivots vector of column interchangements;
662: -- tol tolerance to decide whether a number equals zero;
663: -- index indicator of the current step;
664: -- feasi must be true on entry.
665:
666: -- ON RETURN :
667: -- tableau upper triangular up to the row and column index,
668: -- under the condition that feasi remains true;
669: -- rhs transformed right hand side vector;
670: -- sol new solution, with additional meaningful component
671: -- sol(index);
672: -- cosines eventually permuted vector of cosines:
673: -- cosines(index) is the cosine of the (new) column,
674: -- indicated by index, and the right hand side vector;
675: -- pivots updated vector of column interchangements;
676: -- feasi if true, then a new column which yields a positive
677: -- contribution of the right hand side vector has been
678: -- found, if this was not the case, then feasi is false;
679:
680: col,row,tmp : integer;
681:
682: begin
683: Select_Column(tableau,lastcol,rhs,cosines,index,tol,row,col);
684: if col < index
685: then feasi := false;
686: else feasi := true;
687: if col /= index
688: then tmp := pivots(col); pivots(col) := pivots(index);
689: pivots(index) := tmp;
690: Interchange_Columns(tableau,col,index);
691: Interchange(cosines,col,index);
692: end if;
693: if row /= index
694: then Interchange_Rows(tableau,lastcol,rhs,row,index);
695: end if;
696: Rotate(tableau,lastcol,index,rhs,tol);
697: -- Upper_Triangulate(tableau,lastcol,index,rhs,tol);
698: -- Solve(tableau,rhs,tol,index,sol); -- only needed at end
699: -- feasi := Positive(sol,index,tol); -- double check...
700: end if;
701: end Next;
702:
703: procedure Next ( tableau : in out Matrix; lastcol : in integer;
704: rhs,sol : in out Standard_Floating_Vectors.Vector;
705: pivots : in out Standard_Integer_Vectors.Vector;
706: tol : in double_float;
707: index : in integer; rank : out integer;
708: feasi : in out boolean ) is
709:
710: -- DESCRIPTION :
711: -- Lexicographic enumeration of the candidates, stops when a feasible
712: -- solution is encountered.
713:
714: -- ON ENTRY :
715: -- tableau upper triangular up to the row and column index-1;
716: -- lastcol index of the last column of interest in the tableau;
717: -- rhs right hand side vector;
718: -- sol solution vector for the components 1..index;
719: -- pivots vector of column interchangements;
720: -- tol tolerance to decide whether a number equals zero;
721: -- index indicator of the current step;
722: -- feasi must be true on entry.
723:
724: -- ON RETURN :
725: -- tableau upper triangular up to the row and column index,
726: -- under the condition that feasi remains true;
727: -- rhs transformed right hand side vector;
728: -- sol new solution, with additional meaningful component
729: -- sol(index);
730: -- pivots updated vector of column interchangements;
731: -- rank rank of the current tableau;
732: -- feasi if true, then a new column which yields a positive
733: -- contribution of the right hand side vector has been
734: -- found, if this was not the case, then feasi is false;
735:
736: column,col,row,tmp : integer;
737: stop : boolean := false;
738:
739: begin
740: column := index;
741: loop
742: Select_Column(tableau,lastcol,rhs,index,column,tol,row,col);
743: if col < column
744: then feasi := false; stop := true;
745: else
746: feasi := true;
747: if col /= index
748: then tmp := pivots(col); pivots(col) := pivots(index);
749: pivots(index) := tmp;
750: Interchange_Columns(tableau,col,index);
751: end if;
752: if row /= index
753: then Interchange_Rows(tableau,lastcol,rhs,row,index);
754: end if;
755: Rotate(tableau,lastcol,index,rhs,tol);
756: -- Upper_Triangulate(tableau,lastcol,index,rhs,tol);
757: if (index = lastcol) or else (index = tableau'last(1))
758: or else (Norm(rhs,index+1) < tol)
759: then stop := true; rank := index;
760: else
761: Next(tableau,lastcol,rhs,sol,pivots,tol,index+1,rank,feasi);
762: if feasi
763: then stop := true;
764: else column := col+1;
765: end if;
766: end if;
767: end if;
768: exit when stop;
769: end loop;
770: end Next;
771:
772: procedure Complementary_Slackness
773: ( tableau : in out matrix;
774: rhs : in out Standard_Floating_Vectors.Vector;
775: tol : in double_float;
776: solution : out Standard_Floating_Vectors.Vector;
777: columns : out Standard_Integer_Vectors.Vector;
778: feasible : out boolean ) is
779: begin
780: Complementary_Slackness
781: (tableau,tableau'last(2),rhs,tol,solution,columns,feasible);
782: end Complementary_Slackness;
783:
784: procedure Complementary_Slackness1
785: ( tableau : in out Matrix; lastcol : in integer;
786: rhs : in out Standard_Floating_Vectors.Vector;
787: tol : in double_float;
788: solution : out Standard_Floating_Vectors.Vector;
789: columns : out Standard_Integer_Vectors.Vector;
790: feasible : out boolean ) is
791:
792: -- NOTE :
793: -- This version is a straigth forward search and finds not always
794: -- the feasible solution.
795:
796: feasi : boolean;
797: pivots : Standard_Integer_Vectors.Vector(tableau'first(2)..lastcol);
798: cosis : Standard_Floating_Vectors.Vector(tableau'first(2)..lastcol);
799: sol : Standard_Floating_Vectors.Vector(rhs'range) := (rhs'range => 0.0);
800: rank : integer;
801:
802: begin
803: for i in pivots'range loop -- initialize vector of pivots
804: pivots(i) := i;
805: end loop;
806: Initialize(tableau,lastcol,rhs,sol,tol,pivots,cosis,feasi);
807: if feasi
808: then
809: if Norm(rhs,rhs'first+1) > tol -- check whether residual > 0
810: then
811: if tableau'first(1)+1 > lastcol
812: then feasi := false;
813: else
814: rank := 1;
815: for i in tableau'first(1)+1..tableau'last(1) loop
816: Next(tableau,lastcol,rhs,sol,cosis,pivots,tol,i,feasi);
817: exit when not feasi;
818: exit when (i = lastcol);
819: rank := rank + 1;
820: exit when (Norm(rhs,i+1) < tol);
821: end loop;
822: if feasi
823: then Solve(tableau,rhs,tol,rank,sol);
824: end if;
825: end if;
826: end if;
827: if feasi and then (tableau'last(1) > lastcol)
828: then feasi := (Norm(rhs,lastcol+1) < tol);
829: end if;
830: if feasi and then (Norm(sol) < tol)
831: then feasi := (Norm(rhs) < tol);
832: end if;
833: solution := sol;
834: if tableau'last(1) <= lastcol
835: then
836: for i in tableau'range(1) loop
837: columns(i) := pivots(i);
838: end loop;
839: else
840: for i in tableau'first(2)..lastcol loop
841: columns(i) := pivots(i);
842: end loop;
843: end if;
844: end if;
845: feasible := feasi;
846: end Complementary_Slackness1;
847:
848: procedure Complementary_Slackness
849: ( tableau : in out Matrix; lastcol : in integer;
850: rhs : in out Standard_Floating_Vectors.Vector;
851: tol : in double_float;
852: solution : out Standard_Floating_Vectors.Vector;
853: columns : out Standard_Integer_Vectors.Vector;
854: feasible : out boolean ) is
855:
856: -- NOTE :
857: -- This version enumerates in a lexicographic order all candidates
858: -- for entering the basis and stops when a feasible solution has been
859: -- found.
860:
861: feasi : boolean;
862: pivots : Standard_Integer_Vectors.Vector(tableau'first(2)..lastcol);
863: sol : Standard_Floating_Vectors.Vector(rhs'range) := (rhs'range => 0.0);
864: ind,rank : integer;
865:
866: begin
867: for i in pivots'range loop -- initialize vector of pivots
868: pivots(i) := i;
869: end loop;
870: Initialize(tableau,lastcol,rhs,sol,tol,pivots,feasi);
871: if feasi
872: then
873: if Norm(rhs,rhs'first+1) > tol -- check whether residual > 0
874: then
875: if tableau'first(1)+1 > lastcol
876: then feasi := false;
877: else
878: ind := tableau'first(1)+1;
879: Next(tableau,lastcol,rhs,sol,pivots,tol,ind,rank,feasi);
880: if feasi
881: then Solve(tableau,rhs,tol,rank,sol);
882: end if;
883: end if;
884: end if;
885: if feasi and then (tableau'last(1) > lastcol)
886: then feasi := (Norm(rhs,lastcol+1) < tol);
887: end if;
888: if feasi and then (Norm(sol) < tol)
889: then feasi := (Norm(rhs) < tol);
890: end if;
891: solution := sol;
892: if tableau'last(1) <= lastcol
893: then
894: for i in tableau'range(1) loop
895: columns(i) := pivots(i);
896: end loop;
897: else
898: for i in tableau'first(2)..lastcol loop
899: columns(i) := pivots(i);
900: end loop;
901: end if;
902: end if;
903: feasible := feasi;
904: end Complementary_Slackness;
905:
906: procedure Complementary_Slackness2
907: ( tableau : in out Matrix; lastcol : in integer;
908: rhs : in out Standard_Floating_Vectors.Vector;
909: tol : in double_float;
910: solution : out Standard_Floating_Vectors.Vector;
911: columns : out Standard_Integer_Vectors.Vector;
912: feasible : out boolean ) is
913:
914: -- NOTE :
915: -- This version explicitly uses linear programming.
916:
917: begin
918: Is_In_Cone(tableau,lastcol,rhs,tol,solution,columns,feasible);
919: end Complementary_Slackness2;
920:
921: end Floating_Linear_Inequalities;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to floating_linear_inequalities.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports