00001 #include <math.h>
00002 #include <stddef.h>
00003 #include <stdio.h>
00004 #include <stdlib.h>
00005 #include "f2c.h"
00006 
00007 /* prototype for function at end (from TOMS, via Netlib and f2c) */
00008 
00009 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
00010                     integer *klmd, integer *klm2d, integer *nklmd,
00011                     integer *n2d, real *q,
00012                     integer *kode, real *toler, integer *iter, real *x,
00013                     real *res, real * error, real *cu, integer *iu, integer *s) ;
00014 
00015 /*---------------------------------------------------------------------
00016   Approximately (L1) solve equations
00017 
00018              j=nvec-1
00019     z[i] = SUM        A[j][i] * y[j]   (i=0..ndim-1)
00020              j=0
00021 
00022   for y[j] (j=0..nvec-1), subject to constraints (based on y[j] input)
00023 
00024   If input cony != 0, then you can supply constraints on the values
00025   of the output y[j] by putting values into the input y[j]:
00026 
00027   input y[j] =  0 ==> unconstrained
00028              =  1 ==> y[j] must be non-negative on output
00029              = -1 ==> y[j] must be non-positive on output
00030 
00031   If cony == 0, then the input y[j] is ignored.
00032 
00033   The return value of the function is E = sum |z[i]-SUM[i]| >= 0
00034   if everything worked, and is a negative number if an error occured.
00035 -----------------------------------------------------------------------*/
00036 
00037 float cl1_solve( int ndim, int nvec, float *z, float **A, float *y, int cony )
00038 {
00039    /* loop counters */
00040 
00041    int jj , ii ;
00042 
00043    /* variables for CL1 (types declared in f2c.h) */
00044 
00045    integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00046    real *q , toler , *x , *res , error , *cu ;
00047 
00048    /*-- check inputs --*/
00049 
00050    if( ndim < 1 || nvec < 1 )                         kode = 4 ;
00051    if( A == NULL || y == NULL || z == NULL )          kode = 4 ;
00052    for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00053 
00054    if( kode ){
00055      fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00056      return (float)(-kode) ;
00057    }
00058 
00059    /*-- setup call to CL1 --*/
00060 
00061    k     = ndim ;
00062    l     = 0 ;     /* no linear equality constraints */
00063    m     = 0 ;     /* no linear inequality constraints */
00064    n     = nvec ;
00065 
00066    klmd  = k+l+m ;
00067    klm2d = k+l+m+2 ;
00068    nklmd = n+k+l+m ;
00069    n2d   = n+2 ;
00070 
00071    kode  = (cony != 0) ; /* enforce implicit constraints on x[] */
00072    iter  = 11*klmd ;
00073 
00074    toler = 0.0001 ;
00075    error = 0.0 ;
00076 
00077    /* input/output matrices & vectors */
00078 
00079    q     = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00080    x     = (real *) calloc( 1, sizeof(real) * n2d ) ;
00081    res   = (real *) calloc( 1, sizeof(real) * klmd ) ;
00082 
00083    /* workspaces */
00084 
00085    cu    = (real *)    calloc( 1, sizeof(real) * 2*nklmd ) ;
00086    iu    = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00087    s     = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00088 
00089    /* load matrices & vectors */
00090 
00091    for( jj=0 ; jj < nvec ; jj++ )
00092       for( ii=0 ; ii < ndim ; ii++ )
00093          q[ii+jj*klm2d] = A[jj][ii] ;   /* matrix */
00094 
00095    for( ii=0 ; ii < ndim ; ii++ )
00096       q[ii+nvec*klm2d] = z[ii] ;        /* vector */
00097 
00098    if( cony ){
00099      for( jj=0 ; jj < nvec ; jj++ )     /* constraints on solution */
00100        x[jj] = (y[jj] < 0.0) ? -1.0
00101               :(y[jj] > 0.0) ?  1.0 : 0.0 ;
00102    }
00103 
00104    for( ii=0 ; ii < klmd ; ii++ )       /* no constraints on resids */
00105       res[ii] = 0.0 ;
00106 
00107    /*-- do the work --*/
00108 
00109    cl1_fort( &k, &l, &m, &n,
00110              &klmd, &klm2d, &nklmd, &n2d,
00111              q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00112 
00113    free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
00114 
00115    if( kode != 0 ){
00116      free(x) ;
00117 #if 0
00118      switch( kode ){
00119        case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00120        case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n")     ; break;
00121        case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n")      ; break;
00122       default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n")     ; break;
00123      }
00124 #endif
00125      return (float)(-kode) ;
00126    }
00127 
00128    /*-- copy results into output --*/
00129 
00130    for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00131 
00132    free(x) ; return (float)error ;
00133 }
00134 
00135 /*---------------------------------------------------------------------
00136   Approximately (L1) solve equations
00137 
00138              j=nvec-1
00139     z[i] = SUM        A[j][i] * y[j]   (i=0..ndim-1)
00140              j=0
00141 
00142   for y[j] (j=0..nvec-1), subject to constraints on the signs of y[j]
00143   (based on y[j] input) AND subject to constraints on the signs of
00144   the residuals z[i]-SUM[i] (based on rez[i] input).
00145 
00146   If input cony != 0, then you can supply constraints on the values
00147   of the output y[j] by putting values into the input y[j]:
00148 
00149   input y[j] =  0 ==> unconstrained
00150              =  1 ==> y[j] must be non-negative on output
00151              = -1 ==> y[j] must be non-positive on output
00152 
00153   If cony == 0, then the input y[j] is ignored.
00154 
00155   If input conr != 0, then you can supply constraints on the values
00156   of the residuals z[i]-SUM[i] by putting values into the input rez[i]:
00157 
00158   input rez[i] =  0 ==> unconstrained
00159                =  1 ==> z[i]-SUM[i] must be non-negative
00160                = -1 ==> z[i]-SUM[i] must be non-positive
00161 
00162   If conr == 0, then the input values in rez[] are ignored.
00163 
00164   The outputs are in y[] and rez[].
00165 
00166   The return value of the function is E = sum |z[i]-SUM[i]| >= 0
00167   if everything worked, and is a negative number if an error occured.
00168 -----------------------------------------------------------------------*/
00169 
00170 float cl1_solve_res( int ndim, int nvec, float *z, float **A,
00171                      float *y, int cony , float *rez , int conr )
00172 {
00173    /* loop counters */
00174 
00175    int jj , ii ;
00176 
00177    /* variables for CL1 (types declared in f2c.h) */
00178 
00179    integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00180    real *q , toler , *x , *res , error , *cu ;
00181 
00182    /*-- check inputs --*/
00183 
00184    if( ndim < 1 || nvec < 1 )                         kode = 4 ;
00185    if( A == NULL || y == NULL || z == NULL )          kode = 4 ;
00186    for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00187 
00188    if( kode ){
00189      fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00190      return (float)(-kode) ;
00191    }
00192 
00193    /*-- setup call to CL1 --*/
00194 
00195    k     = ndim ;
00196    l     = 0 ;     /* no linear equality constraints */
00197    m     = 0 ;     /* no linear inequality constraints */
00198    n     = nvec ;
00199 
00200    klmd  = k+l+m ;
00201    klm2d = k+l+m+2 ;
00202    nklmd = n+k+l+m ;
00203    n2d   = n+2 ;
00204 
00205    kode  = (cony != 0 || conr != 0) ; /* enforce implicit constraints on x[] */
00206    iter  = 11*klmd ;
00207 
00208    toler = 0.0001 ;
00209    error = 0.0 ;
00210 
00211    /* input/output matrices & vectors */
00212 
00213    q     = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00214    x     = (real *) calloc( 1, sizeof(real) * n2d ) ;
00215    res   = (real *) calloc( 1, sizeof(real) * klmd ) ;
00216 
00217    /* workspaces */
00218 
00219    cu    = (real *)    calloc( 1, sizeof(real) * 2*nklmd ) ;
00220    iu    = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00221    s     = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00222 
00223    /* load matrices & vectors */
00224 
00225    for( jj=0 ; jj < nvec ; jj++ )
00226       for( ii=0 ; ii < ndim ; ii++ )
00227          q[ii+jj*klm2d] = A[jj][ii] ;   /* matrix */
00228 
00229    for( ii=0 ; ii < ndim ; ii++ )
00230       q[ii+nvec*klm2d] = z[ii] ;        /* vector */
00231 
00232    if( cony ){
00233      for( jj=0 ; jj < nvec ; jj++ )     /* constraints on solution */
00234        x[jj] = (y[jj] < 0.0) ? -1.0
00235               :(y[jj] > 0.0) ?  1.0 : 0.0 ;
00236    }
00237 
00238    if( conr ){
00239      for( ii=0 ; ii < ndim ; ii++ )     /* constraints on resids */
00240        res[ii] = (rez[ii] < 0.0) ? -1.0
00241                 :(rez[ii] > 0.0) ?  1.0 : 0.0 ;
00242    }
00243 
00244    /*-- do the work --*/
00245 
00246    cl1_fort( &k, &l, &m, &n,
00247              &klmd, &klm2d, &nklmd, &n2d,
00248              q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00249 
00250    free(q) ; free(cu) ; free(iu) ; free(s) ;
00251 
00252    if( kode != 0 ){
00253      free(x) ; free(res) ;
00254 #if 0
00255      switch( kode ){
00256        case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00257        case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n")     ; break;
00258        case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n")      ; break;
00259       default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n")     ; break;
00260      }
00261 #endif
00262      return (float)(-kode) ;
00263    }
00264 
00265    /*-- copy results into output --*/
00266 
00267    for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00268 
00269    for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = (float) res[ii] ;
00270 
00271    free(res); free(x); return (float)error;
00272 }
00273 
00274 #if 0
00275 /*---------------------------------------------------------------------
00276   Find a set of coefficients ec[] so that
00277 
00278                               j=nvec-1 {                         }
00279       S[i] =  base_vec[i] + SUM        { ec[j] * extra_vec[j][i] }
00280                               j=0      {                         }
00281 
00282   is non-negative for i=0..ndim-1, and so that
00283 
00284         i=ndim-1 {      }
00285       SUM        { S[i] }
00286         i=0      {      }
00287 
00288   is as small as possible.
00289 
00290   The return value of the function is 0 if everything worked, and
00291   nonzero if an error occured.  If the function worked, then the
00292   coefficients are returned in ec[] (which must be allocated by
00293   the caller).
00294 
00295   Method: uses the CL1 subroutine from the TOMS library;
00296           the f2c translation appears at the end of this file
00297 
00298   N.B.: NOT TESTED
00299 -----------------------------------------------------------------------*/
00300 
00301 int cl1_pos_sum( int ndim , int nvec ,
00302                  float * base_vec , float ** extra_vec , float * ec )
00303 {
00304    /* internal variables */
00305 
00306    int jj , ii ;
00307 
00308    /* variables for CL1 */
00309 
00310    integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode,iter , *iu,*s ;
00311    real *q , toler , *x , *res , error , *cu ;
00312 
00313    /*-- check inputs --*/
00314 
00315    if( ndim < 1 || nvec < 1 )                                 return 4 ;
00316    if( base_vec == NULL || extra_vec == NULL || ec == NULL )  return 4 ;
00317    for( jj=0 ; jj < nvec ; jj++ ) if( extra_vec[jj] == NULL ) return 4 ;
00318 
00319    /*-- setup call to CL1 --*/
00320 
00321    k     = ndim ;
00322    l     = 0 ;
00323    m     = 0 ;
00324    n     = nvec ;
00325 
00326    klmd  = k+l+m ;
00327    klm2d = k+l+m+2 ;
00328    nklmd = n+k+l+m ;
00329    n2d   = n+2 ;
00330 
00331    kode  = 1 ;
00332    iter  = 10*klmd ;
00333 
00334    toler = 0.0001 ;
00335    error = 0.0 ;
00336 
00337    /* input/output matrices & vectors */
00338 
00339    q     = (real *) malloc( sizeof(real) * klm2d*n2d ) ;
00340    x     = (real *) malloc( sizeof(real) * n2d ) ;
00341    res   = (real *) malloc( sizeof(real) * klmd ) ;
00342 
00343    /* workspaces */
00344 
00345    cu    = (real *) malloc( sizeof(real) * 2*nklmd ) ;
00346    iu    = (integer *) malloc( sizeof(integer) * 2*nklmd ) ;
00347    s     = (integer *) malloc( sizeof(integer) * klmd ) ;
00348 
00349    /* load matrices & vectors */
00350 
00351    for( jj=0 ; jj < nvec ; jj++ )
00352       for( ii=0 ; ii < ndim ; ii++ )
00353          q[ii+jj*klm2d] = extra_vec[jj][ii] ;
00354 
00355    for( ii=0 ; ii < ndim ; ii++ )
00356       q[ii+nvec*klm2d] = base_vec[ii] ;
00357 
00358    for( jj=0 ; jj < nvec ; jj++ ) x[jj] = 0.0 ;
00359 
00360    for( ii=0 ; ii < ndim ; ii++ ) res[ii] = 1.0 ;
00361 
00362    /*-- do the work --*/
00363 
00364    cl1_fort( &k, &l, &m, &n,
00365              &klmd, &klm2d, &nklmd, &n2d,
00366              q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00367 
00368    free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
00369 
00370    if( kode != 0 ){ free(x) ; return (int)kode ; }
00371 
00372    /*-- copy results into output --*/
00373 
00374    for( jj=0 ; jj < nvec ; jj++ ) ec[jj] = (float)(-x[jj]) ;
00375 
00376    free(x) ; return 0 ;
00377 }
00378 #endif
00379 
00380 /*======================================================================
00381    Translated from the TOMS library CL1 algorithm
00382 ========================================================================*/
00383 
00384 /* cl1.f -- translated by f2c (version 19970805).
00385    You must link the resulting object file with the libraries:
00386         -lf2c -lm   (in that order)
00387 */
00388 
00389 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
00390         integer *klmd, integer *klm2d, integer *nklmd, integer *n2d, real *q,
00391         integer *kode, real *toler, integer *iter, real *x, real *res, real *
00392         error, real *cu, integer *iu, integer *s)
00393 {
00394     /* System generated locals */
00395     integer q_dim1, q_offset, i__1, i__2;
00396     real r__1;
00397 
00398     /* Local variables */
00399     static integer iimn, nklm;
00400     static real xmin, xmax;
00401     static integer iout, i__, j;
00402     static real z__;
00403     static integer iineg, maxit, n1, n2;
00404     static real pivot;
00405     static integer ia, ii, kk, in, nk, js;
00406     static real sn;
00407     static integer iphase, kforce;
00408     static real zu, zv;
00409     static integer nk1;
00410     static real tpivot;
00411     static integer klm, jmn, nkl, jpn;
00412     static real cuv;
00413     static doublereal sum;
00414     static integer klm1, klm2, nkl1;
00415 
00416 
00417 /* THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX */
00418 /* METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION */
00419 /* TO A K BY N SYSTEM OF LINEAR EQUATIONS */
00420 /*             AX=B */
00421 /* SUBJECT TO L LINEAR EQUALITY CONSTRAINTS */
00422 /*             CX=D */
00423 /* AND M LINEAR INEQUALITY CONSTRAINTS */
00424 /*             EX.LE.F. */
00425 
00426 /* DESCRIPTION OF PARAMETERS */
00427 
00428 /* K      NUMBER OF ROWS OF THE MATRIX A (K.GE.1). */
00429 /* L      NUMBER OF ROWS OF THE MATRIX C (L.GE.0). */
00430 /* M      NUMBER OF ROWS OF THE MATRIX E (M.GE.0). */
00431 /* N      NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1). */
00432 /* KLMD   SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS. */
00433 /* KLM2D  SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS. */
00434 /* NKLMD  SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS. */
00435 /* N2D    SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS */
00436 
00437 /* Q      TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND */
00438 /*        AT LEAST N2D COLUMNS. */
00439 /*        ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS */
00440 /*        B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS */
00441 /*        AND N+1 COLUMNS OF Q AS FOLLOWS */
00442 /*             A B */
00443 /*         Q = C D */
00444 /*             E F */
00445 /*        THESE VALUES ARE DESTROYED BY THE SUBROUTINE. */
00446 
00447 /* KODE   A CODE USED ON ENTRY TO, AND EXIT */
00448 /*        FROM, THE SUBROUTINE. */
00449 /*        ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0. */
00450 /*        HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS */
00451 /*        ARE TO BE INCLUDED IMPLICITLY, RATHER THAN */
00452 /*        EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE */
00453 /*        SHOULD BE SET TO 1, AND THE NONNEGATIVITY */
00454 /*        CONSTRAINTS INCLUDED IN THE ARRAYS X AND */
00455 /*        RES (SEE BELOW). */
00456 /*        ON EXIT, KODE HAS ONE OF THE */
00457 /*        FOLLOWING VALUES */
00458 /*             0- OPTIMAL SOLUTION FOUND, */
00459 /*             1- NO FEASIBLE SOLUTION TO THE */
00460 /*                CONSTRAINTS, */
00461 /*             2- CALCULATIONS TERMINATED */
00462 /*                PREMATURELY DUE TO ROUNDING ERRORS, */
00463 /*             3- MAXIMUM NUMBER OF ITERATIONS REACHED. */
00464 
00465 /* TOLER  A SMALL POSITIVE TOLERANCE. EMPIRICAL */
00466 /*        EVIDENCE SUGGESTS TOLER = 10**(-D*2/3), */
00467 /*        WHERE D REPRESENTS THE NUMBER OF DECIMAL */
00468 /*        DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY, */
00469 /*        THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO */
00470 /*        AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED */
00471 /*        TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY */
00472 /*        NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER. */
00473 
00474 /* ITER   ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON */
00475 /*        THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. */
00476 /*        A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER */
00477 /*        GIVES THE NUMBER OF SIMPLEX ITERATIONS. */
00478 
00479 /* X      ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D. */
00480 /*        ON EXIT THIS ARRAY CONTAINS A */
00481 /*        SOLUTION TO THE L1 PROBLEM. IF KODE=1 */
00482 /*        ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE */
00483 /*        SIMPLE NONNEGATIVITY CONSTRAINTS ON THE */
00484 /*        VARIABLES. THE VALUES -1, 0, OR 1 */
00485 /*        FOR X(J) INDICATE THAT THE J-TH VARIABLE */
00486 /*        IS RESTRICTED TO BE .LE.0, UNRESTRICTED, */
00487 /*        OR .GE.0 RESPECTIVELY. */
00488 
00489 /* RES    ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD. */
00490 /*        ON EXIT THIS CONTAINS THE RESIDUALS B-AX */
00491 /*        IN THE FIRST K COMPONENTS, D-CX IN THE */
00492 /*        NEXT L COMPONENTS (THESE WILL BE =0),AND */
00493 /*        F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON */
00494 /*        ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE */
00495 /*        NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS */
00496 /*        B-AX. THE VALUES -1, 0, OR 1 FOR RES(I) */
00497 /*        INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS */
00498 /*        RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0 */
00499 /*        RESPECTIVELY. */
00500 
00501 /* ERROR  ON EXIT, THIS GIVES THE MINIMUM SUM OF */
00502 /*        ABSOLUTE VALUES OF THE RESIDUALS. */
00503 
00504 /* CU     A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND */
00505 /*        AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
00506 
00507 /* IU     A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND */
00508 /*        AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
00509 
00510 /* S      INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR */
00511 /*        WORKSPACE. */
00512 
00513 /* IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO */
00514 /* DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY */
00515 /* THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN */
00516 /* EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE */
00517 /* FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS. */
00518 /*     SUBROUTINE COL(V1, V2, XMLT, NOTROW, K) */
00519 /* THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE */
00520 /* VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW). */
00521 /*     DIMENSION V1(K), V2(K) */
00522 /*     KEND = NOTROW - 1 */
00523 /*     KSTART = NOTROW + 1 */
00524 /*     IF (KEND .LT. 1) GO TO 20 */
00525 /*     DO 10 I=1,KEND */
00526 /*        V1(I) = V1(I) + XMLT*V2(I) */
00527 /*  10 CONTINUE */
00528 /*     IF(KSTART .GT. K) GO TO 40 */
00529 /*  20 DO 30 I=KSTART,K */
00530 /*       V1(I) = V1(I) + XMLT*V2(I) */
00531 /*  30 CONTINUE */
00532 /*  40 RETURN */
00533 /*     END */
00534 /* SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR */
00535 /* INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION. */
00536 /* -----------------------------------------------------------------------
00537  */
00538 
00539 /* INITIALIZATION. */
00540 
00541     /* Parameter adjustments */
00542     --s;
00543     --res;
00544     iu -= 3;
00545     cu -= 3;
00546     --x;
00547     q_dim1 = *klm2d;
00548     q_offset = q_dim1 + 1;
00549     q -= q_offset;
00550 
00551     /* Function Body */
00552     maxit = *iter;
00553     n1 = *n + 1;
00554     n2 = *n + 2;
00555     nk = *n + *k;
00556     nk1 = nk + 1;
00557     nkl = nk + *l;
00558     nkl1 = nkl + 1;
00559     klm = *k + *l + *m;
00560     klm1 = klm + 1;
00561     klm2 = klm + 2;
00562     nklm = *n + klm;
00563     kforce = 1;
00564     *iter = 0;
00565     js = 1;
00566     ia = 0;
00567 /* SET UP LABELS IN Q. */
00568     i__1 = *n;
00569     for (j = 1; j <= i__1; ++j) {
00570         q[klm2 + j * q_dim1] = (real) j;
00571 /* L10: */
00572     }
00573     i__1 = klm;
00574     for (i__ = 1; i__ <= i__1; ++i__) {
00575         q[i__ + n2 * q_dim1] = (real) (*n + i__);
00576         if (q[i__ + n1 * q_dim1] >= 0.f) {
00577             goto L30;
00578         }
00579         i__2 = n2;
00580         for (j = 1; j <= i__2; ++j) {
00581             q[i__ + j * q_dim1] = -q[i__ + j * q_dim1];
00582 /* L20: */
00583         }
00584 L30:
00585         ;
00586     }
00587 /* SET UP PHASE 1 COSTS. */
00588     iphase = 2;
00589     i__1 = nklm;
00590     for (j = 1; j <= i__1; ++j) {
00591         cu[(j << 1) + 1] = 0.f;
00592         cu[(j << 1) + 2] = 0.f;
00593         iu[(j << 1) + 1] = 0;
00594         iu[(j << 1) + 2] = 0;
00595 /* L40: */
00596     }
00597     if (*l == 0) {
00598         goto L60;
00599     }
00600     i__1 = nkl;
00601     for (j = nk1; j <= i__1; ++j) {
00602         cu[(j << 1) + 1] = 1.f;
00603         cu[(j << 1) + 2] = 1.f;
00604         iu[(j << 1) + 1] = 1;
00605         iu[(j << 1) + 2] = 1;
00606 /* L50: */
00607     }
00608     iphase = 1;
00609 L60:
00610     if (*m == 0) {
00611         goto L80;
00612     }
00613     i__1 = nklm;
00614     for (j = nkl1; j <= i__1; ++j) {
00615         cu[(j << 1) + 2] = 1.f;
00616         iu[(j << 1) + 2] = 1;
00617         jmn = j - *n;
00618         if (q[jmn + n2 * q_dim1] < 0.f) {
00619             iphase = 1;
00620         }
00621 /* L70: */
00622     }
00623 L80:
00624     if (*kode == 0) {
00625         goto L150;
00626     }
00627     i__1 = *n;
00628     for (j = 1; j <= i__1; ++j) {
00629         if ((r__1 = x[j]) < 0.f) {
00630             goto L90;
00631         } else if (r__1 == 0) {
00632             goto L110;
00633         } else {
00634             goto L100;
00635         }
00636 L90:
00637         cu[(j << 1) + 1] = 1.f;
00638         iu[(j << 1) + 1] = 1;
00639         goto L110;
00640 L100:
00641         cu[(j << 1) + 2] = 1.f;
00642         iu[(j << 1) + 2] = 1;
00643 L110:
00644         ;
00645     }
00646     i__1 = *k;
00647     for (j = 1; j <= i__1; ++j) {
00648         jpn = j + *n;
00649         if ((r__1 = res[j]) < 0.f) {
00650             goto L120;
00651         } else if (r__1 == 0) {
00652             goto L140;
00653         } else {
00654             goto L130;
00655         }
00656 L120:
00657         cu[(jpn << 1) + 1] = 1.f;
00658         iu[(jpn << 1) + 1] = 1;
00659         if (q[j + n2 * q_dim1] > 0.f) {
00660             iphase = 1;
00661         }
00662         goto L140;
00663 L130:
00664         cu[(jpn << 1) + 2] = 1.f;
00665         iu[(jpn << 1) + 2] = 1;
00666         if (q[j + n2 * q_dim1] < 0.f) {
00667             iphase = 1;
00668         }
00669 L140:
00670         ;
00671     }
00672 L150:
00673     if (iphase == 2) {
00674         goto L500;
00675     }
00676 /* COMPUTE THE MARGINAL COSTS. */
00677 L160:
00678     i__1 = n1;
00679     for (j = js; j <= i__1; ++j) {
00680         sum = 0.;
00681         i__2 = klm;
00682         for (i__ = 1; i__ <= i__2; ++i__) {
00683             ii = q[i__ + n2 * q_dim1];
00684             if (ii < 0) {
00685                 goto L170;
00686             }
00687             z__ = cu[(ii << 1) + 1];
00688             goto L180;
00689 L170:
00690             iineg = -ii;
00691             z__ = cu[(iineg << 1) + 2];
00692 L180:
00693             sum += (doublereal) q[i__ + j * q_dim1] * (doublereal) z__;
00694 /* L190: */
00695         }
00696         q[klm1 + j * q_dim1] = sum;
00697 /* L200: */
00698     }
00699     i__1 = *n;
00700     for (j = js; j <= i__1; ++j) {
00701         ii = q[klm2 + j * q_dim1];
00702         if (ii < 0) {
00703             goto L210;
00704         }
00705         z__ = cu[(ii << 1) + 1];
00706         goto L220;
00707 L210:
00708         iineg = -ii;
00709         z__ = cu[(iineg << 1) + 2];
00710 L220:
00711         q[klm1 + j * q_dim1] -= z__;
00712 /* L230: */
00713     }
00714 /* DETERMINE THE VECTOR TO ENTER THE BASIS. */
00715 L240:
00716     xmax = 0.f;
00717     if (js > *n) {
00718         goto L490;
00719     }
00720     i__1 = *n;
00721     for (j = js; j <= i__1; ++j) {
00722         zu = q[klm1 + j * q_dim1];
00723         ii = q[klm2 + j * q_dim1];
00724         if (ii > 0) {
00725             goto L250;
00726         }
00727         ii = -ii;
00728         zv = zu;
00729         zu = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
00730         goto L260;
00731 L250:
00732         zv = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
00733 L260:
00734         if (kforce == 1 && ii > *n) {
00735             goto L280;
00736         }
00737         if (iu[(ii << 1) + 1] == 1) {
00738             goto L270;
00739         }
00740         if (zu <= xmax) {
00741             goto L270;
00742         }
00743         xmax = zu;
00744         in = j;
00745 L270:
00746         if (iu[(ii << 1) + 2] == 1) {
00747             goto L280;
00748         }
00749         if (zv <= xmax) {
00750             goto L280;
00751         }
00752         xmax = zv;
00753         in = j;
00754 L280:
00755         ;
00756     }
00757     if (xmax <= *toler) {
00758         goto L490;
00759     }
00760     if (q[klm1 + in * q_dim1] == xmax) {
00761         goto L300;
00762     }
00763     i__1 = klm2;
00764     for (i__ = 1; i__ <= i__1; ++i__) {
00765         q[i__ + in * q_dim1] = -q[i__ + in * q_dim1];
00766 /* L290: */
00767     }
00768     q[klm1 + in * q_dim1] = xmax;
00769 /* DETERMINE THE VECTOR TO LEAVE THE BASIS. */
00770 L300:
00771     if (iphase == 1 || ia == 0) {
00772         goto L330;
00773     }
00774     xmax = 0.f;
00775     i__1 = ia;
00776     for (i__ = 1; i__ <= i__1; ++i__) {
00777         z__ = (r__1 = q[i__ + in * q_dim1], dabs(r__1));
00778         if (z__ <= xmax) {
00779             goto L310;
00780         }
00781         xmax = z__;
00782         iout = i__;
00783 L310:
00784         ;
00785     }
00786     if (xmax <= *toler) {
00787         goto L330;
00788     }
00789     i__1 = n2;
00790     for (j = 1; j <= i__1; ++j) {
00791         z__ = q[ia + j * q_dim1];
00792         q[ia + j * q_dim1] = q[iout + j * q_dim1];
00793         q[iout + j * q_dim1] = z__;
00794 /* L320: */
00795     }
00796     iout = ia;
00797     --ia;
00798     pivot = q[iout + in * q_dim1];
00799     goto L420;
00800 L330:
00801     kk = 0;
00802     i__1 = klm;
00803     for (i__ = 1; i__ <= i__1; ++i__) {
00804         z__ = q[i__ + in * q_dim1];
00805         if (z__ <= *toler) {
00806             goto L340;
00807         }
00808         ++kk;
00809         res[kk] = q[i__ + n1 * q_dim1] / z__;
00810         s[kk] = i__;
00811 L340:
00812         ;
00813     }
00814 L350:
00815     if (kk > 0) {
00816         goto L360;
00817     }
00818     *kode = 2;
00819     goto L590;
00820 L360:
00821     xmin = res[1];
00822     iout = s[1];
00823     j = 1;
00824     if (kk == 1) {
00825         goto L380;
00826     }
00827     i__1 = kk;
00828     for (i__ = 2; i__ <= i__1; ++i__) {
00829         if (res[i__] >= xmin) {
00830             goto L370;
00831         }
00832         j = i__;
00833         xmin = res[i__];
00834         iout = s[i__];
00835 L370:
00836         ;
00837     }
00838     res[j] = res[kk];
00839     s[j] = s[kk];
00840 L380:
00841     --kk;
00842     pivot = q[iout + in * q_dim1];
00843     ii = q[iout + n2 * q_dim1];
00844     if (iphase == 1) {
00845         goto L400;
00846     }
00847     if (ii < 0) {
00848         goto L390;
00849     }
00850     if (iu[(ii << 1) + 2] == 1) {
00851         goto L420;
00852     }
00853     goto L400;
00854 L390:
00855     iineg = -ii;
00856     if (iu[(iineg << 1) + 1] == 1) {
00857         goto L420;
00858     }
00859 L400:
00860     ii = abs(ii);
00861     cuv = cu[(ii << 1) + 1] + cu[(ii << 1) + 2];
00862     if (q[klm1 + in * q_dim1] - pivot * cuv <= *toler) {
00863         goto L420;
00864     }
00865 /* BYPASS INTERMEDIATE VERTICES. */
00866     i__1 = n1;
00867     for (j = js; j <= i__1; ++j) {
00868         z__ = q[iout + j * q_dim1];
00869         q[klm1 + j * q_dim1] -= z__ * cuv;
00870         q[iout + j * q_dim1] = -z__;
00871 /* L410: */
00872     }
00873     q[iout + n2 * q_dim1] = -q[iout + n2 * q_dim1];
00874     goto L350;
00875 /* GAUSS-JORDAN ELIMINATION. */
00876 L420:
00877     if (*iter < maxit) {
00878         goto L430;
00879     }
00880     *kode = 3;
00881     goto L590;
00882 L430:
00883     ++(*iter);
00884     i__1 = n1;
00885     for (j = js; j <= i__1; ++j) {
00886         if (j != in) {
00887             q[iout + j * q_dim1] /= pivot;
00888         }
00889 /* L440: */
00890     }
00891 /* IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION */
00892 /* SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN */
00893 /* TO AND INCLUDING STATEMENT NUMBER 460 BY.. */
00894 /*     DO 460 J=JS,N1 */
00895 /*        IF(J .EQ. IN) GO TO 460 */
00896 /*        Z = -Q(IOUT,J) */
00897 /*        CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1) */
00898 /* 460 CONTINUE */
00899     i__1 = n1;
00900     for (j = js; j <= i__1; ++j) {
00901         if (j == in) {
00902             goto L460;
00903         }
00904         z__ = -q[iout + j * q_dim1];
00905         i__2 = klm1;
00906         for (i__ = 1; i__ <= i__2; ++i__) {
00907             if (i__ != iout) {
00908                 q[i__ + j * q_dim1] += z__ * q[i__ + in * q_dim1];
00909             }
00910 /* L450: */
00911         }
00912 L460:
00913         ;
00914     }
00915     tpivot = -pivot;
00916     i__1 = klm1;
00917     for (i__ = 1; i__ <= i__1; ++i__) {
00918         if (i__ != iout) {
00919             q[i__ + in * q_dim1] /= tpivot;
00920         }
00921 /* L470: */
00922     }
00923     q[iout + in * q_dim1] = 1.f / pivot;
00924     z__ = q[iout + n2 * q_dim1];
00925     q[iout + n2 * q_dim1] = q[klm2 + in * q_dim1];
00926     q[klm2 + in * q_dim1] = z__;
00927     ii = dabs(z__);
00928     if (iu[(ii << 1) + 1] == 0 || iu[(ii << 1) + 2] == 0) {
00929         goto L240;
00930     }
00931     i__1 = klm2;
00932     for (i__ = 1; i__ <= i__1; ++i__) {
00933         z__ = q[i__ + in * q_dim1];
00934         q[i__ + in * q_dim1] = q[i__ + js * q_dim1];
00935         q[i__ + js * q_dim1] = z__;
00936 /* L480: */
00937     }
00938     ++js;
00939     goto L240;
00940 /* TEST FOR OPTIMALITY. */
00941 L490:
00942     if (kforce == 0) {
00943         goto L580;
00944     }
00945     if (iphase == 1 && q[klm1 + n1 * q_dim1] <= *toler) {
00946         goto L500;
00947     }
00948     kforce = 0;
00949     goto L240;
00950 /* SET UP PHASE 2 COSTS. */
00951 L500:
00952     iphase = 2;
00953     i__1 = nklm;
00954     for (j = 1; j <= i__1; ++j) {
00955         cu[(j << 1) + 1] = 0.f;
00956         cu[(j << 1) + 2] = 0.f;
00957 /* L510: */
00958     }
00959     i__1 = nk;
00960     for (j = n1; j <= i__1; ++j) {
00961         cu[(j << 1) + 1] = 1.f;
00962         cu[(j << 1) + 2] = 1.f;
00963 /* L520: */
00964     }
00965     i__1 = klm;
00966     for (i__ = 1; i__ <= i__1; ++i__) {
00967         ii = q[i__ + n2 * q_dim1];
00968         if (ii > 0) {
00969             goto L530;
00970         }
00971         ii = -ii;
00972         if (iu[(ii << 1) + 2] == 0) {
00973             goto L560;
00974         }
00975         cu[(ii << 1) + 2] = 0.f;
00976         goto L540;
00977 L530:
00978         if (iu[(ii << 1) + 1] == 0) {
00979             goto L560;
00980         }
00981         cu[(ii << 1) + 1] = 0.f;
00982 L540:
00983         ++ia;
00984         i__2 = n2;
00985         for (j = 1; j <= i__2; ++j) {
00986             z__ = q[ia + j * q_dim1];
00987             q[ia + j * q_dim1] = q[i__ + j * q_dim1];
00988             q[i__ + j * q_dim1] = z__;
00989 /* L550: */
00990         }
00991 L560:
00992         ;
00993     }
00994     goto L160;
00995 L570:
00996     if (q[klm1 + n1 * q_dim1] <= *toler) {
00997         goto L500;
00998     }
00999     *kode = 1;
01000     goto L590;
01001 L580:
01002     if (iphase == 1) {
01003         goto L570;
01004     }
01005 /* PREPARE OUTPUT. */
01006     *kode = 0;
01007 L590:
01008     sum = 0.;
01009     i__1 = *n;
01010     for (j = 1; j <= i__1; ++j) {
01011         x[j] = 0.f;
01012 /* L600: */
01013     }
01014     i__1 = klm;
01015     for (i__ = 1; i__ <= i__1; ++i__) {
01016         res[i__] = 0.f;
01017 /* L610: */
01018     }
01019     i__1 = klm;
01020     for (i__ = 1; i__ <= i__1; ++i__) {
01021         ii = q[i__ + n2 * q_dim1];
01022         sn = 1.f;
01023         if (ii > 0) {
01024             goto L620;
01025         }
01026         ii = -ii;
01027         sn = -1.f;
01028 L620:
01029         if (ii > *n) {
01030             goto L630;
01031         }
01032         x[ii] = sn * q[i__ + n1 * q_dim1];
01033         goto L640;
01034 L630:
01035         iimn = ii - *n;
01036         res[iimn] = sn * q[i__ + n1 * q_dim1];
01037         if (ii >= n1 && ii <= nk) {
01038             sum += (doublereal) q[i__ + n1 * q_dim1];
01039         }
01040 L640:
01041         ;
01042     }
01043     *error = sum;
01044     return 0;
01045 } /* cl1_ */
01046 
01047 /*===================================================================
01048       SUBROUTINE CL1(K, L, M, N, KLMD, KLM2D, NKLMD, N2D,
01049      * Q, KODE, TOLER, ITER, X, RES, ERROR, CU, IU, S)
01050 C
01051 C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX
01052 C METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION
01053 C TO A K BY N SYSTEM OF LINEAR EQUATIONS
01054 C             AX=B
01055 C SUBJECT TO L LINEAR EQUALITY CONSTRAINTS
01056 C             CX=D
01057 C AND M LINEAR INEQUALITY CONSTRAINTS
01058 C             EX.LE.F.
01059 C
01060 C DESCRIPTION OF PARAMETERS
01061 C
01062 C K      NUMBER OF ROWS OF THE MATRIX A (K.GE.1).
01063 C L      NUMBER OF ROWS OF THE MATRIX C (L.GE.0).
01064 C M      NUMBER OF ROWS OF THE MATRIX E (M.GE.0).
01065 C N      NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1).
01066 C KLMD   SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS.
01067 C KLM2D  SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS.
01068 C NKLMD  SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS.
01069 C N2D    SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS
01070 C
01071 C Q      TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND
01072 C        AT LEAST N2D COLUMNS.
01073 C        ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS
01074 C        B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS
01075 C        AND N+1 COLUMNS OF Q AS FOLLOWS
01076 C             A B
01077 C         Q = C D
01078 C             E F
01079 C        THESE VALUES ARE DESTROYED BY THE SUBROUTINE.
01080 C
01081 C KODE   A CODE USED ON ENTRY TO, AND EXIT
01082 C        FROM, THE SUBROUTINE.
01083 C        ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0.
01084 C        HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS
01085 C        ARE TO BE INCLUDED IMPLICITLY, RATHER THAN
01086 C        EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE
01087 C        SHOULD BE SET TO 1, AND THE NONNEGATIVITY
01088 C        CONSTRAINTS INCLUDED IN THE ARRAYS X AND
01089 C        RES (SEE BELOW).
01090 C        ON EXIT, KODE HAS ONE OF THE
01091 C        FOLLOWING VALUES
01092 C             0- OPTIMAL SOLUTION FOUND,
01093 C             1- NO FEASIBLE SOLUTION TO THE
01094 C                CONSTRAINTS,
01095 C             2- CALCULATIONS TERMINATED
01096 C                PREMATURELY DUE TO ROUNDING ERRORS,
01097 C             3- MAXIMUM NUMBER OF ITERATIONS REACHED.
01098 C
01099 C TOLER  A SMALL POSITIVE TOLERANCE. EMPIRICAL
01100 C        EVIDENCE SUGGESTS TOLER = 10**(-D*2/3),
01101 C        WHERE D REPRESENTS THE NUMBER OF DECIMAL
01102 C        DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY,
01103 C        THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO
01104 C        AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED
01105 C        TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY
01106 C        NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER.
01107 C
01108 C ITER   ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON
01109 C        THE MAXIMUM NUMBER OF ITERATIONS ALLOWED.
01110 C        A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER
01111 C        GIVES THE NUMBER OF SIMPLEX ITERATIONS.
01112 C
01113 C X      ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D.
01114 C        ON EXIT THIS ARRAY CONTAINS A
01115 C        SOLUTION TO THE L1 PROBLEM. IF KODE=1
01116 C        ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE
01117 C        SIMPLE NONNEGATIVITY CONSTRAINTS ON THE
01118 C        VARIABLES. THE VALUES -1, 0, OR 1
01119 C        FOR X(J) INDICATE THAT THE J-TH VARIABLE
01120 C        IS RESTRICTED TO BE .LE.0, UNRESTRICTED,
01121 C        OR .GE.0 RESPECTIVELY.
01122 C
01123 C RES    ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD.
01124 C        ON EXIT THIS CONTAINS THE RESIDUALS B-AX
01125 C        IN THE FIRST K COMPONENTS, D-CX IN THE
01126 C        NEXT L COMPONENTS (THESE WILL BE =0),AND
01127 C        F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON
01128 C        ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE
01129 C        NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS
01130 C        B-AX. THE VALUES -1, 0, OR 1 FOR RES(I)
01131 C        INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS
01132 C        RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0
01133 C        RESPECTIVELY.
01134 C
01135 C ERROR  ON EXIT, THIS GIVES THE MINIMUM SUM OF
01136 C        ABSOLUTE VALUES OF THE RESIDUALS.
01137 C
01138 C CU     A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND
01139 C        AT LEAST NKLMD COLUMNS USED FOR WORKSPACE.
01140 C
01141 C IU     A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND
01142 C        AT LEAST NKLMD COLUMNS USED FOR WORKSPACE.
01143 C
01144 C S      INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR
01145 C        WORKSPACE.
01146 C
01147 C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO
01148 C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY
01149 C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN
01150 C EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE
01151 C FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS.
01152 C     SUBROUTINE COL(V1, V2, XMLT, NOTROW, K)
01153 C THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE
01154 C VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW).
01155 C     DIMENSION V1(K), V2(K)
01156 C     KEND = NOTROW - 1
01157 C     KSTART = NOTROW + 1
01158 C     IF (KEND .LT. 1) GO TO 20
01159 C     DO 10 I=1,KEND
01160 C        V1(I) = V1(I) + XMLT*V2(I)
01161 C  10 CONTINUE
01162 C     IF(KSTART .GT. K) GO TO 40
01163 C  20 DO 30 I=KSTART,K
01164 C       V1(I) = V1(I) + XMLT*V2(I)
01165 C  30 CONTINUE
01166 C  40 RETURN
01167 C     END
01168 C SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR
01169 C INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION.
01170 C-----------------------------------------------------------------------
01171       DOUBLE PRECISION SUM
01172       DOUBLE PRECISION DBLE
01173       REAL Q, X, Z, CU, SN, ZU, ZV, CUV, RES, XMAX, XMIN,
01174      * ERROR, PIVOT, TOLER, TPIVOT
01175       REAL ABS
01176       INTEGER I, J, K, L, M, N, S, IA, II, IN, IU, JS, KK,
01177      * NK, N1, N2, JMN, JPN, KLM, NKL, NK1, N2D, IIMN,
01178      * IOUT, ITER, KLMD, KLM1, KLM2, KODE, NKLM, NKL1,
01179      * KLM2D, MAXIT, NKLMD, IPHASE, KFORCE, IINEG
01180       INTEGER IABS
01181       DIMENSION Q(KLM2D,N2D), X(N2D), RES(KLMD),
01182      * CU(2,NKLMD), IU(2,NKLMD), S(KLMD)
01183 C
01184 C INITIALIZATION.
01185 C
01186       MAXIT = ITER
01187       N1 = N + 1
01188       N2 = N + 2
01189       NK = N + K
01190       NK1 = NK + 1
01191       NKL = NK + L
01192       NKL1 = NKL + 1
01193       KLM = K + L + M
01194       KLM1 = KLM + 1
01195       KLM2 = KLM + 2
01196       NKLM = N + KLM
01197       KFORCE = 1
01198       ITER = 0
01199       JS = 1
01200       IA = 0
01201 C SET UP LABELS IN Q.
01202       DO 10 J=1,N
01203          Q(KLM2,J) = J
01204    10 CONTINUE
01205       DO 30 I=1,KLM
01206          Q(I,N2) = N + I
01207          IF (Q(I,N1).GE.0.) GO TO 30
01208          DO 20 J=1,N2
01209             Q(I,J) = -Q(I,J)
01210    20    CONTINUE
01211    30 CONTINUE
01212 C SET UP PHASE 1 COSTS.
01213       IPHASE = 2
01214       DO 40 J=1,NKLM
01215          CU(1,J) = 0.
01216          CU(2,J) = 0.
01217          IU(1,J) = 0
01218          IU(2,J) = 0
01219    40 CONTINUE
01220       IF (L.EQ.0) GO TO 60
01221       DO 50 J=NK1,NKL
01222          CU(1,J) = 1.
01223          CU(2,J) = 1.
01224          IU(1,J) = 1
01225          IU(2,J) = 1
01226    50 CONTINUE
01227       IPHASE = 1
01228    60 IF (M.EQ.0) GO TO 80
01229       DO 70 J=NKL1,NKLM
01230          CU(2,J) = 1.
01231          IU(2,J) = 1
01232          JMN = J - N
01233          IF (Q(JMN,N2).LT.0.) IPHASE = 1
01234    70 CONTINUE
01235    80 IF (KODE.EQ.0) GO TO 150
01236       DO 110 J=1,N
01237          IF (X(J)) 90, 110, 100
01238    90    CU(1,J) = 1.
01239          IU(1,J) = 1
01240          GO TO 110
01241   100    CU(2,J) = 1.
01242          IU(2,J) = 1
01243   110 CONTINUE
01244       DO 140 J=1,K
01245          JPN = J + N
01246          IF (RES(J)) 120, 140, 130
01247   120    CU(1,JPN) = 1.
01248          IU(1,JPN) = 1
01249          IF (Q(J,N2).GT.0.0) IPHASE = 1
01250          GO TO 140
01251   130    CU(2,JPN) = 1.
01252          IU(2,JPN) = 1
01253          IF (Q(J,N2).LT.0.0) IPHASE = 1
01254   140 CONTINUE
01255   150 IF (IPHASE.EQ.2) GO TO 500
01256 C COMPUTE THE MARGINAL COSTS.
01257   160 DO 200 J=JS,N1
01258          SUM = 0.D0
01259          DO 190 I=1,KLM
01260             II = Q(I,N2)
01261             IF (II.LT.0) GO TO 170
01262             Z = CU(1,II)
01263             GO TO 180
01264   170       IINEG = -II
01265             Z = CU(2,IINEG)
01266   180       SUM = SUM + DBLE(Q(I,J))*DBLE(Z)
01267   190    CONTINUE
01268          Q(KLM1,J) = SUM
01269   200 CONTINUE
01270       DO 230 J=JS,N
01271          II = Q(KLM2,J)
01272          IF (II.LT.0) GO TO 210
01273          Z = CU(1,II)
01274          GO TO 220
01275   210    IINEG = -II
01276          Z = CU(2,IINEG)
01277   220    Q(KLM1,J) = Q(KLM1,J) - Z
01278   230 CONTINUE
01279 C DETERMINE THE VECTOR TO ENTER THE BASIS.
01280   240 XMAX = 0.
01281       IF (JS.GT.N) GO TO 490
01282       DO 280 J=JS,N
01283          ZU = Q(KLM1,J)
01284          II = Q(KLM2,J)
01285          IF (II.GT.0) GO TO 250
01286          II = -II
01287          ZV = ZU
01288          ZU = -ZU - CU(1,II) - CU(2,II)
01289          GO TO 260
01290   250    ZV = -ZU - CU(1,II) - CU(2,II)
01291   260    IF (KFORCE.EQ.1 .AND. II.GT.N) GO TO 280
01292          IF (IU(1,II).EQ.1) GO TO 270
01293          IF (ZU.LE.XMAX) GO TO 270
01294          XMAX = ZU
01295          IN = J
01296   270    IF (IU(2,II).EQ.1) GO TO 280
01297          IF (ZV.LE.XMAX) GO TO 280
01298          XMAX = ZV
01299          IN = J
01300   280 CONTINUE
01301       IF (XMAX.LE.TOLER) GO TO 490
01302       IF (Q(KLM1,IN).EQ.XMAX) GO TO 300
01303       DO 290 I=1,KLM2
01304          Q(I,IN) = -Q(I,IN)
01305   290 CONTINUE
01306       Q(KLM1,IN) = XMAX
01307 C DETERMINE THE VECTOR TO LEAVE THE BASIS.
01308   300 IF (IPHASE.EQ.1 .OR. IA.EQ.0) GO TO 330
01309       XMAX = 0.
01310       DO 310 I=1,IA
01311          Z = ABS(Q(I,IN))
01312          IF (Z.LE.XMAX) GO TO 310
01313          XMAX = Z
01314          IOUT = I
01315   310 CONTINUE
01316       IF (XMAX.LE.TOLER) GO TO 330
01317       DO 320 J=1,N2
01318          Z = Q(IA,J)
01319          Q(IA,J) = Q(IOUT,J)
01320          Q(IOUT,J) = Z
01321   320 CONTINUE
01322       IOUT = IA
01323       IA = IA - 1
01324       PIVOT = Q(IOUT,IN)
01325       GO TO 420
01326   330 KK = 0
01327       DO 340 I=1,KLM
01328          Z = Q(I,IN)
01329          IF (Z.LE.TOLER) GO TO 340
01330          KK = KK + 1
01331          RES(KK) = Q(I,N1)/Z
01332          S(KK) = I
01333   340 CONTINUE
01334   350 IF (KK.GT.0) GO TO 360
01335       KODE = 2
01336       GO TO 590
01337   360 XMIN = RES(1)
01338       IOUT = S(1)
01339       J = 1
01340       IF (KK.EQ.1) GO TO 380
01341       DO 370 I=2,KK
01342          IF (RES(I).GE.XMIN) GO TO 370
01343          J = I
01344          XMIN = RES(I)
01345          IOUT = S(I)
01346   370 CONTINUE
01347       RES(J) = RES(KK)
01348       S(J) = S(KK)
01349   380 KK = KK - 1
01350       PIVOT = Q(IOUT,IN)
01351       II = Q(IOUT,N2)
01352       IF (IPHASE.EQ.1) GO TO 400
01353       IF (II.LT.0) GO TO 390
01354       IF (IU(2,II).EQ.1) GO TO 420
01355       GO TO 400
01356   390 IINEG = -II
01357       IF (IU(1,IINEG).EQ.1) GO TO 420
01358   400 II = IABS(II)
01359       CUV = CU(1,II) + CU(2,II)
01360       IF (Q(KLM1,IN)-PIVOT*CUV.LE.TOLER) GO TO 420
01361 C BYPASS INTERMEDIATE VERTICES.
01362       DO 410 J=JS,N1
01363          Z = Q(IOUT,J)
01364          Q(KLM1,J) = Q(KLM1,J) - Z*CUV
01365          Q(IOUT,J) = -Z
01366   410 CONTINUE
01367       Q(IOUT,N2) = -Q(IOUT,N2)
01368       GO TO 350
01369 C GAUSS-JORDAN ELIMINATION.
01370   420 IF (ITER.LT.MAXIT) GO TO 430
01371       KODE = 3
01372       GO TO 590
01373   430 ITER = ITER + 1
01374       DO 440 J=JS,N1
01375          IF (J.NE.IN) Q(IOUT,J) = Q(IOUT,J)/PIVOT
01376   440 CONTINUE
01377 C IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION
01378 C SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN
01379 C TO AND INCLUDING STATEMENT NUMBER 460 BY..
01380 C     DO 460 J=JS,N1
01381 C        IF(J .EQ. IN) GO TO 460
01382 C        Z = -Q(IOUT,J)
01383 C        CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1)
01384 C 460 CONTINUE
01385       DO 460 J=JS,N1
01386          IF (J.EQ.IN) GO TO 460
01387          Z = -Q(IOUT,J)
01388          DO 450 I=1,KLM1
01389             IF (I.NE.IOUT) Q(I,J) = Q(I,J) + Z*Q(I,IN)
01390   450    CONTINUE
01391   460 CONTINUE
01392       TPIVOT = -PIVOT
01393       DO 470 I=1,KLM1
01394          IF (I.NE.IOUT) Q(I,IN) = Q(I,IN)/TPIVOT
01395   470 CONTINUE
01396       Q(IOUT,IN) = 1./PIVOT
01397       Z = Q(IOUT,N2)
01398       Q(IOUT,N2) = Q(KLM2,IN)
01399       Q(KLM2,IN) = Z
01400       II = ABS(Z)
01401       IF (IU(1,II).EQ.0 .OR. IU(2,II).EQ.0) GO TO 240
01402       DO 480 I=1,KLM2
01403          Z = Q(I,IN)
01404          Q(I,IN) = Q(I,JS)
01405          Q(I,JS) = Z
01406   480 CONTINUE
01407       JS = JS + 1
01408       GO TO 240
01409 C TEST FOR OPTIMALITY.
01410   490 IF (KFORCE.EQ.0) GO TO 580
01411       IF (IPHASE.EQ.1 .AND. Q(KLM1,N1).LE.TOLER) GO TO 500
01412       KFORCE = 0
01413       GO TO 240
01414 C SET UP PHASE 2 COSTS.
01415   500 IPHASE = 2
01416       DO 510 J=1,NKLM
01417          CU(1,J) = 0.
01418          CU(2,J) = 0.
01419   510 CONTINUE
01420       DO 520 J=N1,NK
01421          CU(1,J) = 1.
01422          CU(2,J) = 1.
01423   520 CONTINUE
01424       DO 560 I=1,KLM
01425          II = Q(I,N2)
01426          IF (II.GT.0) GO TO 530
01427          II = -II
01428          IF (IU(2,II).EQ.0) GO TO 560
01429          CU(2,II) = 0.
01430          GO TO 540
01431   530    IF (IU(1,II).EQ.0) GO TO 560
01432          CU(1,II) = 0.
01433   540    IA = IA + 1
01434          DO 550 J=1,N2
01435             Z = Q(IA,J)
01436             Q(IA,J) = Q(I,J)
01437             Q(I,J) = Z
01438   550    CONTINUE
01439   560 CONTINUE
01440       GO TO 160
01441   570 IF (Q(KLM1,N1).LE.TOLER) GO TO 500
01442       KODE = 1
01443       GO TO 590
01444   580 IF (IPHASE.EQ.1) GO TO 570
01445 C PREPARE OUTPUT.
01446       KODE = 0
01447   590 SUM = 0.D0
01448       DO 600 J=1,N
01449          X(J) = 0.
01450   600 CONTINUE
01451       DO 610 I=1,KLM
01452          RES(I) = 0.
01453   610 CONTINUE
01454       DO 640 I=1,KLM
01455          II = Q(I,N2)
01456          SN = 1.
01457          IF (II.GT.0) GO TO 620
01458          II = -II
01459          SN = -1.
01460   620    IF (II.GT.N) GO TO 630
01461          X(II) = SN*Q(I,N1)
01462          GO TO 640
01463   630    IIMN = II - N
01464          RES(IIMN) = SN*Q(I,N1)
01465          IF (II.GE.N1 .AND. II.LE.NK) SUM = SUM +
01466      *    DBLE(Q(I,N1))
01467   640 CONTINUE
01468       ERROR = SUM
01469       RETURN
01470       END
01471 ================================================================*/