1 #include <math.h>
2 #include <stddef.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include "f2c.h"
6
7 /* prototype for function at end (from TOMS, via Netlib and f2c) */
8
9 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
10 integer *klmd, integer *klm2d, integer *nklmd,
11 integer *n2d, real *q,
12 integer *kode, real *toler, integer *iter, real *x,
13 real *res, real * error, real *cu, integer *iu, integer *s) ;
14
15 /*---------------------------------------------------------------------
16 Approximately (L1) solve equations
17
18 j=nvec-1
19 z[i] = SUM A[j][i] * y[j] (i=0..ndim-1)
20 j=0
21
22 for y[j] (j=0..nvec-1), subject to constraints (based on y[j] input)
23
24 If input cony != 0, then you can supply constraints on the values
25 of the output y[j] by putting values into the input y[j]:
26
27 input y[j] = 0 ==> unconstrained
28 = 1 ==> y[j] must be non-negative on output
29 = -1 ==> y[j] must be non-positive on output
30
31 If cony == 0, then the input y[j] is ignored.
32
33 The return value of the function is E = sum |z[i]-SUM[i]| >= 0
34 if everything worked, and is a negative number if an error occured.
35 -----------------------------------------------------------------------*/
36
cl1_solve(int ndim,int nvec,float * z,float ** A,float * y,int cony)37 float cl1_solve( int ndim, int nvec, float *z, float **A, float *y, int cony )
38 {
39 /* loop counters */
40
41 int jj , ii ;
42
43 /* variables for CL1 (types declared in f2c.h) */
44
45 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
46 real *q , toler , *x , *res , error , *cu ;
47
48 /*-- check inputs --*/
49
50 if( ndim < 1 || nvec < 1 ) kode = 4 ;
51 if( A == NULL || y == NULL || z == NULL ) kode = 4 ;
52 for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
53
54 if( kode ){
55 fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
56 return (float)(-kode) ;
57 }
58
59 /*-- setup call to CL1 --*/
60
61 k = ndim ;
62 l = 0 ; /* no linear equality constraints */
63 m = 0 ; /* no linear inequality constraints */
64 n = nvec ;
65
66 klmd = k+l+m ;
67 klm2d = k+l+m+2 ;
68 nklmd = n+k+l+m ;
69 n2d = n+2 ;
70
71 kode = (cony != 0) ; /* enforce implicit constraints on x[] */
72 iter = 11*klmd ;
73
74 toler = 0.0001 ;
75 error = 0.0 ;
76
77 /* input/output matrices & vectors */
78
79 q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
80 x = (real *) calloc( 1, sizeof(real) * n2d ) ;
81 res = (real *) calloc( 1, sizeof(real) * klmd ) ;
82
83 /* workspaces */
84
85 cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ;
86 iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
87 s = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
88
89 /* load matrices & vectors */
90
91 for( jj=0 ; jj < nvec ; jj++ )
92 for( ii=0 ; ii < ndim ; ii++ )
93 q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */
94
95 for( ii=0 ; ii < ndim ; ii++ )
96 q[ii+nvec*klm2d] = z[ii] ; /* vector */
97
98 if( cony ){
99 for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */
100 x[jj] = (y[jj] < 0.0) ? -1.0
101 :(y[jj] > 0.0) ? 1.0 : 0.0 ;
102 }
103
104 for( ii=0 ; ii < klmd ; ii++ ) /* no constraints on resids */
105 res[ii] = 0.0 ;
106
107 /*-- do the work --*/
108
109 cl1_fort( &k, &l, &m, &n,
110 &klmd, &klm2d, &nklmd, &n2d,
111 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
112
113 free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
114
115 if( kode != 0 ){
116 free(x) ;
117 #if 0
118 switch( kode ){
119 case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
120 case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break;
121 case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break;
122 default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break;
123 }
124 #endif
125 return (float)(-kode) ;
126 }
127
128 /*-- copy results into output --*/
129
130 for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
131
132 free(x) ; return (float)error ;
133 }
134
135 /*---------------------------------------------------------------------
136 Approximately (L1) solve equations
137
138 j=nvec-1
139 z[i] = SUM A[j][i] * y[j] (i=0..ndim-1)
140 j=0
141
142 for y[j] (j=0..nvec-1), subject to constraints on the signs of y[j]
143 (based on y[j] input) AND subject to constraints on the signs of
144 the residuals z[i]-SUM[i] (based on rez[i] input).
145
146 If input cony != 0, then you can supply constraints on the values
147 of the output y[j] by putting values into the input y[j]:
148
149 input y[j] = 0 ==> unconstrained
150 = 1 ==> y[j] must be non-negative on output
151 = -1 ==> y[j] must be non-positive on output
152
153 If cony == 0, then the input y[j] is ignored.
154
155 If input conr != 0, then you can supply constraints on the values
156 of the residuals z[i]-SUM[i] by putting values into the input rez[i]:
157
158 input rez[i] = 0 ==> unconstrained
159 = 1 ==> z[i]-SUM[i] must be non-negative
160 = -1 ==> z[i]-SUM[i] must be non-positive
161
162 If conr == 0, then the input values in rez[] are ignored.
163
164 The outputs are in y[] and rez[].
165
166 The return value of the function is E = sum |z[i]-SUM[i]| >= 0
167 if everything worked, and is a negative number if an error occured.
168 -----------------------------------------------------------------------*/
169
cl1_solve_res(int ndim,int nvec,float * z,float ** A,float * y,int cony,float * rez,int conr)170 float cl1_solve_res( int ndim, int nvec, float *z, float **A,
171 float *y, int cony , float *rez , int conr )
172 {
173 /* loop counters */
174
175 int jj , ii ;
176
177 /* variables for CL1 (types declared in f2c.h) */
178
179 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
180 real *q , toler , *x , *res , error , *cu ;
181
182 /*-- check inputs --*/
183
184 if( ndim < 1 || nvec < 1 ) kode = 4 ;
185 if( A == NULL || y == NULL || z == NULL ) kode = 4 ;
186 for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
187
188 if( kode ){
189 fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
190 return (float)(-kode) ;
191 }
192
193 /*-- setup call to CL1 --*/
194
195 k = ndim ;
196 l = 0 ; /* no linear equality constraints */
197 m = 0 ; /* no linear inequality constraints */
198 n = nvec ;
199
200 klmd = k+l+m ;
201 klm2d = k+l+m+2 ;
202 nklmd = n+k+l+m ;
203 n2d = n+2 ;
204
205 kode = (cony != 0 || conr != 0) ; /* enforce implicit constraints on x[] */
206 iter = 11*klmd ;
207
208 toler = 0.0001 ;
209 error = 0.0 ;
210
211 /* input/output matrices & vectors */
212
213 q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
214 x = (real *) calloc( 1, sizeof(real) * n2d ) ;
215 res = (real *) calloc( 1, sizeof(real) * klmd ) ;
216
217 /* workspaces */
218
219 cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ;
220 iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
221 s = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
222
223 /* load matrices & vectors */
224
225 for( jj=0 ; jj < nvec ; jj++ )
226 for( ii=0 ; ii < ndim ; ii++ )
227 q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */
228
229 for( ii=0 ; ii < ndim ; ii++ )
230 q[ii+nvec*klm2d] = z[ii] ; /* vector */
231
232 if( cony ){
233 for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */
234 x[jj] = (y[jj] < 0.0) ? -1.0
235 :(y[jj] > 0.0) ? 1.0 : 0.0 ;
236 }
237
238 if( conr ){
239 for( ii=0 ; ii < ndim ; ii++ ) /* constraints on resids */
240 res[ii] = (rez[ii] < 0.0) ? -1.0
241 :(rez[ii] > 0.0) ? 1.0 : 0.0 ;
242 }
243
244 /*-- do the work --*/
245
246 cl1_fort( &k, &l, &m, &n,
247 &klmd, &klm2d, &nklmd, &n2d,
248 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
249
250 free(q) ; free(cu) ; free(iu) ; free(s) ;
251
252 if( kode != 0 ){
253 free(x) ; free(res) ;
254 #if 0
255 switch( kode ){
256 case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
257 case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break;
258 case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break;
259 default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break;
260 }
261 #endif
262 return (float)(-kode) ;
263 }
264
265 /*-- copy results into output --*/
266
267 for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
268
269 for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = (float) res[ii] ;
270
271 free(res); free(x); return (float)error;
272 }
273
274 #if 0
275 /*---------------------------------------------------------------------
276 Find a set of coefficients ec[] so that
277
278 j=nvec-1 { }
279 S[i] = base_vec[i] + SUM { ec[j] * extra_vec[j][i] }
280 j=0 { }
281
282 is non-negative for i=0..ndim-1, and so that
283
284 i=ndim-1 { }
285 SUM { S[i] }
286 i=0 { }
287
288 is as small as possible.
289
290 The return value of the function is 0 if everything worked, and
291 nonzero if an error occured. If the function worked, then the
292 coefficients are returned in ec[] (which must be allocated by
293 the caller).
294
295 Method: uses the CL1 subroutine from the TOMS library;
296 the f2c translation appears at the end of this file
297
298 N.B.: NOT TESTED
299 -----------------------------------------------------------------------*/
300
301 int cl1_pos_sum( int ndim , int nvec ,
302 float * base_vec , float ** extra_vec , float * ec )
303 {
304 /* internal variables */
305
306 int jj , ii ;
307
308 /* variables for CL1 */
309
310 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode,iter , *iu,*s ;
311 real *q , toler , *x , *res , error , *cu ;
312
313 /*-- check inputs --*/
314
315 if( ndim < 1 || nvec < 1 ) return 4 ;
316 if( base_vec == NULL || extra_vec == NULL || ec == NULL ) return 4 ;
317 for( jj=0 ; jj < nvec ; jj++ ) if( extra_vec[jj] == NULL ) return 4 ;
318
319 /*-- setup call to CL1 --*/
320
321 k = ndim ;
322 l = 0 ;
323 m = 0 ;
324 n = nvec ;
325
326 klmd = k+l+m ;
327 klm2d = k+l+m+2 ;
328 nklmd = n+k+l+m ;
329 n2d = n+2 ;
330
331 kode = 1 ;
332 iter = 10*klmd ;
333
334 toler = 0.0001 ;
335 error = 0.0 ;
336
337 /* input/output matrices & vectors */
338
339 q = (real *) malloc( sizeof(real) * klm2d*n2d ) ;
340 x = (real *) malloc( sizeof(real) * n2d ) ;
341 res = (real *) malloc( sizeof(real) * klmd ) ;
342
343 /* workspaces */
344
345 cu = (real *) malloc( sizeof(real) * 2*nklmd ) ;
346 iu = (integer *) malloc( sizeof(integer) * 2*nklmd ) ;
347 s = (integer *) malloc( sizeof(integer) * klmd ) ;
348
349 /* load matrices & vectors */
350
351 for( jj=0 ; jj < nvec ; jj++ )
352 for( ii=0 ; ii < ndim ; ii++ )
353 q[ii+jj*klm2d] = extra_vec[jj][ii] ;
354
355 for( ii=0 ; ii < ndim ; ii++ )
356 q[ii+nvec*klm2d] = base_vec[ii] ;
357
358 for( jj=0 ; jj < nvec ; jj++ ) x[jj] = 0.0 ;
359
360 for( ii=0 ; ii < ndim ; ii++ ) res[ii] = 1.0 ;
361
362 /*-- do the work --*/
363
364 cl1_fort( &k, &l, &m, &n,
365 &klmd, &klm2d, &nklmd, &n2d,
366 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
367
368 free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
369
370 if( kode != 0 ){ free(x) ; return (int)kode ; }
371
372 /*-- copy results into output --*/
373
374 for( jj=0 ; jj < nvec ; jj++ ) ec[jj] = (float)(-x[jj]) ;
375
376 free(x) ; return 0 ;
377 }
378 #endif
379
380 /*======================================================================
381 Translated from the TOMS library CL1 algorithm
382 ========================================================================*/
383
384 /* cl1.f -- translated by f2c (version 19970805).
385 You must link the resulting object file with the libraries:
386 -lf2c -lm (in that order)
387 */
388
cl1_fort(integer * k,integer * l,integer * m,integer * n,integer * klmd,integer * klm2d,integer * nklmd,integer * n2d,real * q,integer * kode,real * toler,integer * iter,real * x,real * res,real * error,real * cu,integer * iu,integer * s)389 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
390 integer *klmd, integer *klm2d, integer *nklmd, integer *n2d, real *q,
391 integer *kode, real *toler, integer *iter, real *x, real *res, real *
392 error, real *cu, integer *iu, integer *s)
393 {
394 /* System generated locals */
395 integer q_dim1, q_offset, i__1, i__2;
396 real r__1;
397
398 /* Local variables ['static' removed by RWCox] */
399 integer iimn, nklm;
400 real xmin, xmax;
401 integer iout=0, i__, j;
402 real z__;
403 integer iineg, maxit, n1, n2;
404 real pivot;
405 integer ia, ii, kk, in=0, nk, js;
406 real sn;
407 integer iphase, kforce;
408 real zu, zv;
409 integer nk1;
410 real tpivot;
411 integer klm, jmn, nkl, jpn;
412 real cuv;
413 doublereal sum;
414 integer klm1, klm2, nkl1;
415
416
417 /* THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX */
418 /* METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION */
419 /* TO A K BY N SYSTEM OF LINEAR EQUATIONS */
420 /* AX=B */
421 /* SUBJECT TO L LINEAR EQUALITY CONSTRAINTS */
422 /* CX=D */
423 /* AND M LINEAR INEQUALITY CONSTRAINTS */
424 /* EX.LE.F. */
425
426 /* DESCRIPTION OF PARAMETERS */
427
428 /* K NUMBER OF ROWS OF THE MATRIX A (K.GE.1). */
429 /* L NUMBER OF ROWS OF THE MATRIX C (L.GE.0). */
430 /* M NUMBER OF ROWS OF THE MATRIX E (M.GE.0). */
431 /* N NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1). */
432 /* KLMD SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS. */
433 /* KLM2D SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS. */
434 /* NKLMD SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS. */
435 /* N2D SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS */
436
437 /* Q TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND */
438 /* AT LEAST N2D COLUMNS. */
439 /* ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS */
440 /* B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS */
441 /* AND N+1 COLUMNS OF Q AS FOLLOWS */
442 /* A B */
443 /* Q = C D */
444 /* E F */
445 /* THESE VALUES ARE DESTROYED BY THE SUBROUTINE. */
446
447 /* KODE A CODE USED ON ENTRY TO, AND EXIT */
448 /* FROM, THE SUBROUTINE. */
449 /* ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0. */
450 /* HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS */
451 /* ARE TO BE INCLUDED IMPLICITLY, RATHER THAN */
452 /* EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE */
453 /* SHOULD BE SET TO 1, AND THE NONNEGATIVITY */
454 /* CONSTRAINTS INCLUDED IN THE ARRAYS X AND */
455 /* RES (SEE BELOW). */
456 /* ON EXIT, KODE HAS ONE OF THE */
457 /* FOLLOWING VALUES */
458 /* 0- OPTIMAL SOLUTION FOUND, */
459 /* 1- NO FEASIBLE SOLUTION TO THE */
460 /* CONSTRAINTS, */
461 /* 2- CALCULATIONS TERMINATED */
462 /* PREMATURELY DUE TO ROUNDING ERRORS, */
463 /* 3- MAXIMUM NUMBER OF ITERATIONS REACHED. */
464
465 /* TOLER A SMALL POSITIVE TOLERANCE. EMPIRICAL */
466 /* EVIDENCE SUGGESTS TOLER = 10**(-D*2/3), */
467 /* WHERE D REPRESENTS THE NUMBER OF DECIMAL */
468 /* DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY, */
469 /* THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO */
470 /* AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED */
471 /* TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY */
472 /* NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER. */
473
474 /* ITER ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON */
475 /* THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. */
476 /* A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER */
477 /* GIVES THE NUMBER OF SIMPLEX ITERATIONS. */
478
479 /* X ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D. */
480 /* ON EXIT THIS ARRAY CONTAINS A */
481 /* SOLUTION TO THE L1 PROBLEM. IF KODE=1 */
482 /* ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE */
483 /* SIMPLE NONNEGATIVITY CONSTRAINTS ON THE */
484 /* VARIABLES. THE VALUES -1, 0, OR 1 */
485 /* FOR X(J) INDICATE THAT THE J-TH VARIABLE */
486 /* IS RESTRICTED TO BE .LE.0, UNRESTRICTED, */
487 /* OR .GE.0 RESPECTIVELY. */
488
489 /* RES ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD. */
490 /* ON EXIT THIS CONTAINS THE RESIDUALS B-AX */
491 /* IN THE FIRST K COMPONENTS, D-CX IN THE */
492 /* NEXT L COMPONENTS (THESE WILL BE =0),AND */
493 /* F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON */
494 /* ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE */
495 /* NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS */
496 /* B-AX. THE VALUES -1, 0, OR 1 FOR RES(I) */
497 /* INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS */
498 /* RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0 */
499 /* RESPECTIVELY. */
500
501 /* ERROR ON EXIT, THIS GIVES THE MINIMUM SUM OF */
502 /* ABSOLUTE VALUES OF THE RESIDUALS. */
503
504 /* CU A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND */
505 /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
506
507 /* IU A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND */
508 /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
509
510 /* S INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR */
511 /* WORKSPACE. */
512
513 /* IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO */
514 /* DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY */
515 /* THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN */
516 /* EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE */
517 /* FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS. */
518 /* SUBROUTINE COL(V1, V2, XMLT, NOTROW, K) */
519 /* THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE */
520 /* VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW). */
521 /* DIMENSION V1(K), V2(K) */
522 /* KEND = NOTROW - 1 */
523 /* KSTART = NOTROW + 1 */
524 /* IF (KEND .LT. 1) GO TO 20 */
525 /* DO 10 I=1,KEND */
526 /* V1(I) = V1(I) + XMLT*V2(I) */
527 /* 10 CONTINUE */
528 /* IF(KSTART .GT. K) GO TO 40 */
529 /* 20 DO 30 I=KSTART,K */
530 /* V1(I) = V1(I) + XMLT*V2(I) */
531 /* 30 CONTINUE */
532 /* 40 RETURN */
533 /* END */
534 /* SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR */
535 /* INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION. */
536 /* -----------------------------------------------------------------------
537 */
538
539 /* INITIALIZATION. */
540
541 /* Parameter adjustments */
542 --s;
543 --res;
544 iu -= 3;
545 cu -= 3;
546 --x;
547 q_dim1 = *klm2d;
548 q_offset = q_dim1 + 1;
549 q -= q_offset;
550
551 /* Function Body */
552 maxit = *iter;
553 n1 = *n + 1;
554 n2 = *n + 2;
555 nk = *n + *k;
556 nk1 = nk + 1;
557 nkl = nk + *l;
558 nkl1 = nkl + 1;
559 klm = *k + *l + *m;
560 klm1 = klm + 1;
561 klm2 = klm + 2;
562 nklm = *n + klm;
563 kforce = 1;
564 *iter = 0;
565 js = 1;
566 ia = 0;
567 /* SET UP LABELS IN Q. */
568 i__1 = *n;
569 for (j = 1; j <= i__1; ++j) {
570 q[klm2 + j * q_dim1] = (real) j;
571 /* L10: */
572 }
573 i__1 = klm;
574 for (i__ = 1; i__ <= i__1; ++i__) {
575 q[i__ + n2 * q_dim1] = (real) (*n + i__);
576 if (q[i__ + n1 * q_dim1] >= 0.f) {
577 goto L30;
578 }
579 i__2 = n2;
580 for (j = 1; j <= i__2; ++j) {
581 q[i__ + j * q_dim1] = -q[i__ + j * q_dim1];
582 /* L20: */
583 }
584 L30:
585 ;
586 }
587 /* SET UP PHASE 1 COSTS. */
588 iphase = 2;
589 i__1 = nklm;
590 for (j = 1; j <= i__1; ++j) {
591 cu[(j << 1) + 1] = 0.f;
592 cu[(j << 1) + 2] = 0.f;
593 iu[(j << 1) + 1] = 0;
594 iu[(j << 1) + 2] = 0;
595 /* L40: */
596 }
597 if (*l == 0) {
598 goto L60;
599 }
600 i__1 = nkl;
601 for (j = nk1; j <= i__1; ++j) {
602 cu[(j << 1) + 1] = 1.f;
603 cu[(j << 1) + 2] = 1.f;
604 iu[(j << 1) + 1] = 1;
605 iu[(j << 1) + 2] = 1;
606 /* L50: */
607 }
608 iphase = 1;
609 L60:
610 if (*m == 0) {
611 goto L80;
612 }
613 i__1 = nklm;
614 for (j = nkl1; j <= i__1; ++j) {
615 cu[(j << 1) + 2] = 1.f;
616 iu[(j << 1) + 2] = 1;
617 jmn = j - *n;
618 if (q[jmn + n2 * q_dim1] < 0.f) {
619 iphase = 1;
620 }
621 /* L70: */
622 }
623 L80:
624 if (*kode == 0) {
625 goto L150;
626 }
627 i__1 = *n;
628 for (j = 1; j <= i__1; ++j) {
629 if ((r__1 = x[j]) < 0.f) {
630 goto L90;
631 } else if (r__1 == 0) {
632 goto L110;
633 } else {
634 goto L100;
635 }
636 L90:
637 cu[(j << 1) + 1] = 1.f;
638 iu[(j << 1) + 1] = 1;
639 goto L110;
640 L100:
641 cu[(j << 1) + 2] = 1.f;
642 iu[(j << 1) + 2] = 1;
643 L110:
644 ;
645 }
646 i__1 = *k;
647 for (j = 1; j <= i__1; ++j) {
648 jpn = j + *n;
649 if ((r__1 = res[j]) < 0.f) {
650 goto L120;
651 } else if (r__1 == 0) {
652 goto L140;
653 } else {
654 goto L130;
655 }
656 L120:
657 cu[(jpn << 1) + 1] = 1.f;
658 iu[(jpn << 1) + 1] = 1;
659 if (q[j + n2 * q_dim1] > 0.f) {
660 iphase = 1;
661 }
662 goto L140;
663 L130:
664 cu[(jpn << 1) + 2] = 1.f;
665 iu[(jpn << 1) + 2] = 1;
666 if (q[j + n2 * q_dim1] < 0.f) {
667 iphase = 1;
668 }
669 L140:
670 ;
671 }
672 L150:
673 if (iphase == 2) {
674 goto L500;
675 }
676 /* COMPUTE THE MARGINAL COSTS. */
677 L160:
678 i__1 = n1;
679 for (j = js; j <= i__1; ++j) {
680 sum = 0.;
681 i__2 = klm;
682 for (i__ = 1; i__ <= i__2; ++i__) {
683 ii = q[i__ + n2 * q_dim1];
684 if (ii < 0) {
685 goto L170;
686 }
687 z__ = cu[(ii << 1) + 1];
688 goto L180;
689 L170:
690 iineg = -ii;
691 z__ = cu[(iineg << 1) + 2];
692 L180:
693 sum += (doublereal) q[i__ + j * q_dim1] * (doublereal) z__;
694 /* L190: */
695 }
696 q[klm1 + j * q_dim1] = sum;
697 /* L200: */
698 }
699 i__1 = *n;
700 for (j = js; j <= i__1; ++j) {
701 ii = q[klm2 + j * q_dim1];
702 if (ii < 0) {
703 goto L210;
704 }
705 z__ = cu[(ii << 1) + 1];
706 goto L220;
707 L210:
708 iineg = -ii;
709 z__ = cu[(iineg << 1) + 2];
710 L220:
711 q[klm1 + j * q_dim1] -= z__;
712 /* L230: */
713 }
714 /* DETERMINE THE VECTOR TO ENTER THE BASIS. */
715 L240:
716 xmax = 0.f;
717 if (js > *n) {
718 goto L490;
719 }
720 i__1 = *n;
721 for (j = js; j <= i__1; ++j) {
722 zu = q[klm1 + j * q_dim1];
723 ii = q[klm2 + j * q_dim1];
724 if (ii > 0) {
725 goto L250;
726 }
727 ii = -ii;
728 zv = zu;
729 zu = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
730 goto L260;
731 L250:
732 zv = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
733 L260:
734 if (kforce == 1 && ii > *n) {
735 goto L280;
736 }
737 if (iu[(ii << 1) + 1] == 1) {
738 goto L270;
739 }
740 if (zu <= xmax) {
741 goto L270;
742 }
743 xmax = zu;
744 in = j;
745 L270:
746 if (iu[(ii << 1) + 2] == 1) {
747 goto L280;
748 }
749 if (zv <= xmax) {
750 goto L280;
751 }
752 xmax = zv;
753 in = j;
754 L280:
755 ;
756 }
757 if (xmax <= *toler) {
758 goto L490;
759 }
760 if (q[klm1 + in * q_dim1] == xmax) {
761 goto L300;
762 }
763 i__1 = klm2;
764 for (i__ = 1; i__ <= i__1; ++i__) {
765 q[i__ + in * q_dim1] = -q[i__ + in * q_dim1];
766 /* L290: */
767 }
768 q[klm1 + in * q_dim1] = xmax;
769 /* DETERMINE THE VECTOR TO LEAVE THE BASIS. */
770 L300:
771 if (iphase == 1 || ia == 0) {
772 goto L330;
773 }
774 xmax = 0.f;
775 i__1 = ia;
776 for (i__ = 1; i__ <= i__1; ++i__) {
777 z__ = (r__1 = q[i__ + in * q_dim1], dabs(r__1));
778 if (z__ <= xmax) {
779 goto L310;
780 }
781 xmax = z__;
782 iout = i__;
783 L310:
784 ;
785 }
786 if (xmax <= *toler) {
787 goto L330;
788 }
789 i__1 = n2;
790 for (j = 1; j <= i__1; ++j) {
791 z__ = q[ia + j * q_dim1];
792 q[ia + j * q_dim1] = q[iout + j * q_dim1];
793 q[iout + j * q_dim1] = z__;
794 /* L320: */
795 }
796 iout = ia;
797 --ia;
798 pivot = q[iout + in * q_dim1];
799 goto L420;
800 L330:
801 kk = 0;
802 i__1 = klm;
803 for (i__ = 1; i__ <= i__1; ++i__) {
804 z__ = q[i__ + in * q_dim1];
805 if (z__ <= *toler) {
806 goto L340;
807 }
808 ++kk;
809 res[kk] = q[i__ + n1 * q_dim1] / z__;
810 s[kk] = i__;
811 L340:
812 ;
813 }
814 L350:
815 if (kk > 0) {
816 goto L360;
817 }
818 *kode = 2;
819 goto L590;
820 L360:
821 xmin = res[1];
822 iout = s[1];
823 j = 1;
824 if (kk == 1) {
825 goto L380;
826 }
827 i__1 = kk;
828 for (i__ = 2; i__ <= i__1; ++i__) {
829 if (res[i__] >= xmin) {
830 goto L370;
831 }
832 j = i__;
833 xmin = res[i__];
834 iout = s[i__];
835 L370:
836 ;
837 }
838 res[j] = res[kk];
839 s[j] = s[kk];
840 L380:
841 --kk;
842 pivot = q[iout + in * q_dim1];
843 ii = q[iout + n2 * q_dim1];
844 if (iphase == 1) {
845 goto L400;
846 }
847 if (ii < 0) {
848 goto L390;
849 }
850 if (iu[(ii << 1) + 2] == 1) {
851 goto L420;
852 }
853 goto L400;
854 L390:
855 iineg = -ii;
856 if (iu[(iineg << 1) + 1] == 1) {
857 goto L420;
858 }
859 L400:
860 ii = abs(ii);
861 cuv = cu[(ii << 1) + 1] + cu[(ii << 1) + 2];
862 if (q[klm1 + in * q_dim1] - pivot * cuv <= *toler) {
863 goto L420;
864 }
865 /* BYPASS INTERMEDIATE VERTICES. */
866 i__1 = n1;
867 for (j = js; j <= i__1; ++j) {
868 z__ = q[iout + j * q_dim1];
869 q[klm1 + j * q_dim1] -= z__ * cuv;
870 q[iout + j * q_dim1] = -z__;
871 /* L410: */
872 }
873 q[iout + n2 * q_dim1] = -q[iout + n2 * q_dim1];
874 goto L350;
875 /* GAUSS-JORDAN ELIMINATION. */
876 L420:
877 if (*iter < maxit) {
878 goto L430;
879 }
880 *kode = 3;
881 goto L590;
882 L430:
883 ++(*iter);
884 i__1 = n1;
885 for (j = js; j <= i__1; ++j) {
886 if (j != in) {
887 q[iout + j * q_dim1] /= pivot;
888 }
889 /* L440: */
890 }
891 /* IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION */
892 /* SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN */
893 /* TO AND INCLUDING STATEMENT NUMBER 460 BY.. */
894 /* DO 460 J=JS,N1 */
895 /* IF(J .EQ. IN) GO TO 460 */
896 /* Z = -Q(IOUT,J) */
897 /* CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1) */
898 /* 460 CONTINUE */
899 i__1 = n1;
900 for (j = js; j <= i__1; ++j) {
901 if (j == in) {
902 goto L460;
903 }
904 z__ = -q[iout + j * q_dim1];
905 i__2 = klm1;
906 for (i__ = 1; i__ <= i__2; ++i__) {
907 if (i__ != iout) {
908 q[i__ + j * q_dim1] += z__ * q[i__ + in * q_dim1];
909 }
910 /* L450: */
911 }
912 L460:
913 ;
914 }
915 tpivot = -pivot;
916 i__1 = klm1;
917 for (i__ = 1; i__ <= i__1; ++i__) {
918 if (i__ != iout) {
919 q[i__ + in * q_dim1] /= tpivot;
920 }
921 /* L470: */
922 }
923 q[iout + in * q_dim1] = 1.f / pivot;
924 z__ = q[iout + n2 * q_dim1];
925 q[iout + n2 * q_dim1] = q[klm2 + in * q_dim1];
926 q[klm2 + in * q_dim1] = z__;
927 ii = dabs(z__);
928 if (iu[(ii << 1) + 1] == 0 || iu[(ii << 1) + 2] == 0) {
929 goto L240;
930 }
931 i__1 = klm2;
932 for (i__ = 1; i__ <= i__1; ++i__) {
933 z__ = q[i__ + in * q_dim1];
934 q[i__ + in * q_dim1] = q[i__ + js * q_dim1];
935 q[i__ + js * q_dim1] = z__;
936 /* L480: */
937 }
938 ++js;
939 goto L240;
940 /* TEST FOR OPTIMALITY. */
941 L490:
942 if (kforce == 0) {
943 goto L580;
944 }
945 if (iphase == 1 && q[klm1 + n1 * q_dim1] <= *toler) {
946 goto L500;
947 }
948 kforce = 0;
949 goto L240;
950 /* SET UP PHASE 2 COSTS. */
951 L500:
952 iphase = 2;
953 i__1 = nklm;
954 for (j = 1; j <= i__1; ++j) {
955 cu[(j << 1) + 1] = 0.f;
956 cu[(j << 1) + 2] = 0.f;
957 /* L510: */
958 }
959 i__1 = nk;
960 for (j = n1; j <= i__1; ++j) {
961 cu[(j << 1) + 1] = 1.f;
962 cu[(j << 1) + 2] = 1.f;
963 /* L520: */
964 }
965 i__1 = klm;
966 for (i__ = 1; i__ <= i__1; ++i__) {
967 ii = q[i__ + n2 * q_dim1];
968 if (ii > 0) {
969 goto L530;
970 }
971 ii = -ii;
972 if (iu[(ii << 1) + 2] == 0) {
973 goto L560;
974 }
975 cu[(ii << 1) + 2] = 0.f;
976 goto L540;
977 L530:
978 if (iu[(ii << 1) + 1] == 0) {
979 goto L560;
980 }
981 cu[(ii << 1) + 1] = 0.f;
982 L540:
983 ++ia;
984 i__2 = n2;
985 for (j = 1; j <= i__2; ++j) {
986 z__ = q[ia + j * q_dim1];
987 q[ia + j * q_dim1] = q[i__ + j * q_dim1];
988 q[i__ + j * q_dim1] = z__;
989 /* L550: */
990 }
991 L560:
992 ;
993 }
994 goto L160;
995 L570:
996 if (q[klm1 + n1 * q_dim1] <= *toler) {
997 goto L500;
998 }
999 *kode = 1;
1000 goto L590;
1001 L580:
1002 if (iphase == 1) {
1003 goto L570;
1004 }
1005 /* PREPARE OUTPUT. */
1006 *kode = 0;
1007 L590:
1008 sum = 0.;
1009 i__1 = *n;
1010 for (j = 1; j <= i__1; ++j) {
1011 x[j] = 0.f;
1012 /* L600: */
1013 }
1014 i__1 = klm;
1015 for (i__ = 1; i__ <= i__1; ++i__) {
1016 res[i__] = 0.f;
1017 /* L610: */
1018 }
1019 i__1 = klm;
1020 for (i__ = 1; i__ <= i__1; ++i__) {
1021 ii = q[i__ + n2 * q_dim1];
1022 sn = 1.f;
1023 if (ii > 0) {
1024 goto L620;
1025 }
1026 ii = -ii;
1027 sn = -1.f;
1028 L620:
1029 if (ii > *n) {
1030 goto L630;
1031 }
1032 x[ii] = sn * q[i__ + n1 * q_dim1];
1033 goto L640;
1034 L630:
1035 iimn = ii - *n;
1036 res[iimn] = sn * q[i__ + n1 * q_dim1];
1037 if (ii >= n1 && ii <= nk) {
1038 sum += (doublereal) q[i__ + n1 * q_dim1];
1039 }
1040 L640:
1041 ;
1042 }
1043 *error = sum;
1044 return 0;
1045 } /* cl1_ */
1046
1047 /*===================================================================
1048 SUBROUTINE CL1(K, L, M, N, KLMD, KLM2D, NKLMD, N2D,
1049 * Q, KODE, TOLER, ITER, X, RES, ERROR, CU, IU, S)
1050 C
1051 C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX
1052 C METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION
1053 C TO A K BY N SYSTEM OF LINEAR EQUATIONS
1054 C AX=B
1055 C SUBJECT TO L LINEAR EQUALITY CONSTRAINTS
1056 C CX=D
1057 C AND M LINEAR INEQUALITY CONSTRAINTS
1058 C EX.LE.F.
1059 C
1060 C DESCRIPTION OF PARAMETERS
1061 C
1062 C K NUMBER OF ROWS OF THE MATRIX A (K.GE.1).
1063 C L NUMBER OF ROWS OF THE MATRIX C (L.GE.0).
1064 C M NUMBER OF ROWS OF THE MATRIX E (M.GE.0).
1065 C N NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1).
1066 C KLMD SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS.
1067 C KLM2D SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS.
1068 C NKLMD SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS.
1069 C N2D SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS
1070 C
1071 C Q TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND
1072 C AT LEAST N2D COLUMNS.
1073 C ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS
1074 C B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS
1075 C AND N+1 COLUMNS OF Q AS FOLLOWS
1076 C A B
1077 C Q = C D
1078 C E F
1079 C THESE VALUES ARE DESTROYED BY THE SUBROUTINE.
1080 C
1081 C KODE A CODE USED ON ENTRY TO, AND EXIT
1082 C FROM, THE SUBROUTINE.
1083 C ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0.
1084 C HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS
1085 C ARE TO BE INCLUDED IMPLICITLY, RATHER THAN
1086 C EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE
1087 C SHOULD BE SET TO 1, AND THE NONNEGATIVITY
1088 C CONSTRAINTS INCLUDED IN THE ARRAYS X AND
1089 C RES (SEE BELOW).
1090 C ON EXIT, KODE HAS ONE OF THE
1091 C FOLLOWING VALUES
1092 C 0- OPTIMAL SOLUTION FOUND,
1093 C 1- NO FEASIBLE SOLUTION TO THE
1094 C CONSTRAINTS,
1095 C 2- CALCULATIONS TERMINATED
1096 C PREMATURELY DUE TO ROUNDING ERRORS,
1097 C 3- MAXIMUM NUMBER OF ITERATIONS REACHED.
1098 C
1099 C TOLER A SMALL POSITIVE TOLERANCE. EMPIRICAL
1100 C EVIDENCE SUGGESTS TOLER = 10**(-D*2/3),
1101 C WHERE D REPRESENTS THE NUMBER OF DECIMAL
1102 C DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY,
1103 C THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO
1104 C AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED
1105 C TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY
1106 C NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER.
1107 C
1108 C ITER ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON
1109 C THE MAXIMUM NUMBER OF ITERATIONS ALLOWED.
1110 C A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER
1111 C GIVES THE NUMBER OF SIMPLEX ITERATIONS.
1112 C
1113 C X ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D.
1114 C ON EXIT THIS ARRAY CONTAINS A
1115 C SOLUTION TO THE L1 PROBLEM. IF KODE=1
1116 C ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE
1117 C SIMPLE NONNEGATIVITY CONSTRAINTS ON THE
1118 C VARIABLES. THE VALUES -1, 0, OR 1
1119 C FOR X(J) INDICATE THAT THE J-TH VARIABLE
1120 C IS RESTRICTED TO BE .LE.0, UNRESTRICTED,
1121 C OR .GE.0 RESPECTIVELY.
1122 C
1123 C RES ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD.
1124 C ON EXIT THIS CONTAINS THE RESIDUALS B-AX
1125 C IN THE FIRST K COMPONENTS, D-CX IN THE
1126 C NEXT L COMPONENTS (THESE WILL BE =0),AND
1127 C F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON
1128 C ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE
1129 C NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS
1130 C B-AX. THE VALUES -1, 0, OR 1 FOR RES(I)
1131 C INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS
1132 C RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0
1133 C RESPECTIVELY.
1134 C
1135 C ERROR ON EXIT, THIS GIVES THE MINIMUM SUM OF
1136 C ABSOLUTE VALUES OF THE RESIDUALS.
1137 C
1138 C CU A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND
1139 C AT LEAST NKLMD COLUMNS USED FOR WORKSPACE.
1140 C
1141 C IU A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND
1142 C AT LEAST NKLMD COLUMNS USED FOR WORKSPACE.
1143 C
1144 C S INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR
1145 C WORKSPACE.
1146 C
1147 C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO
1148 C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY
1149 C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN
1150 C EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE
1151 C FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS.
1152 C SUBROUTINE COL(V1, V2, XMLT, NOTROW, K)
1153 C THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE
1154 C VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW).
1155 C DIMENSION V1(K), V2(K)
1156 C KEND = NOTROW - 1
1157 C KSTART = NOTROW + 1
1158 C IF (KEND .LT. 1) GO TO 20
1159 C DO 10 I=1,KEND
1160 C V1(I) = V1(I) + XMLT*V2(I)
1161 C 10 CONTINUE
1162 C IF(KSTART .GT. K) GO TO 40
1163 C 20 DO 30 I=KSTART,K
1164 C V1(I) = V1(I) + XMLT*V2(I)
1165 C 30 CONTINUE
1166 C 40 RETURN
1167 C END
1168 C SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR
1169 C INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION.
1170 C-----------------------------------------------------------------------
1171 DOUBLE PRECISION SUM
1172 DOUBLE PRECISION DBLE
1173 REAL Q, X, Z, CU, SN, ZU, ZV, CUV, RES, XMAX, XMIN,
1174 * ERROR, PIVOT, TOLER, TPIVOT
1175 REAL ABS
1176 INTEGER I, J, K, L, M, N, S, IA, II, IN, IU, JS, KK,
1177 * NK, N1, N2, JMN, JPN, KLM, NKL, NK1, N2D, IIMN,
1178 * IOUT, ITER, KLMD, KLM1, KLM2, KODE, NKLM, NKL1,
1179 * KLM2D, MAXIT, NKLMD, IPHASE, KFORCE, IINEG
1180 INTEGER IABS
1181 DIMENSION Q(KLM2D,N2D), X(N2D), RES(KLMD),
1182 * CU(2,NKLMD), IU(2,NKLMD), S(KLMD)
1183 C
1184 C INITIALIZATION.
1185 C
1186 MAXIT = ITER
1187 N1 = N + 1
1188 N2 = N + 2
1189 NK = N + K
1190 NK1 = NK + 1
1191 NKL = NK + L
1192 NKL1 = NKL + 1
1193 KLM = K + L + M
1194 KLM1 = KLM + 1
1195 KLM2 = KLM + 2
1196 NKLM = N + KLM
1197 KFORCE = 1
1198 ITER = 0
1199 JS = 1
1200 IA = 0
1201 C SET UP LABELS IN Q.
1202 DO 10 J=1,N
1203 Q(KLM2,J) = J
1204 10 CONTINUE
1205 DO 30 I=1,KLM
1206 Q(I,N2) = N + I
1207 IF (Q(I,N1).GE.0.) GO TO 30
1208 DO 20 J=1,N2
1209 Q(I,J) = -Q(I,J)
1210 20 CONTINUE
1211 30 CONTINUE
1212 C SET UP PHASE 1 COSTS.
1213 IPHASE = 2
1214 DO 40 J=1,NKLM
1215 CU(1,J) = 0.
1216 CU(2,J) = 0.
1217 IU(1,J) = 0
1218 IU(2,J) = 0
1219 40 CONTINUE
1220 IF (L.EQ.0) GO TO 60
1221 DO 50 J=NK1,NKL
1222 CU(1,J) = 1.
1223 CU(2,J) = 1.
1224 IU(1,J) = 1
1225 IU(2,J) = 1
1226 50 CONTINUE
1227 IPHASE = 1
1228 60 IF (M.EQ.0) GO TO 80
1229 DO 70 J=NKL1,NKLM
1230 CU(2,J) = 1.
1231 IU(2,J) = 1
1232 JMN = J - N
1233 IF (Q(JMN,N2).LT.0.) IPHASE = 1
1234 70 CONTINUE
1235 80 IF (KODE.EQ.0) GO TO 150
1236 DO 110 J=1,N
1237 IF (X(J)) 90, 110, 100
1238 90 CU(1,J) = 1.
1239 IU(1,J) = 1
1240 GO TO 110
1241 100 CU(2,J) = 1.
1242 IU(2,J) = 1
1243 110 CONTINUE
1244 DO 140 J=1,K
1245 JPN = J + N
1246 IF (RES(J)) 120, 140, 130
1247 120 CU(1,JPN) = 1.
1248 IU(1,JPN) = 1
1249 IF (Q(J,N2).GT.0.0) IPHASE = 1
1250 GO TO 140
1251 130 CU(2,JPN) = 1.
1252 IU(2,JPN) = 1
1253 IF (Q(J,N2).LT.0.0) IPHASE = 1
1254 140 CONTINUE
1255 150 IF (IPHASE.EQ.2) GO TO 500
1256 C COMPUTE THE MARGINAL COSTS.
1257 160 DO 200 J=JS,N1
1258 SUM = 0.D0
1259 DO 190 I=1,KLM
1260 II = Q(I,N2)
1261 IF (II.LT.0) GO TO 170
1262 Z = CU(1,II)
1263 GO TO 180
1264 170 IINEG = -II
1265 Z = CU(2,IINEG)
1266 180 SUM = SUM + DBLE(Q(I,J))*DBLE(Z)
1267 190 CONTINUE
1268 Q(KLM1,J) = SUM
1269 200 CONTINUE
1270 DO 230 J=JS,N
1271 II = Q(KLM2,J)
1272 IF (II.LT.0) GO TO 210
1273 Z = CU(1,II)
1274 GO TO 220
1275 210 IINEG = -II
1276 Z = CU(2,IINEG)
1277 220 Q(KLM1,J) = Q(KLM1,J) - Z
1278 230 CONTINUE
1279 C DETERMINE THE VECTOR TO ENTER THE BASIS.
1280 240 XMAX = 0.
1281 IF (JS.GT.N) GO TO 490
1282 DO 280 J=JS,N
1283 ZU = Q(KLM1,J)
1284 II = Q(KLM2,J)
1285 IF (II.GT.0) GO TO 250
1286 II = -II
1287 ZV = ZU
1288 ZU = -ZU - CU(1,II) - CU(2,II)
1289 GO TO 260
1290 250 ZV = -ZU - CU(1,II) - CU(2,II)
1291 260 IF (KFORCE.EQ.1 .AND. II.GT.N) GO TO 280
1292 IF (IU(1,II).EQ.1) GO TO 270
1293 IF (ZU.LE.XMAX) GO TO 270
1294 XMAX = ZU
1295 IN = J
1296 270 IF (IU(2,II).EQ.1) GO TO 280
1297 IF (ZV.LE.XMAX) GO TO 280
1298 XMAX = ZV
1299 IN = J
1300 280 CONTINUE
1301 IF (XMAX.LE.TOLER) GO TO 490
1302 IF (Q(KLM1,IN).EQ.XMAX) GO TO 300
1303 DO 290 I=1,KLM2
1304 Q(I,IN) = -Q(I,IN)
1305 290 CONTINUE
1306 Q(KLM1,IN) = XMAX
1307 C DETERMINE THE VECTOR TO LEAVE THE BASIS.
1308 300 IF (IPHASE.EQ.1 .OR. IA.EQ.0) GO TO 330
1309 XMAX = 0.
1310 DO 310 I=1,IA
1311 Z = ABS(Q(I,IN))
1312 IF (Z.LE.XMAX) GO TO 310
1313 XMAX = Z
1314 IOUT = I
1315 310 CONTINUE
1316 IF (XMAX.LE.TOLER) GO TO 330
1317 DO 320 J=1,N2
1318 Z = Q(IA,J)
1319 Q(IA,J) = Q(IOUT,J)
1320 Q(IOUT,J) = Z
1321 320 CONTINUE
1322 IOUT = IA
1323 IA = IA - 1
1324 PIVOT = Q(IOUT,IN)
1325 GO TO 420
1326 330 KK = 0
1327 DO 340 I=1,KLM
1328 Z = Q(I,IN)
1329 IF (Z.LE.TOLER) GO TO 340
1330 KK = KK + 1
1331 RES(KK) = Q(I,N1)/Z
1332 S(KK) = I
1333 340 CONTINUE
1334 350 IF (KK.GT.0) GO TO 360
1335 KODE = 2
1336 GO TO 590
1337 360 XMIN = RES(1)
1338 IOUT = S(1)
1339 J = 1
1340 IF (KK.EQ.1) GO TO 380
1341 DO 370 I=2,KK
1342 IF (RES(I).GE.XMIN) GO TO 370
1343 J = I
1344 XMIN = RES(I)
1345 IOUT = S(I)
1346 370 CONTINUE
1347 RES(J) = RES(KK)
1348 S(J) = S(KK)
1349 380 KK = KK - 1
1350 PIVOT = Q(IOUT,IN)
1351 II = Q(IOUT,N2)
1352 IF (IPHASE.EQ.1) GO TO 400
1353 IF (II.LT.0) GO TO 390
1354 IF (IU(2,II).EQ.1) GO TO 420
1355 GO TO 400
1356 390 IINEG = -II
1357 IF (IU(1,IINEG).EQ.1) GO TO 420
1358 400 II = IABS(II)
1359 CUV = CU(1,II) + CU(2,II)
1360 IF (Q(KLM1,IN)-PIVOT*CUV.LE.TOLER) GO TO 420
1361 C BYPASS INTERMEDIATE VERTICES.
1362 DO 410 J=JS,N1
1363 Z = Q(IOUT,J)
1364 Q(KLM1,J) = Q(KLM1,J) - Z*CUV
1365 Q(IOUT,J) = -Z
1366 410 CONTINUE
1367 Q(IOUT,N2) = -Q(IOUT,N2)
1368 GO TO 350
1369 C GAUSS-JORDAN ELIMINATION.
1370 420 IF (ITER.LT.MAXIT) GO TO 430
1371 KODE = 3
1372 GO TO 590
1373 430 ITER = ITER + 1
1374 DO 440 J=JS,N1
1375 IF (J.NE.IN) Q(IOUT,J) = Q(IOUT,J)/PIVOT
1376 440 CONTINUE
1377 C IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION
1378 C SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN
1379 C TO AND INCLUDING STATEMENT NUMBER 460 BY..
1380 C DO 460 J=JS,N1
1381 C IF(J .EQ. IN) GO TO 460
1382 C Z = -Q(IOUT,J)
1383 C CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1)
1384 C 460 CONTINUE
1385 DO 460 J=JS,N1
1386 IF (J.EQ.IN) GO TO 460
1387 Z = -Q(IOUT,J)
1388 DO 450 I=1,KLM1
1389 IF (I.NE.IOUT) Q(I,J) = Q(I,J) + Z*Q(I,IN)
1390 450 CONTINUE
1391 460 CONTINUE
1392 TPIVOT = -PIVOT
1393 DO 470 I=1,KLM1
1394 IF (I.NE.IOUT) Q(I,IN) = Q(I,IN)/TPIVOT
1395 470 CONTINUE
1396 Q(IOUT,IN) = 1./PIVOT
1397 Z = Q(IOUT,N2)
1398 Q(IOUT,N2) = Q(KLM2,IN)
1399 Q(KLM2,IN) = Z
1400 II = ABS(Z)
1401 IF (IU(1,II).EQ.0 .OR. IU(2,II).EQ.0) GO TO 240
1402 DO 480 I=1,KLM2
1403 Z = Q(I,IN)
1404 Q(I,IN) = Q(I,JS)
1405 Q(I,JS) = Z
1406 480 CONTINUE
1407 JS = JS + 1
1408 GO TO 240
1409 C TEST FOR OPTIMALITY.
1410 490 IF (KFORCE.EQ.0) GO TO 580
1411 IF (IPHASE.EQ.1 .AND. Q(KLM1,N1).LE.TOLER) GO TO 500
1412 KFORCE = 0
1413 GO TO 240
1414 C SET UP PHASE 2 COSTS.
1415 500 IPHASE = 2
1416 DO 510 J=1,NKLM
1417 CU(1,J) = 0.
1418 CU(2,J) = 0.
1419 510 CONTINUE
1420 DO 520 J=N1,NK
1421 CU(1,J) = 1.
1422 CU(2,J) = 1.
1423 520 CONTINUE
1424 DO 560 I=1,KLM
1425 II = Q(I,N2)
1426 IF (II.GT.0) GO TO 530
1427 II = -II
1428 IF (IU(2,II).EQ.0) GO TO 560
1429 CU(2,II) = 0.
1430 GO TO 540
1431 530 IF (IU(1,II).EQ.0) GO TO 560
1432 CU(1,II) = 0.
1433 540 IA = IA + 1
1434 DO 550 J=1,N2
1435 Z = Q(IA,J)
1436 Q(IA,J) = Q(I,J)
1437 Q(I,J) = Z
1438 550 CONTINUE
1439 560 CONTINUE
1440 GO TO 160
1441 570 IF (Q(KLM1,N1).LE.TOLER) GO TO 500
1442 KODE = 1
1443 GO TO 590
1444 580 IF (IPHASE.EQ.1) GO TO 570
1445 C PREPARE OUTPUT.
1446 KODE = 0
1447 590 SUM = 0.D0
1448 DO 600 J=1,N
1449 X(J) = 0.
1450 600 CONTINUE
1451 DO 610 I=1,KLM
1452 RES(I) = 0.
1453 610 CONTINUE
1454 DO 640 I=1,KLM
1455 II = Q(I,N2)
1456 SN = 1.
1457 IF (II.GT.0) GO TO 620
1458 II = -II
1459 SN = -1.
1460 620 IF (II.GT.N) GO TO 630
1461 X(II) = SN*Q(I,N1)
1462 GO TO 640
1463 630 IIMN = II - N
1464 RES(IIMN) = SN*Q(I,N1)
1465 IF (II.GE.N1 .AND. II.LE.NK) SUM = SUM +
1466 * DBLE(Q(I,N1))
1467 640 CONTINUE
1468 ERROR = SUM
1469 RETURN
1470 END
1471 ================================================================*/
1472