#include #include #include #include #include "f2c.h" /* prototype for function at end (from TOMS, via Netlib and f2c) */ static int 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) ; /*--------------------------------------------------------------------- Approximately (L1) solve equations j=nvec-1 z[i] = SUM A[j][i] * y[j] (i=0..ndim-1) j=0 for y[j] (j=0..nvec-1), subject to constraints (based on y[j] input) If input cony != 0, then you can supply constraints on the values of the output y[j] by putting values into the input y[j]: input y[j] = 0 ==> unconstrained = 1 ==> y[j] must be non-negative on output = -1 ==> y[j] must be non-positive on output If cony == 0, then the input y[j] is ignored. The return value of the function is E = sum |z[i]-SUM[i]| >= 0 if everything worked, and is a negative number if an error occured. -----------------------------------------------------------------------*/ float cl1_solve( int ndim, int nvec, float *z, float **A, float *y, int cony ) { /* loop counters */ int jj , ii ; /* variables for CL1 (types declared in f2c.h) */ integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ; real *q , toler , *x , *res , error , *cu ; /*-- check inputs --*/ if( ndim < 1 || nvec < 1 ) kode = 4 ; if( A == NULL || y == NULL || z == NULL ) kode = 4 ; for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ; if( kode ){ fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ; return (float)(-kode) ; } /*-- setup call to CL1 --*/ k = ndim ; l = 0 ; /* no linear equality constraints */ m = 0 ; /* no linear inequality constraints */ n = nvec ; klmd = k+l+m ; klm2d = k+l+m+2 ; nklmd = n+k+l+m ; n2d = n+2 ; kode = (cony != 0) ; /* enforce implicit constraints on x[] */ iter = 11*klmd ; toler = 0.0001 ; error = 0.0 ; /* input/output matrices & vectors */ q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ; x = (real *) calloc( 1, sizeof(real) * n2d ) ; res = (real *) calloc( 1, sizeof(real) * klmd ) ; /* workspaces */ cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ; iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ; s = (integer *) calloc( 1, sizeof(integer) * klmd ) ; /* load matrices & vectors */ for( jj=0 ; jj < nvec ; jj++ ) for( ii=0 ; ii < ndim ; ii++ ) q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */ for( ii=0 ; ii < ndim ; ii++ ) q[ii+nvec*klm2d] = z[ii] ; /* vector */ if( cony ){ for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */ x[jj] = (y[jj] < 0.0) ? -1.0 :(y[jj] > 0.0) ? 1.0 : 0.0 ; } for( ii=0 ; ii < klmd ; ii++ ) /* no constraints on resids */ res[ii] = 0.0 ; /*-- do the work --*/ cl1_fort( &k, &l, &m, &n, &klmd, &klm2d, &nklmd, &n2d, q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ; free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ; if( kode != 0 ){ free(x) ; #if 0 switch( kode ){ case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break; case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break; case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break; default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break; } #endif return (float)(-kode) ; } /*-- copy results into output --*/ for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ; free(x) ; return (float)error ; } /*--------------------------------------------------------------------- Approximately (L1) solve equations j=nvec-1 z[i] = SUM A[j][i] * y[j] (i=0..ndim-1) j=0 for y[j] (j=0..nvec-1), subject to constraints on the signs of y[j] (based on y[j] input) AND subject to constraints on the signs of the residuals z[i]-SUM[i] (based on rez[i] input). If input cony != 0, then you can supply constraints on the values of the output y[j] by putting values into the input y[j]: input y[j] = 0 ==> unconstrained = 1 ==> y[j] must be non-negative on output = -1 ==> y[j] must be non-positive on output If cony == 0, then the input y[j] is ignored. If input conr != 0, then you can supply constraints on the values of the residuals z[i]-SUM[i] by putting values into the input rez[i]: input rez[i] = 0 ==> unconstrained = 1 ==> z[i]-SUM[i] must be non-negative = -1 ==> z[i]-SUM[i] must be non-positive If conr == 0, then the input values in rez[] are ignored. The outputs are in y[] and rez[]. The return value of the function is E = sum |z[i]-SUM[i]| >= 0 if everything worked, and is a negative number if an error occured. -----------------------------------------------------------------------*/ float cl1_solve_res( int ndim, int nvec, float *z, float **A, float *y, int cony , float *rez , int conr ) { /* loop counters */ int jj , ii ; /* variables for CL1 (types declared in f2c.h) */ integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ; real *q , toler , *x , *res , error , *cu ; /*-- check inputs --*/ if( ndim < 1 || nvec < 1 ) kode = 4 ; if( A == NULL || y == NULL || z == NULL ) kode = 4 ; for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ; if( kode ){ fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ; return (float)(-kode) ; } /*-- setup call to CL1 --*/ k = ndim ; l = 0 ; /* no linear equality constraints */ m = 0 ; /* no linear inequality constraints */ n = nvec ; klmd = k+l+m ; klm2d = k+l+m+2 ; nklmd = n+k+l+m ; n2d = n+2 ; kode = (cony != 0 || conr != 0) ; /* enforce implicit constraints on x[] */ iter = 11*klmd ; toler = 0.0001 ; error = 0.0 ; /* input/output matrices & vectors */ q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ; x = (real *) calloc( 1, sizeof(real) * n2d ) ; res = (real *) calloc( 1, sizeof(real) * klmd ) ; /* workspaces */ cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ; iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ; s = (integer *) calloc( 1, sizeof(integer) * klmd ) ; /* load matrices & vectors */ for( jj=0 ; jj < nvec ; jj++ ) for( ii=0 ; ii < ndim ; ii++ ) q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */ for( ii=0 ; ii < ndim ; ii++ ) q[ii+nvec*klm2d] = z[ii] ; /* vector */ if( cony ){ for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */ x[jj] = (y[jj] < 0.0) ? -1.0 :(y[jj] > 0.0) ? 1.0 : 0.0 ; } if( conr ){ for( ii=0 ; ii < ndim ; ii++ ) /* constraints on resids */ res[ii] = (rez[ii] < 0.0) ? -1.0 :(rez[ii] > 0.0) ? 1.0 : 0.0 ; } /*-- do the work --*/ cl1_fort( &k, &l, &m, &n, &klmd, &klm2d, &nklmd, &n2d, q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ; free(q) ; free(cu) ; free(iu) ; free(s) ; if( kode != 0 ){ free(x) ; free(res) ; #if 0 switch( kode ){ case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break; case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break; case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break; default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break; } #endif return (float)(-kode) ; } /*-- copy results into output --*/ for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ; for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = (float) res[ii] ; free(res); free(x); return (float)error; } #if 0 /*--------------------------------------------------------------------- Find a set of coefficients ec[] so that j=nvec-1 { } S[i] = base_vec[i] + SUM { ec[j] * extra_vec[j][i] } j=0 { } is non-negative for i=0..ndim-1, and so that i=ndim-1 { } SUM { S[i] } i=0 { } is as small as possible. The return value of the function is 0 if everything worked, and nonzero if an error occured. If the function worked, then the coefficients are returned in ec[] (which must be allocated by the caller). Method: uses the CL1 subroutine from the TOMS library; the f2c translation appears at the end of this file N.B.: NOT TESTED -----------------------------------------------------------------------*/ int cl1_pos_sum( int ndim , int nvec , float * base_vec , float ** extra_vec , float * ec ) { /* internal variables */ int jj , ii ; /* variables for CL1 */ integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode,iter , *iu,*s ; real *q , toler , *x , *res , error , *cu ; /*-- check inputs --*/ if( ndim < 1 || nvec < 1 ) return 4 ; if( base_vec == NULL || extra_vec == NULL || ec == NULL ) return 4 ; for( jj=0 ; jj < nvec ; jj++ ) if( extra_vec[jj] == NULL ) return 4 ; /*-- setup call to CL1 --*/ k = ndim ; l = 0 ; m = 0 ; n = nvec ; klmd = k+l+m ; klm2d = k+l+m+2 ; nklmd = n+k+l+m ; n2d = n+2 ; kode = 1 ; iter = 10*klmd ; toler = 0.0001 ; error = 0.0 ; /* input/output matrices & vectors */ q = (real *) malloc( sizeof(real) * klm2d*n2d ) ; x = (real *) malloc( sizeof(real) * n2d ) ; res = (real *) malloc( sizeof(real) * klmd ) ; /* workspaces */ cu = (real *) malloc( sizeof(real) * 2*nklmd ) ; iu = (integer *) malloc( sizeof(integer) * 2*nklmd ) ; s = (integer *) malloc( sizeof(integer) * klmd ) ; /* load matrices & vectors */ for( jj=0 ; jj < nvec ; jj++ ) for( ii=0 ; ii < ndim ; ii++ ) q[ii+jj*klm2d] = extra_vec[jj][ii] ; for( ii=0 ; ii < ndim ; ii++ ) q[ii+nvec*klm2d] = base_vec[ii] ; for( jj=0 ; jj < nvec ; jj++ ) x[jj] = 0.0 ; for( ii=0 ; ii < ndim ; ii++ ) res[ii] = 1.0 ; /*-- do the work --*/ cl1_fort( &k, &l, &m, &n, &klmd, &klm2d, &nklmd, &n2d, q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ; free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ; if( kode != 0 ){ free(x) ; return (int)kode ; } /*-- copy results into output --*/ for( jj=0 ; jj < nvec ; jj++ ) ec[jj] = (float)(-x[jj]) ; free(x) ; return 0 ; } #endif /*====================================================================== Translated from the TOMS library CL1 algorithm ========================================================================*/ /* cl1.f -- translated by f2c (version 19970805). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ static int 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) { /* System generated locals */ integer q_dim1, q_offset, i__1, i__2; real r__1; /* Local variables */ static integer iimn, nklm; static real xmin, xmax; static integer iout, i__, j; static real z__; static integer iineg, maxit, n1, n2; static real pivot; static integer ia, ii, kk, in, nk, js; static real sn; static integer iphase, kforce; static real zu, zv; static integer nk1; static real tpivot; static integer klm, jmn, nkl, jpn; static real cuv; static doublereal sum; static integer klm1, klm2, nkl1; /* THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX */ /* METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION */ /* TO A K BY N SYSTEM OF LINEAR EQUATIONS */ /* AX=B */ /* SUBJECT TO L LINEAR EQUALITY CONSTRAINTS */ /* CX=D */ /* AND M LINEAR INEQUALITY CONSTRAINTS */ /* EX.LE.F. */ /* DESCRIPTION OF PARAMETERS */ /* K NUMBER OF ROWS OF THE MATRIX A (K.GE.1). */ /* L NUMBER OF ROWS OF THE MATRIX C (L.GE.0). */ /* M NUMBER OF ROWS OF THE MATRIX E (M.GE.0). */ /* N NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1). */ /* KLMD SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS. */ /* KLM2D SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS. */ /* NKLMD SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS. */ /* N2D SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS */ /* Q TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND */ /* AT LEAST N2D COLUMNS. */ /* ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS */ /* B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS */ /* AND N+1 COLUMNS OF Q AS FOLLOWS */ /* A B */ /* Q = C D */ /* E F */ /* THESE VALUES ARE DESTROYED BY THE SUBROUTINE. */ /* KODE A CODE USED ON ENTRY TO, AND EXIT */ /* FROM, THE SUBROUTINE. */ /* ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0. */ /* HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS */ /* ARE TO BE INCLUDED IMPLICITLY, RATHER THAN */ /* EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE */ /* SHOULD BE SET TO 1, AND THE NONNEGATIVITY */ /* CONSTRAINTS INCLUDED IN THE ARRAYS X AND */ /* RES (SEE BELOW). */ /* ON EXIT, KODE HAS ONE OF THE */ /* FOLLOWING VALUES */ /* 0- OPTIMAL SOLUTION FOUND, */ /* 1- NO FEASIBLE SOLUTION TO THE */ /* CONSTRAINTS, */ /* 2- CALCULATIONS TERMINATED */ /* PREMATURELY DUE TO ROUNDING ERRORS, */ /* 3- MAXIMUM NUMBER OF ITERATIONS REACHED. */ /* TOLER A SMALL POSITIVE TOLERANCE. EMPIRICAL */ /* EVIDENCE SUGGESTS TOLER = 10**(-D*2/3), */ /* WHERE D REPRESENTS THE NUMBER OF DECIMAL */ /* DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY, */ /* THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO */ /* AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED */ /* TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY */ /* NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER. */ /* ITER ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON */ /* THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. */ /* A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER */ /* GIVES THE NUMBER OF SIMPLEX ITERATIONS. */ /* X ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D. */ /* ON EXIT THIS ARRAY CONTAINS A */ /* SOLUTION TO THE L1 PROBLEM. IF KODE=1 */ /* ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE */ /* SIMPLE NONNEGATIVITY CONSTRAINTS ON THE */ /* VARIABLES. THE VALUES -1, 0, OR 1 */ /* FOR X(J) INDICATE THAT THE J-TH VARIABLE */ /* IS RESTRICTED TO BE .LE.0, UNRESTRICTED, */ /* OR .GE.0 RESPECTIVELY. */ /* RES ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD. */ /* ON EXIT THIS CONTAINS THE RESIDUALS B-AX */ /* IN THE FIRST K COMPONENTS, D-CX IN THE */ /* NEXT L COMPONENTS (THESE WILL BE =0),AND */ /* F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON */ /* ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE */ /* NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS */ /* B-AX. THE VALUES -1, 0, OR 1 FOR RES(I) */ /* INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS */ /* RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0 */ /* RESPECTIVELY. */ /* ERROR ON EXIT, THIS GIVES THE MINIMUM SUM OF */ /* ABSOLUTE VALUES OF THE RESIDUALS. */ /* CU A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND */ /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */ /* IU A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND */ /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */ /* S INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR */ /* WORKSPACE. */ /* IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO */ /* DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY */ /* THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN */ /* EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE */ /* FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS. */ /* SUBROUTINE COL(V1, V2, XMLT, NOTROW, K) */ /* THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE */ /* VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW). */ /* DIMENSION V1(K), V2(K) */ /* KEND = NOTROW - 1 */ /* KSTART = NOTROW + 1 */ /* IF (KEND .LT. 1) GO TO 20 */ /* DO 10 I=1,KEND */ /* V1(I) = V1(I) + XMLT*V2(I) */ /* 10 CONTINUE */ /* IF(KSTART .GT. K) GO TO 40 */ /* 20 DO 30 I=KSTART,K */ /* V1(I) = V1(I) + XMLT*V2(I) */ /* 30 CONTINUE */ /* 40 RETURN */ /* END */ /* SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR */ /* INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION. */ /* ----------------------------------------------------------------------- */ /* INITIALIZATION. */ /* Parameter adjustments */ --s; --res; iu -= 3; cu -= 3; --x; q_dim1 = *klm2d; q_offset = q_dim1 + 1; q -= q_offset; /* Function Body */ maxit = *iter; n1 = *n + 1; n2 = *n + 2; nk = *n + *k; nk1 = nk + 1; nkl = nk + *l; nkl1 = nkl + 1; klm = *k + *l + *m; klm1 = klm + 1; klm2 = klm + 2; nklm = *n + klm; kforce = 1; *iter = 0; js = 1; ia = 0; /* SET UP LABELS IN Q. */ i__1 = *n; for (j = 1; j <= i__1; ++j) { q[klm2 + j * q_dim1] = (real) j; /* L10: */ } i__1 = klm; for (i__ = 1; i__ <= i__1; ++i__) { q[i__ + n2 * q_dim1] = (real) (*n + i__); if (q[i__ + n1 * q_dim1] >= 0.f) { goto L30; } i__2 = n2; for (j = 1; j <= i__2; ++j) { q[i__ + j * q_dim1] = -q[i__ + j * q_dim1]; /* L20: */ } L30: ; } /* SET UP PHASE 1 COSTS. */ iphase = 2; i__1 = nklm; for (j = 1; j <= i__1; ++j) { cu[(j << 1) + 1] = 0.f; cu[(j << 1) + 2] = 0.f; iu[(j << 1) + 1] = 0; iu[(j << 1) + 2] = 0; /* L40: */ } if (*l == 0) { goto L60; } i__1 = nkl; for (j = nk1; j <= i__1; ++j) { cu[(j << 1) + 1] = 1.f; cu[(j << 1) + 2] = 1.f; iu[(j << 1) + 1] = 1; iu[(j << 1) + 2] = 1; /* L50: */ } iphase = 1; L60: if (*m == 0) { goto L80; } i__1 = nklm; for (j = nkl1; j <= i__1; ++j) { cu[(j << 1) + 2] = 1.f; iu[(j << 1) + 2] = 1; jmn = j - *n; if (q[jmn + n2 * q_dim1] < 0.f) { iphase = 1; } /* L70: */ } L80: if (*kode == 0) { goto L150; } i__1 = *n; for (j = 1; j <= i__1; ++j) { if ((r__1 = x[j]) < 0.f) { goto L90; } else if (r__1 == 0) { goto L110; } else { goto L100; } L90: cu[(j << 1) + 1] = 1.f; iu[(j << 1) + 1] = 1; goto L110; L100: cu[(j << 1) + 2] = 1.f; iu[(j << 1) + 2] = 1; L110: ; } i__1 = *k; for (j = 1; j <= i__1; ++j) { jpn = j + *n; if ((r__1 = res[j]) < 0.f) { goto L120; } else if (r__1 == 0) { goto L140; } else { goto L130; } L120: cu[(jpn << 1) + 1] = 1.f; iu[(jpn << 1) + 1] = 1; if (q[j + n2 * q_dim1] > 0.f) { iphase = 1; } goto L140; L130: cu[(jpn << 1) + 2] = 1.f; iu[(jpn << 1) + 2] = 1; if (q[j + n2 * q_dim1] < 0.f) { iphase = 1; } L140: ; } L150: if (iphase == 2) { goto L500; } /* COMPUTE THE MARGINAL COSTS. */ L160: i__1 = n1; for (j = js; j <= i__1; ++j) { sum = 0.; i__2 = klm; for (i__ = 1; i__ <= i__2; ++i__) { ii = q[i__ + n2 * q_dim1]; if (ii < 0) { goto L170; } z__ = cu[(ii << 1) + 1]; goto L180; L170: iineg = -ii; z__ = cu[(iineg << 1) + 2]; L180: sum += (doublereal) q[i__ + j * q_dim1] * (doublereal) z__; /* L190: */ } q[klm1 + j * q_dim1] = sum; /* L200: */ } i__1 = *n; for (j = js; j <= i__1; ++j) { ii = q[klm2 + j * q_dim1]; if (ii < 0) { goto L210; } z__ = cu[(ii << 1) + 1]; goto L220; L210: iineg = -ii; z__ = cu[(iineg << 1) + 2]; L220: q[klm1 + j * q_dim1] -= z__; /* L230: */ } /* DETERMINE THE VECTOR TO ENTER THE BASIS. */ L240: xmax = 0.f; if (js > *n) { goto L490; } i__1 = *n; for (j = js; j <= i__1; ++j) { zu = q[klm1 + j * q_dim1]; ii = q[klm2 + j * q_dim1]; if (ii > 0) { goto L250; } ii = -ii; zv = zu; zu = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2]; goto L260; L250: zv = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2]; L260: if (kforce == 1 && ii > *n) { goto L280; } if (iu[(ii << 1) + 1] == 1) { goto L270; } if (zu <= xmax) { goto L270; } xmax = zu; in = j; L270: if (iu[(ii << 1) + 2] == 1) { goto L280; } if (zv <= xmax) { goto L280; } xmax = zv; in = j; L280: ; } if (xmax <= *toler) { goto L490; } if (q[klm1 + in * q_dim1] == xmax) { goto L300; } i__1 = klm2; for (i__ = 1; i__ <= i__1; ++i__) { q[i__ + in * q_dim1] = -q[i__ + in * q_dim1]; /* L290: */ } q[klm1 + in * q_dim1] = xmax; /* DETERMINE THE VECTOR TO LEAVE THE BASIS. */ L300: if (iphase == 1 || ia == 0) { goto L330; } xmax = 0.f; i__1 = ia; for (i__ = 1; i__ <= i__1; ++i__) { z__ = (r__1 = q[i__ + in * q_dim1], dabs(r__1)); if (z__ <= xmax) { goto L310; } xmax = z__; iout = i__; L310: ; } if (xmax <= *toler) { goto L330; } i__1 = n2; for (j = 1; j <= i__1; ++j) { z__ = q[ia + j * q_dim1]; q[ia + j * q_dim1] = q[iout + j * q_dim1]; q[iout + j * q_dim1] = z__; /* L320: */ } iout = ia; --ia; pivot = q[iout + in * q_dim1]; goto L420; L330: kk = 0; i__1 = klm; for (i__ = 1; i__ <= i__1; ++i__) { z__ = q[i__ + in * q_dim1]; if (z__ <= *toler) { goto L340; } ++kk; res[kk] = q[i__ + n1 * q_dim1] / z__; s[kk] = i__; L340: ; } L350: if (kk > 0) { goto L360; } *kode = 2; goto L590; L360: xmin = res[1]; iout = s[1]; j = 1; if (kk == 1) { goto L380; } i__1 = kk; for (i__ = 2; i__ <= i__1; ++i__) { if (res[i__] >= xmin) { goto L370; } j = i__; xmin = res[i__]; iout = s[i__]; L370: ; } res[j] = res[kk]; s[j] = s[kk]; L380: --kk; pivot = q[iout + in * q_dim1]; ii = q[iout + n2 * q_dim1]; if (iphase == 1) { goto L400; } if (ii < 0) { goto L390; } if (iu[(ii << 1) + 2] == 1) { goto L420; } goto L400; L390: iineg = -ii; if (iu[(iineg << 1) + 1] == 1) { goto L420; } L400: ii = abs(ii); cuv = cu[(ii << 1) + 1] + cu[(ii << 1) + 2]; if (q[klm1 + in * q_dim1] - pivot * cuv <= *toler) { goto L420; } /* BYPASS INTERMEDIATE VERTICES. */ i__1 = n1; for (j = js; j <= i__1; ++j) { z__ = q[iout + j * q_dim1]; q[klm1 + j * q_dim1] -= z__ * cuv; q[iout + j * q_dim1] = -z__; /* L410: */ } q[iout + n2 * q_dim1] = -q[iout + n2 * q_dim1]; goto L350; /* GAUSS-JORDAN ELIMINATION. */ L420: if (*iter < maxit) { goto L430; } *kode = 3; goto L590; L430: ++(*iter); i__1 = n1; for (j = js; j <= i__1; ++j) { if (j != in) { q[iout + j * q_dim1] /= pivot; } /* L440: */ } /* IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION */ /* SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN */ /* TO AND INCLUDING STATEMENT NUMBER 460 BY.. */ /* DO 460 J=JS,N1 */ /* IF(J .EQ. IN) GO TO 460 */ /* Z = -Q(IOUT,J) */ /* CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1) */ /* 460 CONTINUE */ i__1 = n1; for (j = js; j <= i__1; ++j) { if (j == in) { goto L460; } z__ = -q[iout + j * q_dim1]; i__2 = klm1; for (i__ = 1; i__ <= i__2; ++i__) { if (i__ != iout) { q[i__ + j * q_dim1] += z__ * q[i__ + in * q_dim1]; } /* L450: */ } L460: ; } tpivot = -pivot; i__1 = klm1; for (i__ = 1; i__ <= i__1; ++i__) { if (i__ != iout) { q[i__ + in * q_dim1] /= tpivot; } /* L470: */ } q[iout + in * q_dim1] = 1.f / pivot; z__ = q[iout + n2 * q_dim1]; q[iout + n2 * q_dim1] = q[klm2 + in * q_dim1]; q[klm2 + in * q_dim1] = z__; ii = dabs(z__); if (iu[(ii << 1) + 1] == 0 || iu[(ii << 1) + 2] == 0) { goto L240; } i__1 = klm2; for (i__ = 1; i__ <= i__1; ++i__) { z__ = q[i__ + in * q_dim1]; q[i__ + in * q_dim1] = q[i__ + js * q_dim1]; q[i__ + js * q_dim1] = z__; /* L480: */ } ++js; goto L240; /* TEST FOR OPTIMALITY. */ L490: if (kforce == 0) { goto L580; } if (iphase == 1 && q[klm1 + n1 * q_dim1] <= *toler) { goto L500; } kforce = 0; goto L240; /* SET UP PHASE 2 COSTS. */ L500: iphase = 2; i__1 = nklm; for (j = 1; j <= i__1; ++j) { cu[(j << 1) + 1] = 0.f; cu[(j << 1) + 2] = 0.f; /* L510: */ } i__1 = nk; for (j = n1; j <= i__1; ++j) { cu[(j << 1) + 1] = 1.f; cu[(j << 1) + 2] = 1.f; /* L520: */ } i__1 = klm; for (i__ = 1; i__ <= i__1; ++i__) { ii = q[i__ + n2 * q_dim1]; if (ii > 0) { goto L530; } ii = -ii; if (iu[(ii << 1) + 2] == 0) { goto L560; } cu[(ii << 1) + 2] = 0.f; goto L540; L530: if (iu[(ii << 1) + 1] == 0) { goto L560; } cu[(ii << 1) + 1] = 0.f; L540: ++ia; i__2 = n2; for (j = 1; j <= i__2; ++j) { z__ = q[ia + j * q_dim1]; q[ia + j * q_dim1] = q[i__ + j * q_dim1]; q[i__ + j * q_dim1] = z__; /* L550: */ } L560: ; } goto L160; L570: if (q[klm1 + n1 * q_dim1] <= *toler) { goto L500; } *kode = 1; goto L590; L580: if (iphase == 1) { goto L570; } /* PREPARE OUTPUT. */ *kode = 0; L590: sum = 0.; i__1 = *n; for (j = 1; j <= i__1; ++j) { x[j] = 0.f; /* L600: */ } i__1 = klm; for (i__ = 1; i__ <= i__1; ++i__) { res[i__] = 0.f; /* L610: */ } i__1 = klm; for (i__ = 1; i__ <= i__1; ++i__) { ii = q[i__ + n2 * q_dim1]; sn = 1.f; if (ii > 0) { goto L620; } ii = -ii; sn = -1.f; L620: if (ii > *n) { goto L630; } x[ii] = sn * q[i__ + n1 * q_dim1]; goto L640; L630: iimn = ii - *n; res[iimn] = sn * q[i__ + n1 * q_dim1]; if (ii >= n1 && ii <= nk) { sum += (doublereal) q[i__ + n1 * q_dim1]; } L640: ; } *error = sum; return 0; } /* cl1_ */ /*=================================================================== SUBROUTINE CL1(K, L, M, N, KLMD, KLM2D, NKLMD, N2D, * Q, KODE, TOLER, ITER, X, RES, ERROR, CU, IU, S) C C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX C METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION C TO A K BY N SYSTEM OF LINEAR EQUATIONS C AX=B C SUBJECT TO L LINEAR EQUALITY CONSTRAINTS C CX=D C AND M LINEAR INEQUALITY CONSTRAINTS C EX.LE.F. C C DESCRIPTION OF PARAMETERS C C K NUMBER OF ROWS OF THE MATRIX A (K.GE.1). C L NUMBER OF ROWS OF THE MATRIX C (L.GE.0). C M NUMBER OF ROWS OF THE MATRIX E (M.GE.0). C N NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1). C KLMD SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS. C KLM2D SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS. C NKLMD SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS. C N2D SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS C C Q TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND C AT LEAST N2D COLUMNS. C ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS C B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS C AND N+1 COLUMNS OF Q AS FOLLOWS C A B C Q = C D C E F C THESE VALUES ARE DESTROYED BY THE SUBROUTINE. C C KODE A CODE USED ON ENTRY TO, AND EXIT C FROM, THE SUBROUTINE. C ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0. C HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS C ARE TO BE INCLUDED IMPLICITLY, RATHER THAN C EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE C SHOULD BE SET TO 1, AND THE NONNEGATIVITY C CONSTRAINTS INCLUDED IN THE ARRAYS X AND C RES (SEE BELOW). C ON EXIT, KODE HAS ONE OF THE C FOLLOWING VALUES C 0- OPTIMAL SOLUTION FOUND, C 1- NO FEASIBLE SOLUTION TO THE C CONSTRAINTS, C 2- CALCULATIONS TERMINATED C PREMATURELY DUE TO ROUNDING ERRORS, C 3- MAXIMUM NUMBER OF ITERATIONS REACHED. C C TOLER A SMALL POSITIVE TOLERANCE. EMPIRICAL C EVIDENCE SUGGESTS TOLER = 10**(-D*2/3), C WHERE D REPRESENTS THE NUMBER OF DECIMAL C DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY, C THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO C AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED C TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY C NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER. C C ITER ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON C THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER C GIVES THE NUMBER OF SIMPLEX ITERATIONS. C C X ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D. C ON EXIT THIS ARRAY CONTAINS A C SOLUTION TO THE L1 PROBLEM. IF KODE=1 C ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE C SIMPLE NONNEGATIVITY CONSTRAINTS ON THE C VARIABLES. THE VALUES -1, 0, OR 1 C FOR X(J) INDICATE THAT THE J-TH VARIABLE C IS RESTRICTED TO BE .LE.0, UNRESTRICTED, C OR .GE.0 RESPECTIVELY. C C RES ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD. C ON EXIT THIS CONTAINS THE RESIDUALS B-AX C IN THE FIRST K COMPONENTS, D-CX IN THE C NEXT L COMPONENTS (THESE WILL BE =0),AND C F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON C ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE C NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS C B-AX. THE VALUES -1, 0, OR 1 FOR RES(I) C INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS C RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0 C RESPECTIVELY. C C ERROR ON EXIT, THIS GIVES THE MINIMUM SUM OF C ABSOLUTE VALUES OF THE RESIDUALS. C C CU A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND C AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. C C IU A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND C AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. C C S INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR C WORKSPACE. C C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN C EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE C FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS. C SUBROUTINE COL(V1, V2, XMLT, NOTROW, K) C THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE C VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW). C DIMENSION V1(K), V2(K) C KEND = NOTROW - 1 C KSTART = NOTROW + 1 C IF (KEND .LT. 1) GO TO 20 C DO 10 I=1,KEND C V1(I) = V1(I) + XMLT*V2(I) C 10 CONTINUE C IF(KSTART .GT. K) GO TO 40 C 20 DO 30 I=KSTART,K C V1(I) = V1(I) + XMLT*V2(I) C 30 CONTINUE C 40 RETURN C END C SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR C INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION. C----------------------------------------------------------------------- DOUBLE PRECISION SUM DOUBLE PRECISION DBLE REAL Q, X, Z, CU, SN, ZU, ZV, CUV, RES, XMAX, XMIN, * ERROR, PIVOT, TOLER, TPIVOT REAL ABS INTEGER I, J, K, L, M, N, S, IA, II, IN, IU, JS, KK, * NK, N1, N2, JMN, JPN, KLM, NKL, NK1, N2D, IIMN, * IOUT, ITER, KLMD, KLM1, KLM2, KODE, NKLM, NKL1, * KLM2D, MAXIT, NKLMD, IPHASE, KFORCE, IINEG INTEGER IABS DIMENSION Q(KLM2D,N2D), X(N2D), RES(KLMD), * CU(2,NKLMD), IU(2,NKLMD), S(KLMD) C C INITIALIZATION. C MAXIT = ITER N1 = N + 1 N2 = N + 2 NK = N + K NK1 = NK + 1 NKL = NK + L NKL1 = NKL + 1 KLM = K + L + M KLM1 = KLM + 1 KLM2 = KLM + 2 NKLM = N + KLM KFORCE = 1 ITER = 0 JS = 1 IA = 0 C SET UP LABELS IN Q. DO 10 J=1,N Q(KLM2,J) = J 10 CONTINUE DO 30 I=1,KLM Q(I,N2) = N + I IF (Q(I,N1).GE.0.) GO TO 30 DO 20 J=1,N2 Q(I,J) = -Q(I,J) 20 CONTINUE 30 CONTINUE C SET UP PHASE 1 COSTS. IPHASE = 2 DO 40 J=1,NKLM CU(1,J) = 0. CU(2,J) = 0. IU(1,J) = 0 IU(2,J) = 0 40 CONTINUE IF (L.EQ.0) GO TO 60 DO 50 J=NK1,NKL CU(1,J) = 1. CU(2,J) = 1. IU(1,J) = 1 IU(2,J) = 1 50 CONTINUE IPHASE = 1 60 IF (M.EQ.0) GO TO 80 DO 70 J=NKL1,NKLM CU(2,J) = 1. IU(2,J) = 1 JMN = J - N IF (Q(JMN,N2).LT.0.) IPHASE = 1 70 CONTINUE 80 IF (KODE.EQ.0) GO TO 150 DO 110 J=1,N IF (X(J)) 90, 110, 100 90 CU(1,J) = 1. IU(1,J) = 1 GO TO 110 100 CU(2,J) = 1. IU(2,J) = 1 110 CONTINUE DO 140 J=1,K JPN = J + N IF (RES(J)) 120, 140, 130 120 CU(1,JPN) = 1. IU(1,JPN) = 1 IF (Q(J,N2).GT.0.0) IPHASE = 1 GO TO 140 130 CU(2,JPN) = 1. IU(2,JPN) = 1 IF (Q(J,N2).LT.0.0) IPHASE = 1 140 CONTINUE 150 IF (IPHASE.EQ.2) GO TO 500 C COMPUTE THE MARGINAL COSTS. 160 DO 200 J=JS,N1 SUM = 0.D0 DO 190 I=1,KLM II = Q(I,N2) IF (II.LT.0) GO TO 170 Z = CU(1,II) GO TO 180 170 IINEG = -II Z = CU(2,IINEG) 180 SUM = SUM + DBLE(Q(I,J))*DBLE(Z) 190 CONTINUE Q(KLM1,J) = SUM 200 CONTINUE DO 230 J=JS,N II = Q(KLM2,J) IF (II.LT.0) GO TO 210 Z = CU(1,II) GO TO 220 210 IINEG = -II Z = CU(2,IINEG) 220 Q(KLM1,J) = Q(KLM1,J) - Z 230 CONTINUE C DETERMINE THE VECTOR TO ENTER THE BASIS. 240 XMAX = 0. IF (JS.GT.N) GO TO 490 DO 280 J=JS,N ZU = Q(KLM1,J) II = Q(KLM2,J) IF (II.GT.0) GO TO 250 II = -II ZV = ZU ZU = -ZU - CU(1,II) - CU(2,II) GO TO 260 250 ZV = -ZU - CU(1,II) - CU(2,II) 260 IF (KFORCE.EQ.1 .AND. II.GT.N) GO TO 280 IF (IU(1,II).EQ.1) GO TO 270 IF (ZU.LE.XMAX) GO TO 270 XMAX = ZU IN = J 270 IF (IU(2,II).EQ.1) GO TO 280 IF (ZV.LE.XMAX) GO TO 280 XMAX = ZV IN = J 280 CONTINUE IF (XMAX.LE.TOLER) GO TO 490 IF (Q(KLM1,IN).EQ.XMAX) GO TO 300 DO 290 I=1,KLM2 Q(I,IN) = -Q(I,IN) 290 CONTINUE Q(KLM1,IN) = XMAX C DETERMINE THE VECTOR TO LEAVE THE BASIS. 300 IF (IPHASE.EQ.1 .OR. IA.EQ.0) GO TO 330 XMAX = 0. DO 310 I=1,IA Z = ABS(Q(I,IN)) IF (Z.LE.XMAX) GO TO 310 XMAX = Z IOUT = I 310 CONTINUE IF (XMAX.LE.TOLER) GO TO 330 DO 320 J=1,N2 Z = Q(IA,J) Q(IA,J) = Q(IOUT,J) Q(IOUT,J) = Z 320 CONTINUE IOUT = IA IA = IA - 1 PIVOT = Q(IOUT,IN) GO TO 420 330 KK = 0 DO 340 I=1,KLM Z = Q(I,IN) IF (Z.LE.TOLER) GO TO 340 KK = KK + 1 RES(KK) = Q(I,N1)/Z S(KK) = I 340 CONTINUE 350 IF (KK.GT.0) GO TO 360 KODE = 2 GO TO 590 360 XMIN = RES(1) IOUT = S(1) J = 1 IF (KK.EQ.1) GO TO 380 DO 370 I=2,KK IF (RES(I).GE.XMIN) GO TO 370 J = I XMIN = RES(I) IOUT = S(I) 370 CONTINUE RES(J) = RES(KK) S(J) = S(KK) 380 KK = KK - 1 PIVOT = Q(IOUT,IN) II = Q(IOUT,N2) IF (IPHASE.EQ.1) GO TO 400 IF (II.LT.0) GO TO 390 IF (IU(2,II).EQ.1) GO TO 420 GO TO 400 390 IINEG = -II IF (IU(1,IINEG).EQ.1) GO TO 420 400 II = IABS(II) CUV = CU(1,II) + CU(2,II) IF (Q(KLM1,IN)-PIVOT*CUV.LE.TOLER) GO TO 420 C BYPASS INTERMEDIATE VERTICES. DO 410 J=JS,N1 Z = Q(IOUT,J) Q(KLM1,J) = Q(KLM1,J) - Z*CUV Q(IOUT,J) = -Z 410 CONTINUE Q(IOUT,N2) = -Q(IOUT,N2) GO TO 350 C GAUSS-JORDAN ELIMINATION. 420 IF (ITER.LT.MAXIT) GO TO 430 KODE = 3 GO TO 590 430 ITER = ITER + 1 DO 440 J=JS,N1 IF (J.NE.IN) Q(IOUT,J) = Q(IOUT,J)/PIVOT 440 CONTINUE C IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION C SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN C TO AND INCLUDING STATEMENT NUMBER 460 BY.. C DO 460 J=JS,N1 C IF(J .EQ. IN) GO TO 460 C Z = -Q(IOUT,J) C CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1) C 460 CONTINUE DO 460 J=JS,N1 IF (J.EQ.IN) GO TO 460 Z = -Q(IOUT,J) DO 450 I=1,KLM1 IF (I.NE.IOUT) Q(I,J) = Q(I,J) + Z*Q(I,IN) 450 CONTINUE 460 CONTINUE TPIVOT = -PIVOT DO 470 I=1,KLM1 IF (I.NE.IOUT) Q(I,IN) = Q(I,IN)/TPIVOT 470 CONTINUE Q(IOUT,IN) = 1./PIVOT Z = Q(IOUT,N2) Q(IOUT,N2) = Q(KLM2,IN) Q(KLM2,IN) = Z II = ABS(Z) IF (IU(1,II).EQ.0 .OR. IU(2,II).EQ.0) GO TO 240 DO 480 I=1,KLM2 Z = Q(I,IN) Q(I,IN) = Q(I,JS) Q(I,JS) = Z 480 CONTINUE JS = JS + 1 GO TO 240 C TEST FOR OPTIMALITY. 490 IF (KFORCE.EQ.0) GO TO 580 IF (IPHASE.EQ.1 .AND. Q(KLM1,N1).LE.TOLER) GO TO 500 KFORCE = 0 GO TO 240 C SET UP PHASE 2 COSTS. 500 IPHASE = 2 DO 510 J=1,NKLM CU(1,J) = 0. CU(2,J) = 0. 510 CONTINUE DO 520 J=N1,NK CU(1,J) = 1. CU(2,J) = 1. 520 CONTINUE DO 560 I=1,KLM II = Q(I,N2) IF (II.GT.0) GO TO 530 II = -II IF (IU(2,II).EQ.0) GO TO 560 CU(2,II) = 0. GO TO 540 530 IF (IU(1,II).EQ.0) GO TO 560 CU(1,II) = 0. 540 IA = IA + 1 DO 550 J=1,N2 Z = Q(IA,J) Q(IA,J) = Q(I,J) Q(I,J) = Z 550 CONTINUE 560 CONTINUE GO TO 160 570 IF (Q(KLM1,N1).LE.TOLER) GO TO 500 KODE = 1 GO TO 590 580 IF (IPHASE.EQ.1) GO TO 570 C PREPARE OUTPUT. KODE = 0 590 SUM = 0.D0 DO 600 J=1,N X(J) = 0. 600 CONTINUE DO 610 I=1,KLM RES(I) = 0. 610 CONTINUE DO 640 I=1,KLM II = Q(I,N2) SN = 1. IF (II.GT.0) GO TO 620 II = -II SN = -1. 620 IF (II.GT.N) GO TO 630 X(II) = SN*Q(I,N1) GO TO 640 630 IIMN = II - N RES(IIMN) = SN*Q(I,N1) IF (II.GE.N1 .AND. II.LE.NK) SUM = SUM + * DBLE(Q(I,N1)) 640 CONTINUE ERROR = SUM RETURN END ================================================================*/