[BACK]Return to Rsimp.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / Ti

Annotation of OpenXM/src/Ti/Rsimp.c, Revision 1.1.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>