Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/generic_integer_linear_solvers.adb, Revision 1.1
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;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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