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>
![[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