Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/floating_linear_inequalities.adb, Revision 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