Annotation of OpenXM/src/Ti/Rsimp.c, Revision 1.1
1.1 ! maekawa 1: /*
! 2: ** Rsimp.c Birk Huber, 4/99
! 3: ** -- definitions of linear programming data structure and
! 4: ** -- basic implementation of revised simplex method.
! 5: **
! 6: **
! 7: ** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
! 8: ** copyright (c) 1999 Birk Huber
! 9: **
! 10: */
! 11: #include <stdio.h>
! 12: #include <math.h>
! 13: #define RSIMP_H 1
! 14: #include "matrices.h"
! 15: #include "Rsimp.h"
! 16: #define verbose 0
! 17: #define zero_tol RS_zt
! 18:
! 19: int LP_MAX_N=0,LP_MAX_M=0;
! 20: int LP_N=0, LP_M=0;
! 21: double **LP_A=0;
! 22: double *LP_B=0;
! 23: double *LP_C=0;
! 24: double *LP_X=0;
! 25: int *LP_Basis=0;
! 26: int *LP_NonBasis=0;
! 27: double **LP_Q=0;
! 28: double **LP_R=0;
! 29: double *LP_t1=0;
! 30: double *LP_t2=0;
! 31: double RS_zt=0.0000001;
! 32:
! 33:
! 34: #define basis0(j) (basis[j])
! 35: #define nonbasis0(j) (nonbasis[j])
! 36: #define AN0(i,j) (A[nonbasis0(j)][i])
! 37: #define AB0(i,j) (A[basis0(j)][i])
! 38: #define CB0(i) (c[basis0(i)])
! 39: #define CN0(i) (c[nonbasis0(i)])
! 40: #define XB0(i) (x[basis0(i)])
! 41: #define XN0(i) (x[nonbasis0(i)])
! 42: #define Y0(i) (t1[i])
! 43: #define W0(i) (t2[i])
! 44: #define D0(i) (t2[i])
! 45: #define R0(i,j) (R[j][i])
! 46: #define Q0(i,j) (Q[i][j])
! 47:
! 48:
! 49: /*
! 50: ** void LP_free_space():
! 51: ** deallocate space for LP data structures.
! 52: ** set all LP globals to 0.
! 53: **
! 54: */
! 55: void LP_free_space(){
! 56: free_matrix(LP_A);
! 57: free_matrix(LP_Q);
! 58: free_matrix(LP_R);
! 59: free_vector(LP_B);
! 60: free_vector(LP_C);
! 61: free_vector(LP_X);
! 62: free_ivector(LP_Basis);
! 63: free_ivector(LP_NonBasis);
! 64: free_vector(LP_t1);
! 65: free_vector(LP_t2);
! 66: LP_M=LP_N=LP_MAX_M=LP_N=0;
! 67: LP_A=LP_Q=LP_R=0;
! 68: LP_B =LP_C=LP_X=LP_t1=LP_t2=0;
! 69: LP_Basis=LP_NonBasis = 0;
! 70: }
! 71:
! 72: /*
! 73: ** void LP_get_space(int M, int N):
! 74: ** allocate space for LP data structures.
! 75: **
! 76: */
! 77: void LP_get_space(int M, int N){
! 78: if (M>LP_MAX_M || N>LP_MAX_N){
! 79: LP_free_space();
! 80: LP_A=new_matrix(M,N);
! 81: LP_Q=new_matrix(M,M);
! 82: LP_R=new_matrix(M,N);
! 83: LP_B = new_vector(M);
! 84: LP_C = new_vector(N);
! 85: LP_X = new_vector(N);
! 86: LP_Basis=new_ivector(M);
! 87: LP_NonBasis = new_ivector(N-M);
! 88: LP_t1=new_vector(N);
! 89: LP_t2=new_vector(N);
! 90: LP_MAX_M=M;
! 91: LP_MAX_N=N;
! 92: }
! 93: LP_M=M;
! 94: LP_N=N;
! 95: }
! 96:
! 97: /*
! 98: ** void Print_LP():
! 99: ** print LP data structures to stdout.
! 100: **
! 101: */
! 102: void Print_LP(){
! 103: int i,j;
! 104: fprintf(stdout,"LP_C \n ");
! 105: for(i=0;i<LP_N;i++)fprintf(stdout," %f ;", LP_C[i] );
! 106: fprintf(stdout,"\nLP_A | LP_B\n");
! 107: for(j=0;j<LP_M;j++) {
! 108: for(i=0;i<LP_N;i++) fprintf(stdout," %f ",LP_A(j,i));
! 109: fprintf(stdout,"%f \n",LP_B[j]);
! 110: }
! 111: fprintf(stdout, "LP_X \n ");
! 112: for(i=0;i<LP_N;i++) fprintf(stdout," %f ",LP_X[i]);
! 113: fprintf(stdout, "\nLP_Basis \n ");
! 114: for(i=0;i<LP_M;i++) fprintf(stdout," %d ", LP_Basis[i]);
! 115: fprintf(stdout,"\nLP_LP_NonBasis \n ");
! 116: for(i=0;i<LP_N-LP_M;i++) fprintf(stdout," %d ",LP_NonBasis[i]);
! 117: fprintf(stdout,"\n\n");
! 118: }
! 119:
! 120: /*
! 121: **
! 122: **int Rsimp(m,n,A,b,c,x,basis,nonbasis,R,Q,t1,t2):
! 123: ** revised simplex method (Using Bland's rule)
! 124: ** and a qr factorization to solve the linear equations
! 125: **
! 126: ** Adapted from algorithms presented in
! 127: ** Linear Approximations and Extensions
! 128: ** (theory and algorithms)
! 129: ** Fang & Puthenpura
! 130: ** Prentice Hall, Engelwood Cliffs NJ (1993)
! 131: ** and
! 132: ** Linear Programming
! 133: ** Chvatal
! 134: ** Freeman and Company, New York, 1983
! 135: **
! 136: ** (developed first in Octave, many thanks to the author)
! 137: **
! 138: **
! 139: ** Solve the problem
! 140: ** minimize C'x,
! 141: ** subject to A*x=b, x>=0
! 142: ** for x,c,b n-vectors, and A an m,n matrix with full row rank
! 143: **
! 144: ** Assumptions:
! 145: ** A mxn matrix with full row rank.
! 146: ** b an m matrix.
! 147: ** c an n-vector.
! 148: ** x an n-vector holding a basic feasible solution,
! 149: ** basis m-vector holding indices of the basic variables in x
! 150: ** nonbasis n-m vector holding the indices not appearing in x.
! 151: **
! 152: ** Returns:
! 153: ** LP_FAIL if algorithm doesn't terminate.
! 154: ** LP_UNBD if problem is unbounded
! 155: ** LP_OPT if optimum found
! 156: ** efects:
! 157: ** A,b,c unchanged.
! 158: ** x basis, nonbasis, hold info describing last basic feasible
! 159: ** solution.
! 160: ** Q,R hold qrdecomp of last basis matrix.
! 161: ** t1,t2 undefined.
! 162: **
! 163: **
! 164: */
! 165:
! 166: int Rsimp(int m, int n, double **A, double *b, double *c,
! 167: double *x, int *basis, int *nonbasis,
! 168: double **R, double **Q, double *t1, double *t2){
! 169: int i,j,k,l,q,qv;
! 170: int max_steps=20;
! 171: double r,a,at;
! 172: int ll;
! 173: void GQR(int,int,double**,double**);
! 174: max_steps=4*n;
! 175:
! 176: for(k=0; k<=max_steps;k++){
! 177: /*
! 178: ++ Step 0) load new basis matrix and factor it
! 179: */
! 180: for(i=0;i<m;i++)for(j=0;j<m;j++)R0(i,j)=AB0(i,j);
! 181: GQR(m,m,Q,R);
! 182:
! 183: /*
! 184: ++ Step 1) solving system B'*w=c(basis)
! 185: ++ a) forward solve R'*y=c(basis)
! 186: */
! 187: for(i=0;i<m;i++){
! 188: Y0(i)=0.0;
! 189: for(j=0;j<i;j++)Y0(i)+=R0(j,i)*Y0(j);
! 190: if (R0(i,i)!=0.0) Y0(i)=(CB0(i)-Y0(i))/R0(i,i);
! 191: else {
! 192: printf("Warning Singular Matrix Found\n");
! 193: return LP_FAIL;
! 194: }
! 195: }
! 196: /*
! 197: ++ b) find w=Q*y
! 198: ++ note: B'*w=(Q*R)'*Q*y= R'*(Q'*Q)*y=R'*y=c(basis)
! 199: */
! 200: for(i=0;i<m;i++){
! 201: W0(i)=0.0;
! 202: for(j=0;j<m;j++)W0(i)+=Q0(i,j)*Y0(j);
! 203: }
! 204:
! 205: /*
! 206: ++ Step 2)find entering variable,
! 207: ++ (use lexicographically first variable with negative reduced cost)
! 208: */
! 209: q=n;
! 210: for(i=0;i<n-m;i++){
! 211: /* calculate reduced cost */
! 212: r=CN0(i);
! 213: for(j=0;j<m;j++) r-=W0(j)*AN0(j,i);
! 214: if (r<-zero_tol && (q==n || nonbasis0(i)<nonbasis0(q))) q=i;
! 215: }
! 216:
! 217: /*
! 218: ++ if ratios were all nonnegative current solution is optimal
! 219: */
! 220: if (q==n){
! 221: if (verbose>0) printf("optimal solution found in %d iterations\n",k);
! 222: return LP_OPT;
! 223: }
! 224: /*
! 225: ++ Step 3)Calculate translation direction for q entering
! 226: ++ by solving system B*d=-A(:,nonbasis(q));
! 227: ++ a) let y=-Q'*A(:,nonbasis(q));
! 228: */
! 229: for(i=0;i<m;i++){
! 230: Y0(i)=0.0;
! 231: for(j=0;j<m;j++) Y0(i)-=Q0(j,i)*AN0(j,q);
! 232: }
! 233:
! 234: /*
! 235: ++ b) back solve Rd=y (d=R\y)
! 236: ++ note B*d= Q*R*d=Q*y=Q*-Q'*A(:nonbasis(q))=-A(:,nonbasis(q))
! 237: */
! 238: for(i=m-1;i>=0;i--){
! 239: D0(i)=0.0;
! 240: for(j=m-1;j>=i+1;j--)D0(i)+=R0(i,j)*D0(j);
! 241: if (R0(i,i)!=0.0) D0(i)=(Y0(i)-D0(i))/R0(i,i);
! 242: else {
! 243: printf("Warning Singular Matrix Found\n");
! 244: return LP_FAIL;
! 245: }
! 246: }
! 247: /*
! 248: ++ Step 4 Choose leaving variable
! 249: ++ (first variable to become negative, by moving in direction D)
! 250: ++ (if none become negative, then objective function unbounded)
! 251: */
! 252: a=0;
! 253: l=-1;
! 254: for(i=0;i<m;i++){
! 255: if (D0(i)<-zero_tol){
! 256: at=-1*XB0(i)/D0(i);
! 257: if (l==-1 || at<a){ a=at; l=i;}
! 258: }
! 259: }
! 260: if (l==-1){
! 261: if (verbose>0){
! 262: printf("Objective function Unbounded (%d iterations)\n",k);
! 263: }
! 264: return LP_UNBD;
! 265: }
! 266: /*
! 267: ++ Step 5) Update solution and basis data
! 268: */
! 269: XN0(q)=a;
! 270: for(j=0;j<m;j++) XB0(j)+=a*D0(j);
! 271: XB0(l)=0.0; /* enforce strict zeroness of nonbasis variables */
! 272: qv=nonbasis0(q);
! 273: nonbasis0(q)=basis0(l);
! 274: basis0(l)=qv;
! 275: }
! 276: if (verbose>=0){
! 277: printf("Simplex Algorithm did not Terminate in %d iterations\n",k);
! 278: }
! 279: return LP_FAIL;
! 280: }
! 281:
! 282:
! 283: /*
! 284: **void GQR(int r, int c, double **Q, double **R):
! 285: ** Use givens rotations on R to bring it into triangular form.
! 286: ** Store orthogonal matrix needed to bring R to triangular form in Q.
! 287: ** Assume R is an rxc matrix and Q is rxr.
! 288: ** (Question: is r>=c needed ?)
! 289: **
! 290: */
! 291: void GQR(int r, int c, double **Q, double **R){
! 292: int i,j,k;
! 293: double s,s1,s2;
! 294: double t1,t2;
! 295:
! 296: for(i=0;i<r;i++){
! 297: for(k=0;k<r;k++)Q0(i,k)=0.0;
! 298: Q0(i,i)=1.0;
! 299: }
! 300:
! 301: for (i=0;i<c;i++)
! 302: for (k=i+1;k<r;k++)
! 303: /* performing givens rotations to zero A[k][i] */
! 304: if (R0(k,i)!=0){
! 305: s=sqrt(R0(i,i)*R0(i,i)+R0(k,i)*R0(k,i));
! 306: s1=R0(i,i)/s;
! 307: s2=R0(k,i)/s;
! 308: for(j=0;j<c;j++) {
! 309: t1=R0(i,j);
! 310: t2=R0(k,j);
! 311: R0(i,j)=s1*t1+s2*t2;
! 312: R0(k,j)=-s2*t1+s1*t2;
! 313: }
! 314: /* actually doing givens row rotations on Q */
! 315: for(j=0;j<r;j++){
! 316: t1=Q0(j,i);
! 317: t2=Q0(j,k);
! 318: Q0(j,i)=s1*t1+s2*t2;
! 319: Q0(j,k)=-s2*t1+s1*t2;
! 320: }
! 321: }
! 322:
! 323: }
! 324:
! 325:
! 326:
! 327:
! 328:
! 329:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Rsimp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / Ti