![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic_integer_linear_solvers.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 / Matrices|
Return to generic_integer_linear_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 maekawa 1: package body Generic_Integer_Linear_Solvers is
2:
3: use Common_Divisors;
4:
5: -- SCALERS :
6:
7: function Divisors ( a : Matrix ) return Vectors.Vector is
8:
9: -- DESCRIPTION :
10: -- Returns a vector containing the gcd of the elements of each row.
11:
12: v : Vectors.Vector(a'range(1));
13:
14: begin
15: for i in a'range(1) loop
16: v(i) := a(i,a'first(2));
17: if not Equal(v(i),one)
18: then for j in (a'first(2)+1)..a'last(2) loop
19: v(i) := gcd(v(i),a(i,j));
20: exit when Equal(v(i),one);
21: end loop;
22: end if;
23: end loop;
24: return v;
25: end Divisors;
26:
27: function Scale ( a : Matrix ) return Matrix is
28:
29: v : Vectors.Vector(a'range(1)) := Divisors(a);
30: b : Matrix(a'range(1),a'range(2));
31:
32: begin
33: for i in b'range(1) loop
34: if (Equal(v(i),one) or Equal(v(i),zero))
35: then for j in b'range(2) loop
36: b(i,j) := a(i,j);
37: end loop;
38: else for j in b'range(2) loop
39: b(i,j) := a(i,j) / v(i);
40: end loop;
41: end if;
42: end loop;
43: return b;
44: end Scale;
45:
46: procedure Scale ( a : in out Matrix; v : out Vectors.Vector ) is
47:
48: dv : Vectors.Vector(a'range(1)) := Divisors(a);
49:
50: begin
51: for i in a'range(1) loop
52: if (Equal(dv(i),one) or Equal(dv(i),zero))
53: then null;
54: else for j in a'range(2) loop
55: a(i,j) := a(i,j) / dv(i);
56: end loop;
57: end if;
58: end loop;
59: v := dv;
60: end Scale;
61:
62: procedure Scale ( a : in out Matrix ) is
63:
64: v : Vectors.Vector(a'range(1)) := Divisors(a);
65:
66: begin
67: for i in a'range(1) loop
68: if (Equal(v(i),one) or Equal(v(i),zero))
69: then null;
70: else for j in a'range(2) loop
71: a(i,j) := a(i,j) / v(i);
72: end loop;
73: end if;
74: end loop;
75: end Scale;
76:
77: procedure Scale ( a : in out Matrix; row,col : in integer ) is
78:
79: g : number := a(row,col);
80:
81: begin
82: if not Equal(g,one)
83: then for l in (col+1)..a'last(2) loop
84: g := gcd(g,a(row,l));
85: exit when Equal(g,one);
86: end loop;
87: end if;
88: if (not Equal(g,zero)) and (not Equal(g,one))
89: then for l in row..a'last(2) loop
90: a(row,l) := a(row,l)/g;
91: end loop;
92: end if;
93: end Scale;
94:
95: -- STATIC TRIANGULATORS :
96:
97: procedure Upper_Triangulate ( l : out Matrix; a : in out Matrix ) is
98:
99: row,pivot : integer;
100: temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
101: ll : Matrix(a'range(1),a'range(1));
102:
103: begin
104: for i in ll'range(1) loop
105: for j in ll'range(2) loop
106: Copy(zero,ll(i,j));
107: end loop;
108: Copy(one,ll(i,i));
109: end loop;
110: row := a'first(1);
111: for j in a'first(2)..a'last(2) loop
112: pivot := row-1; -- find pivot
113: for i in row..a'last(1) loop
114: if not Equal(a(i,j),zero)
115: then pivot := i;
116: exit;
117: end if;
118: end loop;
119: if pivot >= row
120: then if pivot /= row -- interchange if necessary
121: then for k in a'range(2) loop
122: Copy(a(row,k),temp);
123: Copy(a(pivot,k),a(row,k));
124: Copy(temp,a(pivot,k));
125: end loop;
126: for k in ll'range(2) loop
127: Copy(ll(row,k),temp);
128: Copy(ll(pivot,k),ll(row,k));
129: Copy(temp,ll(pivot,k));
130: end loop;
131: end if;
132: for i in (row+1)..a'last(1) loop -- make zeroes
133: if not Equal(a(i,j),zero)
134: then gcd(a(row,j),a(i,j),ka,lb,d);
135: aa := a(row,j)/d; bb := a(i,j)/d;
136: if Equal(aa,bb) and then Equal(ka,zero)
137: then Copy(lb,ka); Copy(zero,lb);
138: end if;
139: if Equal(aa,-bb) and then Equal(ka,zero)
140: then ka := -lb; Copy(zero,lb);
141: end if;
142: for k in j..a'last(2) loop
143: Copy(a(row,k),a_rowk); Clear(a(row,k));
144: -- a_rowk := a(row,k);
145: Copy(a(i,k),a_ik); Clear(a(i,k)); -- a_ik := a(i,k);
146: a(row,k) := ka*a_rowk + lb*a_ik;
147: a(i,k) := -bb*a_rowk + aa*a_ik;
148: end loop;
149: for k in ll'range(2) loop
150: a_rowk := ll(row,k);
151: a_ik := ll(i,k);
152: ll(row,k) := ka*a_rowk + lb*a_ik;
153: ll(i,k) := -bb*a_rowk + aa*a_ik;
154: end loop;
155: end if;
156: end loop;
157: row := row + 1;
158: end if;
159: exit when row > a'last(1);
160: end loop;
161: l := ll;
162: end Upper_Triangulate;
163:
164: procedure Upper_Triangulate ( a : in out Matrix ) is
165:
166: row,pivot : integer;
167: temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
168:
169: begin
170: row := a'first(1);
171: for j in a'first(2)..a'last(2) loop
172: pivot := row-1; -- find pivot
173: for i in row..a'last(1) loop
174: if not Equal(a(i,j),zero)
175: then pivot := i;
176: exit;
177: end if;
178: end loop;
179: if pivot >= row
180: then if pivot /= row -- interchange if necessary
181: then for k in a'range(2) loop
182: Copy(a(row,k),temp);
183: Copy(a(pivot,k),a(row,k));
184: Copy(temp,a(pivot,k));
185: end loop;
186: end if;
187: for i in (row+1)..a'last(1) loop -- make zeroes
188: if not Equal(a(i,j),zero)
189: then gcd(a(row,j),a(i,j),ka,lb,d);
190: aa := a(row,j)/d; bb := a(i,j)/d;
191: if Equal(aa,bb) and then Equal(ka,zero)
192: then Copy(lb,ka); Copy(zero,lb);
193: end if;
194: if Equal(aa,-bb) and then Equal(ka,zero)
195: then ka := -lb; Copy(zero,lb);
196: end if;
197: for k in j..a'last(2) loop
198: Copy(a(row,k),a_rowk); Clear(a(row,k));
199: --a_rowk := a(row,k);
200: Copy(a(i,k),a_ik); Clear(a(i,k)); --a_ik := a(i,k);
201: a(row,k) := ka*a_rowk + lb*a_ik;
202: a(i,k) := -bb*a_rowk + aa*a_ik;
203: end loop;
204: end if;
205: end loop;
206: row := row + 1;
207: end if;
208: exit when row > a'last(1);
209: end loop;
210: end Upper_Triangulate;
211:
212: procedure Upper_Triangulate
213: ( a : in out Matrix;
214: ipvt : in out Standard_Integer_Vectors.Vector ) is
215:
216: row,pivot,tmppiv : integer;
217: temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
218:
219: begin
220: row := a'first(1);
221: for j in a'first(2)..a'last(2) loop
222: pivot := row-1; -- find pivot
223: for i in row..a'last(1) loop
224: if not Equal(a(i,j),zero)
225: then pivot := i;
226: exit;
227: end if;
228: end loop;
229: if pivot >= row
230: then if pivot /= row -- interchange if necessary
231: then for k in a'range(2) loop
232: Copy(a(row,k),temp);
233: Copy(a(pivot,k),a(row,k));
234: Copy(temp,a(pivot,k));
235: end loop;
236: tmppiv := ipvt(row);
237: ipvt(row) := ipvt(pivot);
238: ipvt(pivot) := tmppiv;
239: end if;
240: for i in (row+1)..a'last(1) loop -- make zeroes
241: if not Equal(a(i,j),zero)
242: then gcd(a(row,j),a(i,j),ka,lb,d);
243: aa := a(row,j)/d; bb := a(i,j)/d;
244: if Equal(aa,bb) and then Equal(ka,zero)
245: then Copy(lb,ka); Copy(zero,lb);
246: end if;
247: if Equal(aa,-bb) and then Equal(ka,zero)
248: then ka := -lb; Copy(zero,lb);
249: end if;
250: for k in j..a'last(2) loop
251: Copy(a(row,k),a_rowk); Clear(a(row,k));
252: -- a_rowk := a(row,k);
253: Copy(a(i,k),a_ik); Clear(a(i,k)); -- a_ik := a(i,k);
254: a(row,k) := ka*a_rowk + lb*a_ik;
255: a(i,k) := -bb*a_rowk + aa*a_ik;
256: end loop;
257: end if;
258: end loop;
259: row := row + 1;
260: end if;
261: exit when row > a'last(1);
262: end loop;
263: end Upper_Triangulate;
264:
265: procedure Lower_Triangulate ( a : in out Matrix; u : out Matrix ) is
266:
267: column,pivot : integer;
268: temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
269: uu : Matrix(a'range(2),a'range(2));
270:
271: begin
272:
273: for i in uu'range(1) loop
274: for j in uu'range(2) loop
275: Copy(zero,uu(i,j));
276: end loop;
277: Copy(one,uu(i,i));
278: end loop;
279: column := a'first(2);
280: for i in a'first(1)..a'last(1) loop
281: pivot := column-1; -- find pivot
282: for j in column..a'last(2) loop
283: if not Equal(a(i,j),zero)
284: then pivot := j;
285: exit;
286: end if;
287: end loop;
288: if pivot >= column
289: then if pivot /= column -- interchange if necessary
290: then for k in a'range(1) loop
291: Copy(a(k,column),temp);
292: Copy(a(k,pivot),a(k,column));
293: Copy(temp,a(k,pivot));
294: end loop;
295: for k in uu'range(1) loop
296: Copy(uu(k,column),temp);
297: Copy(uu(k,pivot),uu(k,column));
298: Copy(temp,uu(k,pivot));
299: end loop;
300: end if;
301: for j in (column+1)..a'last(2) loop -- make zeroes
302: if not Equal(a(i,j),zero)
303: then gcd(a(i,column),a(i,j),ka,lb,d);
304: aa := a(i,column)/d; bb := a(i,j)/d;
305: if Equal(aa,bb) and then Equal(ka,zero)
306: then Copy(lb,ka); Copy(zero,lb);
307: end if;
308: if Equal(aa,-bb) and then Equal(ka,zero)
309: then ka := -lb; Copy(zero,lb);
310: end if;
311: for k in i..a'last(1) loop
312: Copy(a(k,column),a_kcolumn); Clear(a(k,column));
313: -- a_kcolumn := a(k,column);
314: Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
315: a(k,column) := a_kcolumn*ka + a_kj*lb;
316: a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
317: end loop;
318: for k in uu'range(1) loop
319: Copy(uu(k,column),a_kcolumn); Clear(uu(k,column));
320: -- a_kcolumn := uu(k,column);
321: Copy(uu(k,j),a_kj); Clear(u(k,j)); -- a_kj := uu(k,j);
322: uu(k,column) := a_kcolumn*ka + a_kj*lb;
323: uu(k,j) := a_kcolumn*(-bb) + a_kj*aa;
324: end loop;
325: end if;
326: end loop;
327: column := column + 1;
328: end if;
329: exit when column > a'last(2);
330: end loop;
331: u := uu;
332: end Lower_Triangulate;
333:
334: procedure Lower_Triangulate ( a : in out Matrix ) is
335:
336: column,pivot : integer;
337: temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
338:
339: begin
340: column := a'first(2);
341: for i in a'first(1)..a'last(1) loop
342: pivot := column-1; -- find pivot
343: for j in column..a'last(2) loop
344: if not Equal(a(i,j),zero)
345: then pivot := j;
346: exit;
347: end if;
348: end loop;
349: if pivot >= column
350: then if pivot /= column -- interchange if necessary
351: then for k in a'range(1) loop
352: Copy(a(k,column),temp);
353: Copy(a(k,pivot),a(k,column));
354: Copy(temp,a(k,pivot));
355: end loop;
356: end if;
357: for j in (column+1)..a'last(2) loop -- make zeroes
358: if not Equal(a(i,j),zero)
359: then gcd(a(i,column),a(i,j),ka,lb,d);
360: aa := a(i,column)/d; bb := a(i,j)/d;
361: if Equal(aa,bb) and then Equal(ka,zero)
362: then Copy(lb,ka); Copy(zero,lb);
363: end if;
364: if Equal(aa,-bb) and then Equal(ka,zero)
365: then ka := -lb; Copy(zero,lb);
366: end if;
367: for k in i..a'last(1) loop
368: Copy(a(k,column),a_kcolumn); Clear(a(k,column));
369: -- a_kcolumn := a(k,column);
370: Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
371: a(k,column) := a_kcolumn*ka + a_kj*lb;
372: a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
373: end loop;
374: end if;
375: end loop;
376: column := column + 1;
377: end if;
378: exit when column > a'last(2);
379: end loop;
380: end Lower_Triangulate;
381:
382: procedure Lower_Triangulate
383: ( a : in out Matrix;
384: ipvt : in out Standard_Integer_Vectors.Vector ) is
385:
386: column,pivot,tmppiv : integer;
387: temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
388:
389: begin
390: column := a'first(2);
391: for i in a'first(1)..a'last(1) loop
392: pivot := column-1; -- find pivot
393: for j in column..a'last(2) loop
394: if not Equal(a(i,j),zero)
395: then pivot := j;
396: exit;
397: end if;
398: end loop;
399: if pivot >= column
400: then if pivot /= column -- interchange if necessary
401: then for k in a'range(1) loop
402: Copy(a(k,column),temp);
403: Copy(a(k,pivot),a(k,column));
404: Copy(temp,a(k,pivot));
405: end loop;
406: tmppiv := ipvt(column);
407: ipvt(column) := ipvt(pivot);
408: ipvt(pivot) := tmppiv;
409: end if;
410: for j in (column+1)..a'last(2) loop -- make zeroes
411: if not Equal(a(i,j),zero)
412: then gcd(a(i,column),a(i,j),ka,lb,d);
413: aa := a(i,column)/d; bb := a(i,j)/d;
414: if Equal(aa,bb) and then Equal(ka,zero)
415: then Copy(lb,ka); Copy(zero,lb);
416: end if;
417: if Equal(aa,-bb) and then Equal(ka,zero)
418: then ka := -lb; Copy(zero,lb);
419: end if;
420: for k in i..a'last(1) loop
421: Copy(a(k,column),a_kcolumn); Clear(a(k,column));
422: -- a_kcolumn := a(k,column);
423: Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
424: a(k,column) := a_kcolumn*ka + a_kj*lb;
425: a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
426: end loop;
427: end if;
428: end loop;
429: column := column + 1;
430: end if;
431: exit when column > a'last(2);
432: end loop;
433: end Lower_Triangulate;
434:
435: -- SELECTORS :
436:
437: function Det ( a : Matrix ) return number is
438:
439: -- NOTE :
440: -- The triangulation is implemented independently to keep track
441: -- of row interchanges.
442:
443: res : number := one;
444: m : matrix(a'range(1),a'range(2));
445: row,pivot : integer;
446: temp,aa,bb,ka,lb,d,m_rowk,m_ik : number;
447:
448: begin
449: Copy(a,m); -- triangulate m
450: row := m'first(1);
451: for j in m'first(2)..m'last(2) loop
452: pivot := row-1; -- find pivot
453: for i in row..m'last(1) loop
454: if not Equal(m(i,j),zero)
455: then pivot := i;
456: exit;
457: end if;
458: end loop;
459: if pivot >= row
460: then if pivot /= row -- interchange if necessary
461: then for k in m'range(2) loop
462: Copy(m(row,k),temp);
463: Copy(m(pivot,k),m(row,k));
464: Copy(temp,m(pivot,k));
465: end loop;
466: Min(res);
467: end if;
468: for i in (row+1)..m'last(1) loop -- make zeroes
469: if not Equal(m(i,j),zero)
470: then gcd(m(row,j),m(i,j),ka,lb,d);
471: aa := m(row,j)/d; bb := m(i,j)/d;
472: if Equal(aa,bb) and then Equal(ka,zero)
473: then Copy(lb,ka); Copy(zero,lb);
474: end if;
475: if Equal(aa,-bb) and then Equal(ka,zero)
476: then ka := -lb; Copy(zero,lb);
477: end if;
478: for k in j..m'last(2) loop
479: Copy(m(row,k),m_rowk); Clear(m(row,k));
480: -- m_rowk := m(row,k);
481: Copy(m(i,k),m_ik); Clear(m(i,k)); -- m_ik := m(i,k);
482: m(row,k) := ka*m_rowk + lb*m_ik;
483: m(i,k) := -bb*m_rowk + aa*m_ik;
484: end loop;
485: end if;
486: end loop;
487: row := row + 1;
488: end if;
489: exit when row > m'last(1);
490: end loop;
491: for k in m'range(1) loop
492: Mul(res,m(k,k));
493: end loop;
494: Clear(m);
495: return res;
496: end Det;
497:
498: function Per ( i,n : natural; a : matrix;
499: kk : Standard_Integer_Vectors.Vector ) return number is
500: begin
501: if i = n+1
502: then return one;
503: else declare
504: res,acc : number;
505: kkk : Standard_Integer_Vectors.Vector(kk'range) := kk;
506: begin
507: Copy(zero,res);
508: for j in kk'range loop
509: if not Equal(a(i,j),zero) and then (kk(j) /= 0)
510: then kkk(j) := kkk(j) - 1;
511: acc := Per(i+1,n,a,kkk);
512: Add(res,acc);
513: Clear(acc);
514: kkk(j) := kkk(j) + 1;
515: end if;
516: end loop;
517: return res;
518: end;
519: end if;
520: end Per;
521:
522: function Per ( i,n : natural; a : matrix;
523: kk : Standard_Integer_Vectors.Vector; max : number )
524: return number is
525: begin
526: if i = n+1
527: then return one;
528: else declare
529: res,acc : number;
530: kkk : Standard_Integer_Vectors.Vector(kk'range) := kk;
531: begin
532: Copy(zero,res);
533: for j in kk'range loop
534: if not Equal(a(i,j),zero) and then (kk(j) /= 0)
535: then kkk(j) := kkk(j) - 1;
536: acc := a(i,j)*Per(i+1,n,a,kkk,max);
537: Add(res,acc);
538: Clear(acc);
539: kkk(j) := kkk(j) + 1;
540: end if;
541: exit when (res > max) or Equal(res,max);
542: end loop;
543: return res;
544: end;
545: end if;
546: end Per;
547:
548: function Per ( a : matrix; k : Standard_Integer_Vectors.Vector )
549: return number is
550:
551: -- ALGORITHM :
552: -- Row expansion without memory, as developed by C.W. Wampler,
553: -- see `Bezout Number Calculations for Multi-Homogeneous Polynomial
554: -- Systems', Appl. Math. Comput. 51:(2-3), 143-157, 1992.
555:
556: begin
557: return Per(1,a'last(1),a,k);
558: end Per;
559:
560: function Per ( a : matrix; k : Standard_Integer_Vectors.Vector;
561: max : number ) return number is
562:
563: -- ALGORITHM :
564: -- Row expansion without memory, as developed by C.W. Wampler,
565: -- see `Bezout Number Calculations for Multi-Homogeneous Polynomial
566: -- Systems', Appl. Math. Comput. 51:(2-3), 143-157, 1992.
567:
568: begin
569: return Per(1,a'last(1),a,k,max);
570: end Per;
571:
572: function Rank ( a : Matrix ) return natural is
573:
574: res : natural := 0;
575: m : Matrix(a'range(1),a'range(2));
576: column : integer;
577:
578: begin
579: Copy(a,m);
580: Upper_Triangulate(m);
581: -- compute the length of chain of nonzero elements in m :
582: -- search first nonzero element in first row of m :
583: column := m'first(2)-1;
584: for k in m'range(2) loop
585: if not Equal(m(m'first(1),k),zero)
586: then column := k;
587: end if;
588: exit when (column = k);
589: end loop;
590: if column < m'first(2)
591: then return 0; -- all elements of m are zero
592: else for k in m'range(1) loop
593: exit when column > m'last(2);
594: if not Equal(m(k,column),zero)
595: then res := res + 1;
596: else -- search for next nonzero element on row k :
597: for l in column+1..m'last(2) loop
598: if not Equal(m(k,l),zero)
599: then column := l;
600: res := res + 1;
601: end if;
602: exit when (column = l);
603: end loop;
604: end if;
605: column := column + 1;
606: end loop;
607: end if;
608: Clear(m);
609: return res;
610: end Rank;
611:
612: -- DYNAMIC TRIANGULATOR :
613:
614: procedure Triangulate ( l : in integer; m : in out matrix;
615: ipvt : in out Standard_Integer_Vectors.vector;
616: piv : out integer ) is
617:
618: -- DESCRIPTION :
619: -- Updates lth row of m such that m remains upper triangular.
620:
621: tmp,a,b,lcmab,faca,facb : number;
622: pivot,index,tmppiv : integer;
623:
624: begin
625: Switch(ipvt,l,m); -- pivoting for previous unknowns
626: index := 1; -- update : make l-1 zeroes in row l
627: for k in 1..(l-1) loop
628: if (not Equal(m(l,index),zero))
629: and then (not Equal(m(k,index),zero)) -- make m(l,index) zero
630: then a := m(k,index); b := m(l,index);
631: lcmab := lcm(a,b);
632: if lcmab < zero then lcmab := -lcmab; end if;
633: facb := lcmab/b; faca := lcmab/a;
634: if facb > zero
635: then for i in index..m'last(2) loop
636: m(l,i) := facb*m(l,i) - faca*m(k,i);
637: end loop;
638: else for i in index..m'last(2) loop
639: m(l,i) := -facb*m(l,i) + faca*m(k,i);
640: end loop;
641: end if;
642: end if;
643: if not Equal(m(k,index),zero)
644: then index := index + 1;
645: end if;
646: end loop;
647: pivot := 0; -- search pivot
648: for k in l..m'last(2)-1 loop
649: if not Equal(m(l,k),zero)
650: then pivot := k;
651: end if;
652: exit when pivot /= 0;
653: end loop;
654: if pivot > l
655: then for k in 1..l loop -- interchange columns l and pivot
656: Copy(m(k,l),tmp);
657: Copy(m(k,pivot),m(k,l));
658: Copy(tmp,m(k,pivot));
659: end loop;
660: tmppiv := ipvt(l);
661: ipvt(l) := ipvt(pivot);
662: ipvt(pivot) := tmppiv;
663: end if;
664: piv := pivot;
665: end Triangulate;
666:
667: procedure Switch ( ipvt : in Standard_Integer_Vectors.vector;
668: index : in integer; m : in out matrix ) is
669:
670: tmpv : Vectors.Vector(m'range(2));
671:
672: begin
673: for k in tmpv'range loop
674: Copy(m(index,k),tmpv(k));
675: end loop;
676: for k in tmpv'range loop
677: Copy(tmpv(ipvt(k)),m(index,k));
678: end loop;
679: Vectors.Clear(tmpv);
680: end Switch;
681:
682: procedure Switch ( ipvt : in Standard_Integer_Vectors.vector;
683: first,last : in integer; m : in out matrix) is
684:
685: tmpv : Vectors.Vector(m'range(2));
686:
687: begin
688: for index in first..last loop
689: for k in tmpv'range loop
690: Copy(m(index,k),tmpv(k));
691: end loop;
692: for k in tmpv'range loop
693: Copy(tmpv(ipvt(k)),m(index,k));
694: end loop;
695: Vectors.Clear(tmpv);
696: end loop;
697: end Switch;
698:
699: procedure Switch ( l,pivot,index : in integer; m : in out matrix ) is
700:
701: tmp : number;
702:
703: begin
704: if l /= pivot
705: then Copy(m(index,l),tmp);
706: Copy(m(index,pivot),m(index,l));
707: Copy(tmp,m(index,pivot));
708: end if;
709: end Switch;
710:
711: procedure Switch ( l,pivot : in integer;
712: first,last : in integer; m : in out matrix ) is
713:
714: tmp : number;
715:
716: begin
717: if l /= pivot
718: then for index in first..last loop
719: Copy(m(index,l),tmp);
720: Copy(m(index,pivot),m(index,l));
721: Copy(tmp,m(index,pivot));
722: end loop;
723: end if;
724: end Switch;
725:
726: -- SOLVERS :
727:
728: function Check0 ( a : Matrix; x : Vectors.Vector ) return boolean is
729:
730: -- DESCRIPTION :
731: -- Returns true if x is a solution of the system a*x = 0.
732:
733: tmp : Vectors.Vector(a'range(1));
734:
735: begin
736: tmp := a*x;
737: for i in tmp'range loop
738: if not Equal(tmp(i),zero)
739: then Vectors.Clear(tmp); return false;
740: end if;
741: end loop;
742: Vectors.Clear(tmp);
743: return true;
744: end Check0;
745:
746: procedure Solve0 ( a : in Matrix; x : in out Vectors.Vector ) is
747:
748: -- ALGORITHM :
749: -- An intermediate, generating matrix tmp will be constructed,
750: -- such that
751: -- 1) the solution x to tmp*x = 0 is the same of a*x = 0;
752: -- 2) tmp(i,i) and tmp(i,ind) are the only nonzero entries.
753: -- Before this construction, it will be checked whether there
754: -- exists a zero column and the index ind must be determined.
755: -- After the definition of tmp, the back substitution process
756: -- yields a solution.
757:
758: piv,ind : integer;
759: tmp : Matrix(a'first(1)..(a'last(1)+1),a'range(2));
760: res : Vectors.Vector(x'range);
761: pivots : Standard_Integer_Vectors.Vector(x'range);
762: zero_column : boolean;
763:
764: begin
765: -- initialization of tmp,ind and pivots :
766: for i in tmp'range(1) loop
767: for j in tmp'range(2) loop
768: Copy(zero,tmp(i,j));
769: end loop;
770: end loop;
771: for i in pivots'range loop
772: pivots(i) := i;
773: end loop;
774: ind := x'first(1)-1;
775: for i in a'range(1) loop
776: piv := pivots'first-1;
777: for j in a'range(2) loop
778: if not Equal(a(i,j),zero)
779: then piv := pivots(j);
780: pivots(j) := pivots(i);
781: pivots(i) := piv;
782: exit;
783: end if;
784: end loop;
785: zero_column := true;
786: for j in a'first(1)..i loop
787: Copy(a(j,pivots(i)),tmp(j,i));
788: if zero_column and then not Equal(tmp(j,i),zero)
789: then zero_column := false;
790: end if;
791: end loop;
792: if piv < pivots'first or else zero_column or else Equal(tmp(i,i),zero)
793: then ind := i; exit;
794: end if;
795: end loop;
796: if zero_column
797: then for i in x'range loop
798: Copy(zero,x(i));
799: end loop;
800: Copy(one,x(ind));
801: elsif (ind < x'first(1)) and (a'last(1) >= a'last(2))
802: then for i in x'range loop
803: Copy(zero,x(i));
804: end loop;
805: else
806: if ind < x'first(1)
807: then ind := a'last(1)+1;
808: for j in tmp'range(2) loop
809: Copy(zero,tmp(ind,j));
810: end loop;
811: zero_column := true;
812: for j in a'range(1) loop
813: Copy(a(j,pivots(ind)),tmp(j,ind));
814: if zero_column and then not Equal(tmp(j,ind),zero)
815: then zero_column := false;
816: end if;
817: end loop;
818: end if;
819: if zero_column
820: then for i in x'range loop
821: Copy(zero,x(i));
822: end loop;
823: Copy(one,x(ind));
824: else
825: -- construct generating matrix :
826: for i in reverse (tmp'first(2)+1)..(ind-1) loop -- i = column
827: for k in tmp'first(1)..(i-1) loop
828: if not Equal(tmp(k,i),zero) -- make tmp(k,i) zero
829: then declare
830: aa,bb,d : number;
831: begin
832: d := gcd(tmp(i,i),tmp(k,i));
833: aa := tmp(i,i)/d; bb := tmp(k,i)/d;
834: for l in k..(i-1) loop
835: Mul(tmp(k,l),aa);
836: end loop;
837: Copy(zero,tmp(k,i)); --aa*tmp(k,i) - bb*tmp(i,i);
838: tmp(k,ind) := aa*tmp(k,ind) - bb*tmp(i,ind);
839: end;
840: Scale(tmp,k,k); -- to avoid numeric_error
841: end if;
842: end loop; -- upper half of ith colum consists of zero entries
843: end loop;
844: -- generate x by back substitution :
845: x(ind) := tmp(x'first,x'first);
846: for i in (x'first+1)..(ind-1) loop
847: if not Equal(tmp(i,i),zero)
848: then x(ind) := lcm(tmp(i,i),x(ind));
849: end if;
850: end loop;
851: for i in x'first..(ind-1) loop
852: if Equal(tmp(i,i),zero)
853: then Copy(zero,x(i));
854: else x(i) := -(tmp(i,ind)*x(ind))/tmp(i,i);
855: end if;
856: end loop;
857: end if;
858: end if;
859: for i in res'range loop -- take pivots into account
860: Copy(zero,res(i));
861: end loop;
862: for i in x'first..ind loop
863: Copy(x(i),res(pivots(i)));
864: end loop;
865: Vectors.Copy(res,x);
866: end Solve0;
867:
868: procedure Solve1 ( a : in Matrix; x : in out Vectors.Vector;
869: b : in Vectors.Vector; fail : out boolean ) is
870:
871: acc : number;
872:
873: begin
874: fail := false;
875: for i in reverse x'range loop
876: Copy(b(i),x(i));
877: for j in (i+1)..x'last loop
878: acc := a(i,j)*x(j);
879: Sub(x(i),acc);
880: Clear(acc);
881: end loop;
882: if not Equal(x(i),zero) and then not Equal(a(i,i),zero)
883: then acc := Rmd(x(i),a(i,i));
884: if Equal(acc,zero)
885: then Div(x(i),a(i,i));
886: else fail := true;
887: end if;
888: Clear(acc);
889: if fail
890: then Vectors.Clear(x); return;
891: end if;
892: end if;
893: end loop;
894: end Solve1;
895:
896: procedure Solve1 ( a : in Matrix; b : in out Vectors.Vector;
897: fail : out boolean ) is
898:
899: acc : number;
900:
901: begin
902: fail := false;
903: for i in reverse b'range loop
904: for j in (i+1)..b'last loop
905: acc := a(i,j)*b(j);
906: Sub(b(i),acc);
907: Clear(acc);
908: end loop;
909: if not Equal(b(i),zero) and then not Equal(a(i,i),zero)
910: then acc := Rmd(b(i),a(i,i));
911: if Equal(acc,zero)
912: then Div(b(i),a(i,i));
913: else fail := true;
914: end if;
915: Clear(acc);
916: if fail
917: then Vectors.Clear(b); return;
918: end if;
919: end if;
920: end loop;
921: end Solve1;
922:
923: function Solve ( m : Matrix; ipvt : Standard_Integer_Vectors.Vector )
924: return Vectors.Vector is
925:
926: x,res : Vectors.Vector(ipvt'range);
927: a : matrix(m'first(1)..m'last(1)-1,m'range(2));
928: ind : integer := a'first(1); -- index for the current row number
929: cnt0 : natural := 0; -- counts the number of zero rows
930:
931: begin
932: for k in a'range(1) loop
933: if not Equal(m(k,k),zero) -- otherwise : skip zero row !
934: then for l in a'range(2) loop
935: Copy(m(k,l),a(ind,l));
936: end loop;
937: ind := ind + 1;
938: else for l in a'range(2) loop
939: Copy(m(k,l),a(a'last(1) - cnt0,l));
940: end loop;
941: cnt0 := cnt0 + 1;
942: end if;
943: end loop;
944: for i in x'range loop
945: Copy(zero,x(i));
946: end loop;
947: Solve0(a,x);
948: for k in res'range loop
949: res(ipvt(k)) := x(k);
950: end loop;
951: if res(res'last) < zero
952: then Vectors.Min(res);
953: end if;
954: Clear(a);
955: return res;
956: end Solve;
957:
958: end Generic_Integer_Linear_Solvers;