PLearn 0.1
qld.cc
Go to the documentation of this file.
00001 
00002 /* ************************************************************/
00003 
00004 /*  -- translated by f2c (version of 22 July 1992  22:54:52).
00005  */
00006 
00007 /* umd
00008    Must include math.h before f2c.h - f2c does a #define abs.
00009    (Thanks go to Martin Wauchope for providing this correction)
00010    We manually included AT&T's f2c.h in this source file, i.e.
00011    it does not have to be present separately in order to compile.
00012 */
00013 #include <math.h>
00014 
00015 /* CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
00016    !!!! NOTICE !!!! 
00017 
00018    1. The routines contained in this file are due to Prof. K.Schittkowski 
00019    of the University of Bayreuth, Germany (modification of routines
00020    due to Prof. MJD Powell at the University of Cambridge).  They can 
00021    be freely distributed.
00022 
00023    2. A few minor modifications were performed at the University of
00024    Maryland. They are marked in the code by "umd".
00025 
00026    A.L. Tits, J.L. Zhou, and
00027    Craig Lawrence
00028    University of Maryland
00029 
00030    ***********************************************************************
00031 
00032 
00033 
00034    SOLUTION OF QUADRATIC PROGRAMMING PROBLEMS
00035 
00036 
00037 
00038    QL0001 SOLVES THE QUADRATIC PROGRAMMING PROBLEM
00039 
00040    MINIMIZE        .5*X'*C*X + D'*X 
00041    SUBJECT TO      A(J)*X  +  B(J)   =  0  ,  J=1,...,ME
00042    A(J)*X  +  B(J)  >=  0  ,  J=ME+1,...,M
00043    XL  <=  X  <=  XU
00044 
00045    HERE C MUST BE AN N BY N SYMMETRIC AND POSITIVE MATRIX, D AN N-DIMENSIONAL
00046    VECTOR, A AN M BY N MATRIX AND B AN M-DIMENSIONAL VECTOR. THE ABOVE 
00047    SITUATION IS INDICATED BY IWAR(1)=1. ALTERNATIVELY, I.E. IF IWAR(1)=0,
00048    THE OBJECTIVE FUNCTION MATRIX CAN ALSO BE PROVIDED IN FACTORIZED FORM.
00049    IN THIS CASE, C IS AN UPPER TRIANGULAR MATRIX.
00050 
00051    THE SUBROUTINE REORGANIZES SOME DATA SO THAT THE PROBLEM CAN BE SOLVED
00052    BY A MODIFICATION OF AN ALGORITHM PROPOSED BY POWELL (1983).
00053 
00054 
00055    USAGE:
00056 
00057    QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,X,U,IOUT,IFAIL,IPRINT, 
00058    WAR,LWAR,IWAR,LIWAR) 
00059 
00060 
00061    DEFINITION OF THE PARAMETERS: 
00062 
00063    M :        TOTAL NUMBER OF CONSTRAINTS. 
00064    ME :       NUMBER OF EQUALITY CONSTRAINTS.
00065    MMAX :     ROW DIMENSION OF A. MMAX MUST BE AT LEAST ONE AND GREATER 
00066    THAN M.
00067    N :        NUMBER OF VARIABLES.
00068    NMAX :     ROW DIMENSION OF C. NMAX MUST BE GREATER OR EQUAL TO N.
00069    MNN :      MUST BE EQUAL TO M + N + N. 
00070    C(NMAX,NMAX): OBJECTIVE FUNCTION MATRIX WHICH SHOULD BE SYMMETRIC AND
00071    POSITIVE DEFINITE. IF IWAR(1) = 0, C IS SUPPOSED TO BE THE
00072    CHOLESKEY-FACTOR OF ANOTHER MATRIX, I.E. C IS UPPER
00073    TRIANGULAR.
00074    D(NMAX) :  CONTAINS THE CONSTANT VECTOR OF THE OBJECTIVE FUNCTION.
00075    A(MMAX,NMAX): CONTAINS THE DATA MATRIX OF THE LINEAR CONSTRAINTS. 
00076    B(MMAX) :  CONTAINS THE CONSTANT DATA OF THE LINEAR CONSTRAINTS. 
00077    XL(N),XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR THE VARIABLES.
00078    X(N) :     ON RETURN, X CONTAINS THE OPTIMAL SOLUTION VECTOR.
00079    U(MNN) :   ON RETURN, U CONTAINS THE LAGRANGE MULTIPLIERS. THE FIRST 
00080    M POSITIONS ARE RESERVED FOR THE MULTIPLIERS OF THE M
00081    LINEAR CONSTRAINTS AND THE SUBSEQUENT ONES FOR THE
00082    MULTIPLIERS OF THE LOWER AND UPPER BOUNDS. ON SUCCESSFUL 
00083    TERMINATION, ALL VALUES OF U WITH RESPECT TO INEQUALITIES 
00084    AND BOUNDS SHOULD BE GREATER OR EQUAL TO ZERO.
00085    IOUT :     INTEGER INDICATING THE DESIRED OUTPUT UNIT NUMBER, I.E.
00086    ALL WRITE-STATEMENTS START WITH 'WRITE(IOUT,... '.
00087    IFAIL :    SHOWS THE TERMINATION REASON.
00088    IFAIL = 0 :   SUCCESSFUL RETURN.
00089    IFAIL = 1 :   TOO MANY ITERATIONS (MORE THAN 40*(N+M)).
00090    IFAIL = 2 :   ACCURACY INSUFFICIENT TO SATISFY CONVERGENCE
00091    CRITERION.
00092    IFAIL = 5 :   LENGTH OF A WORKING ARRAY IS TOO SHORT.
00093    IFAIL > 10 :  THE CONSTRAINTS ARE INCONSISTENT.
00094    IPRINT :   OUTPUT CONTROL.
00095    IPRINT = 0 :  NO OUTPUT OF QL0001.
00096    IPRINT > 0 :  BRIEF OUTPUT IN ERROR CASES.
00097    WAR(LWAR) : REAL WORKING ARRAY. THE LENGTH LWAR SHOULD BE GRATER THAN
00098    3*NMAX*NMAX/2 + 10*NMAX + 2*MMAX.
00099    IWAR(LIWAR): INTEGER WORKING ARRAY. THE LENGTH LIWAR SHOULD BE AT
00100    LEAST N.
00101    IF IWAR(1)=0 INITIALLY, THEN THE CHOLESKY DECOMPOSITION
00102    WHICH IS REQUIRED BY THE DUAL ALGORITHM TO GET THE FIRST 
00103    UNCONSTRAINED MINIMUM OF THE OBJECTIVE FUNCTION, IS
00104    PERFORMED INTERNALLY. OTHERWISE, I.E. IF IWAR(1)=1, THEN 
00105    IT IS ASSUMED THAT THE USER PROVIDES THE INITIAL FAC-
00106    TORIZATION BY HIMSELF AND STORES IT IN THE UPPER TRIAN-
00107    GULAR PART OF THE ARRAY C.
00108 
00109    A NAMED COMMON-BLOCK  /CMACHE/EPS   MUST BE PROVIDED BY THE USER,
00110    WHERE EPS DEFINES A GUESS FOR THE UNDERLYING MACHINE PRECISION.
00111 
00112 
00113    AUTHOR:    K. SCHITTKOWSKI,
00114    MATHEMATISCHES INSTITUT,
00115    UNIVERSITAET BAYREUTH,
00116    8580 BAYREUTH,
00117    GERMANY, F.R.
00118 
00119 
00120    VERSION:   1.4  (MARCH, 1987)
00121 */
00122 /* f2c.h  --  Standard Fortran to C header file */
00123 
00128 #ifndef F2C_INCLUDE
00129 #define F2C_INCLUDE
00130 
00131 typedef int integer;
00132 typedef char *address;
00133 typedef short int shortint;
00134 typedef float real;
00135 typedef double doublereal;
00136 typedef struct { real r, i; } complex;
00137 typedef struct { doublereal r, i; } doublecomplex;
00138 typedef long int logical;
00139 typedef short int shortlogical;
00140 
00141 #define TRUE_ (1)
00142 #define FALSE_ (0)
00143 
00144 /* Extern is for use with -E */
00145 #ifndef Extern
00146 #define Extern extern
00147 #endif
00148 
00149 /* I/O stuff */
00150 
00151 #ifdef f2c_i2
00152 /* for -i2 */
00153 typedef short flag;
00154 typedef short ftnlen;
00155 typedef short ftnint;
00156 #else
00157 typedef long flag;
00158 typedef long ftnlen;
00159 typedef long ftnint;
00160 #endif
00161 
00162 /*external read, write*/
00163 typedef struct
00164 {       flag cierr;
00165     ftnint ciunit;
00166     flag ciend;
00167     char *cifmt;
00168     ftnint cirec;
00169 } cilist;
00170 
00171 /*internal read, write*/
00172 typedef struct
00173 {       flag icierr;
00174     char *iciunit;
00175     flag iciend;
00176     char *icifmt;
00177     ftnint icirlen;
00178     ftnint icirnum;
00179 } icilist;
00180 
00181 /*open*/
00182 typedef struct
00183 {       flag oerr;
00184     ftnint ounit;
00185     char *ofnm;
00186     ftnlen ofnmlen;
00187     char *osta;
00188     char *oacc;
00189     char *ofm;
00190     ftnint orl;
00191     char *oblnk;
00192 } olist;
00193 
00194 /*close*/
00195 typedef struct
00196 {       flag cerr;
00197     ftnint cunit;
00198     char *csta;
00199 } cllist;
00200 
00201 /*rewind, backspace, endfile*/
00202 typedef struct
00203 {       flag aerr;
00204     ftnint aunit;
00205 } alist;
00206 
00207 /* inquire */
00208 typedef struct
00209 {       flag inerr;
00210     ftnint inunit;
00211     char *infile;
00212     ftnlen infilen;
00213     ftnint      *inex;  /*parameters in standard's order*/
00214     ftnint      *inopen;
00215     ftnint      *innum;
00216     ftnint      *innamed;
00217     char        *inname;
00218     ftnlen      innamlen;
00219     char        *inacc;
00220     ftnlen      inacclen;
00221     char        *inseq;
00222     ftnlen      inseqlen;
00223     char        *indir;
00224     ftnlen      indirlen;
00225     char        *infmt;
00226     ftnlen      infmtlen;
00227     char        *inform;
00228     ftnint      informlen;
00229     char        *inunf;
00230     ftnlen      inunflen;
00231     ftnint      *inrecl;
00232     ftnint      *innrec;
00233     char        *inblank;
00234     ftnlen      inblanklen;
00235 } inlist;
00236 
00237 #define VOID void
00238 
00239 union Multitype {       /* for multiple entry points */
00240     shortint h;
00241     integer i;
00242     real r;
00243     doublereal d;
00244     complex c;
00245     doublecomplex z;
00246 };
00247 
00248 typedef union Multitype Multitype;
00249 
00250 typedef long Long;
00251 
00252 struct Vardesc {        /* for Namelist */
00253     char *name;
00254     char *addr;
00255     Long *dims;
00256     int  type;
00257 };
00258 typedef struct Vardesc Vardesc;
00259 
00260 struct Namelist {
00261     char *name;
00262     Vardesc **vars;
00263     int nvars;
00264 };
00265 typedef struct Namelist Namelist;
00266 
00267 #define abs(x) ((x) >= 0 ? (x) : -(x))
00268 #define dabs(x) (doublereal)abs(x)
00269 #define min(a,b) ((a) <= (b) ? (a) : (b))
00270 #define max(a,b) ((a) >= (b) ? (a) : (b))
00271 #define dmin(a,b) (doublereal)min(a,b)
00272 #define dmax(a,b) (doublereal)max(a,b)
00273 
00274 /* procedure parameter types for -A and -C++ */
00275 
00276 #define F2C_proc_par_types 1
00277 #ifdef __cplusplus
00278 typedef int /* Unknown procedure type */ (*U_fp)(...);
00279 typedef shortint (*J_fp)(...);
00280 typedef integer (*I_fp)(...);
00281 typedef real (*R_fp)(...);
00282 typedef doublereal (*D_fp)(...), (*E_fp)(...);
00283 typedef /* Complex */ VOID (*C_fp)(...);
00284 typedef /* Double Complex */ VOID (*Z_fp)(...);
00285 typedef logical (*L_fp)(...);
00286 typedef shortlogical (*K_fp)(...);
00287 typedef /* Character */ VOID (*H_fp)(...);
00288 typedef /* Subroutine */ int (*S_fp)(...);
00289 #else
00290 typedef int /* Unknown procedure type */ (*U_fp)();
00291 typedef shortint (*J_fp)();
00292 typedef integer (*I_fp)();
00293 typedef real (*R_fp)();
00294 typedef doublereal (*D_fp)(), (*E_fp)();
00295 typedef /* Complex */ VOID (*C_fp)();
00296 typedef /* Double Complex */ VOID (*Z_fp)();
00297 typedef logical (*L_fp)();
00298 typedef shortlogical (*K_fp)();
00299 typedef /* Character */ VOID (*H_fp)();
00300 typedef /* Subroutine */ int (*S_fp)();
00301 #endif
00302 /* E_fp is for real functions when -R is not specified */
00303 typedef VOID C_f;       /* complex function */
00304 typedef VOID H_f;       /* character function */
00305 typedef VOID Z_f;       /* double complex function */
00306 typedef doublereal E_f; /* real function with -R not specified */
00307 
00308 /* undef any lower-case symbols that your C compiler predefines, e.g.: */
00309 
00310 #ifndef Skip_f2c_Undefs
00311 #undef cray
00312 #undef gcos
00313 #undef mc68010
00314 #undef mc68020
00315 #undef mips
00316 #undef pdp11
00317 #undef sgi
00318 #undef sparc
00319 #undef sun
00320 #undef sun2
00321 #undef sun3
00322 #undef sun4
00323 #undef u370
00324 #undef u3b
00325 #undef u3b2
00326 #undef u3b5
00327 #undef unix
00328 #undef vax
00329 #endif
00330 #endif
00331 
00332 
00333 
00334 /* Common Block Declarations */
00335 
00336 struct {
00337     doublereal eps;
00338 } cmache_;
00339 
00340 #define cmache_1 cmache_
00341 
00342 /* Table of constant values */
00343 
00344 static integer c__1 = 1;
00345 
00346 /* umd */
00347 /*
00348   ql0002_ is declared here to provide ANSI C compliance.
00349   (Thanks got to Martin Wauchope for providing this correction)
00350 */
00351 #ifdef __STDC__
00352 
00353 int ql0002_(integer *n,integer *m,integer *meq,integer *mmax,
00354             integer *mn,integer *mnn,integer *nmax,
00355             logical *lql,
00356             doublereal *a,doublereal *b,doublereal *grad,
00357             doublereal *g,doublereal *xl,doublereal *xu,doublereal *x,
00358             integer *nact,integer *iact,integer *maxit,
00359             doublereal *vsmall,
00360             integer *info,
00361             doublereal *diag, doublereal *w,
00362             integer *lw);
00363 #else
00364 int ql0002_();
00365 #endif
00366 /* umd */
00367 /*
00368   When the fortran code was f2c converted, the use of fortran COMMON
00369   blocks was no longer available. Thus an additional variable, eps1,
00370   was added to the parameter list to account for this.
00371 */
00372 /* umd */
00373 /* 
00374    Two alternative definitions are provided in order to give ANSI
00375    compliance.
00376 */
00377 #ifdef __STDC__
00378 int ql0001_(int *m,int *me,int *mmax,int *n,int *nmax,int *mnn,
00379             double *c,double *d,double *a,double *b,double *xl,
00380             double *xu,double *x,double *u,int *iout,int *ifail,
00381             int *iprint,double *war,int *lwar,int *iwar,int *liwar,
00382             double *eps1)
00383 #else
00384 /* Subroutine */ 
00385     int ql0001_(m, me, mmax, n, nmax, mnn, c, d, a, b, xl, xu, x,
00386                 u, iout, ifail, iprint, war, lwar, iwar, liwar, eps1)
00387     integer *m, *me, *mmax, *n, *nmax, *mnn;
00388 doublereal *c, *d, *a, *b, *xl, *xu, *x, *u;
00389 integer *iout, *ifail, *iprint;
00390 doublereal *war;
00391 integer *lwar, *iwar, *liwar;
00392 doublereal *eps1;
00393 #endif
00394 {
00395     /* Format strings */
00396     static char fmt_1000[] = "(/8x,\002***QL: MATRIX G WAS ENLARGED\002,i3\
00397 ,\002-TIMES BY UNITMATRIX\002)";
00398     static char fmt_1100[] = "(/8x,\002***QL: CONSTRAINT \002,i5,\002 NOT CO\
00399 NSISTENT TO \002,/,(10x,10i5))";
00400     static char fmt_1200[] = "(/8x,\002***QL: LWAR TOO SMALL\002)";
00401     static char fmt_1210[] = "(/8x,\002***QL: LIWAR TOO SMALL\002)";
00402     static char fmt_1220[] = "(/8x,\002***QL: MNN TOO SMALL\002)";
00403     static char fmt_1300[] = "(/8x,\002***QL: TOO MANY ITERATIONS (MORE THA\
00404 N\002,i6,\002)\002)";
00405     static char fmt_1400[] = "(/8x,\002***QL: ACCURACY INSUFFICIENT TO ATTAI\
00406 N CONVERGENCE\002)";
00407 
00408     /* System generated locals */
00409     integer c_dim1, c_offset, a_dim1, a_offset, i__1, i__2;
00410 
00411     /* Builtin functions */
00412 /*    integer s_wsfe(), do_fio(), e_wsfe(); */
00413 
00414     /* Local variables */
00415     static doublereal diag;
00416     /* extern int ql0002_(); */
00417     static integer nact, info;
00418     static doublereal zero;
00419     static integer i, j, idiag, maxit;
00420     static doublereal qpeps;
00421     static integer in, mn, lw;
00422     static doublereal ten;
00423     static logical lql;
00424     static integer inw1, inw2;
00425 
00426     /* Fortran I/O blocks */
00427     static cilist io___16 = { 0, 0, 0, fmt_1000, 0 };
00428     static cilist io___18 = { 0, 0, 0, fmt_1100, 0 };
00429     static cilist io___19 = { 0, 0, 0, fmt_1200, 0 };
00430     static cilist io___20 = { 0, 0, 0, fmt_1210, 0 };
00431     static cilist io___21 = { 0, 0, 0, fmt_1220, 0 };
00432     static cilist io___22 = { 0, 0, 0, fmt_1300, 0 };
00433     static cilist io___23 = { 0, 0, 0, fmt_1400, 0 };
00434 
00435 
00436 
00437 
00438 
00439 /*     INTRINSIC FUNCTIONS:  DSQRT */
00440 
00441     /* Parameter adjustments */
00442     --iwar;
00443     --war;
00444     --u;
00445     --x;
00446     --xu;
00447     --xl;
00448     --b;
00449     a_dim1 = *mmax;
00450     a_offset = a_dim1 + 1;
00451     a -= a_offset;
00452     --d;
00453     c_dim1 = *nmax;
00454     c_offset = c_dim1 + 1;
00455     c -= c_offset;
00456 
00457     /* Function Body */
00458     cmache_1.eps = *eps1;
00459 
00460 /*     CONSTANT DATA */
00461 
00462 /* ################################################################# */
00463 
00464     if (fabs(c[*nmax + *nmax * c_dim1]) == 0.e0) {
00465         c[*nmax + *nmax * c_dim1] = cmache_1.eps;
00466     }
00467 
00468 /* umd */
00469 /*  This prevents a subsequent more major modification of the Hessian */
00470 /*  matrix in the important case when a minmax problem (yielding a */
00471 /*  singular Hessian matrix) is being solved. */
00472 /*                                 ----UMCP, April 1991, Jian L. Zhou */
00473 /* ################################################################# */
00474 
00475     lql = FALSE_;
00476     if (iwar[1] == 1) {
00477         lql = TRUE_;
00478     }
00479     zero = 0.;
00480     ten = 10.;
00481     maxit = (*m + *n) * 40;
00482     qpeps = cmache_1.eps;
00483     inw1 = 1;
00484     inw2 = inw1 + *mmax;
00485 
00486 /*     PREPARE PROBLEM DATA FOR EXECUTION */
00487 
00488     if (*m <= 0) {
00489         goto L20;
00490     }
00491     in = inw1;
00492     i__1 = *m;
00493     for (j = 1; j <= i__1; ++j) {
00494         war[in] = -b[j];
00495 /* L10: */
00496         ++in;
00497     }
00498  L20:
00499     lw = *nmax * 3 * *nmax / 2 + *nmax * 10 + *m;
00500     if (inw2 + lw > *lwar) {
00501         goto L80;
00502     }
00503     if (*liwar < *n) {
00504         goto L81;
00505     }
00506     if (*mnn < *m + *n + *n) {
00507         goto L82;
00508     }
00509     mn = *m + *n;
00510 
00511 /*     CALL OF QL0002 */
00512 
00513     ql0002_(n, m, me, mmax, &mn, mnn, nmax, &lql, &a[a_offset], &war[inw1], &
00514             d[1], &c[c_offset], &xl[1], &xu[1], &x[1], &nact, &iwar[1], &
00515             maxit, &qpeps, &info, &diag, &war[inw2], &lw);
00516 
00517 /*     TEST OF MATRIX CORRECTIONS */
00518 
00519     *ifail = 0;
00520     if (info == 1) {
00521         goto L40;
00522     }
00523     if (info == 2) {
00524         goto L90;
00525     }
00526     idiag = 0;
00527     if (diag > zero && diag < 1e3) {
00528         idiag = (integer) diag;
00529     }
00530 /*
00531   if (*iprint > 0 && idiag > 0) {
00532   io___16.ciunit = *iout;
00533   s_wsfe(&io___16);
00534   do_fio(&c__1, (char *)&idiag, (ftnlen)sizeof(integer));
00535   e_wsfe();
00536   }
00537 */
00538     if (info < 0) {
00539         goto L70;
00540     }
00541 
00542 /*     REORDER MULTIPLIER */
00543 
00544     i__1 = *mnn;
00545     for (j = 1; j <= i__1; ++j) {
00546 /* L50: */
00547         u[j] = zero;
00548     }
00549     in = inw2 - 1;
00550     if (nact == 0) {
00551         goto L30;
00552     }
00553     i__1 = nact;
00554     for (i = 1; i <= i__1; ++i) {
00555         j = iwar[i];
00556         u[j] = war[in + i];
00557 /* L60: */
00558     }
00559  L30:
00560     return 0;
00561 
00562 /*     ERROR MESSAGES */
00563 
00564  L70:
00565     *ifail = -info + 10;
00566 /*
00567   if (*iprint > 0 && nact > 0) {
00568   io___18.ciunit = *iout;
00569   s_wsfe(&io___18);
00570   i__1 = -info;
00571   do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer));
00572   i__2 = nact;
00573   for (i = 1; i <= i__2; ++i) {
00574   do_fio(&c__1, (char *)&iwar[i], (ftnlen)sizeof(integer));
00575   }
00576   e_wsfe();
00577   }
00578 */
00579     return 0;
00580  L80:
00581     *ifail = 5;
00582 /*
00583   if (*iprint > 0) {
00584   io___19.ciunit = *iout;
00585   s_wsfe(&io___19);
00586   e_wsfe();
00587   }
00588 */
00589     return 0;
00590  L81:
00591     *ifail = 5;
00592 /*
00593   if (*iprint > 0) {
00594   io___20.ciunit = *iout;
00595   s_wsfe(&io___20);
00596   e_wsfe();
00597   }
00598 */
00599     return 0;
00600  L82:
00601     *ifail = 5;
00602 /*
00603   if (*iprint > 0) {
00604   io___21.ciunit = *iout;
00605   s_wsfe(&io___21);
00606   e_wsfe();
00607   }
00608 */
00609     return 0;
00610  L40:
00611     *ifail = 1;
00612 /*
00613   if (*iprint > 0) {
00614   io___22.ciunit = *iout;
00615   s_wsfe(&io___22);
00616   do_fio(&c__1, (char *)&maxit, (ftnlen)sizeof(integer));
00617   e_wsfe();
00618   }
00619 */
00620     return 0;
00621  L90:
00622     *ifail = 2;
00623 /*
00624   if (*iprint > 0) {
00625   io___23.ciunit = *iout;
00626   s_wsfe(&io___23);
00627   e_wsfe();
00628   }
00629 */
00630     return 0;
00631 
00632 /*     FORMAT-INSTRUCTIONS */
00633 
00634 } /* ql0001_ */
00635 
00636 
00637 /* umd
00638    Two alternative definitions are provided in order to give ANSI
00639    compliance.
00640    (Thanks got to Martin Wauchope for providing this correction)
00641 */
00642 #ifdef __STDC__
00643 int ql0002_(integer *n,integer *m,integer *meq,integer *mmax,
00644             integer *mn,integer *mnn,integer *nmax,
00645             logical *lql,
00646             doublereal *a,doublereal *b,doublereal *grad,
00647             doublereal *g,doublereal *xl,doublereal *xu,doublereal *x,
00648             integer *nact,integer *iact,integer *maxit,
00649             doublereal *vsmall,
00650             integer *info,
00651             doublereal *diag, doublereal *w,
00652             integer *lw)
00653 #else
00654     /* Subroutine */ int ql0002_(n, m, meq, mmax, mn, mnn, nmax, lql, a, b, grad, 
00655                                  g, xl, xu, x, nact, iact, maxit, vsmall, info, diag, w, lw)
00656     integer *n, *m, *meq, *mmax, *mn, *mnn, *nmax;
00657 logical *lql;
00658 doublereal *a, *b, *grad, *g, *xl, *xu, *x;
00659 integer *nact, *iact, *maxit;
00660 doublereal *vsmall;
00661 integer *info;
00662 doublereal *diag, *w;
00663 integer *lw;
00664 #endif
00665 {
00666     /* System generated locals */
00667     integer a_dim1, a_offset, g_dim1, g_offset, i__1, i__2, i__3, i__4;
00668     doublereal d__1, d__2, d__3, d__4;
00669 
00670     /* Builtin functions */
00671     /* umd */
00672     /* double sqrt();    */
00673 
00674     /* Local variables */
00675     static doublereal onha, xmag, suma, sumb, sumc, temp, step, zero;
00676     static integer iwwn;
00677     static doublereal sumx, sumy;
00678     static integer i, j, k;
00679     static doublereal fdiff;
00680     static integer iflag, jflag, kflag, lflag;
00681     static doublereal diagr;
00682     static integer ifinc, kfinc, jfinc, mflag, nflag;
00683     static doublereal vfact, tempa;
00684     static integer iterc, itref;
00685     static doublereal cvmax, ratio, xmagr;
00686     static integer kdrop;
00687     static logical lower;
00688     static integer knext, k1;
00689     static doublereal ga, gb;
00690     static integer ia, id;
00691     static doublereal fdiffa;
00692     static integer ii, il, kk, jl, ip, ir, nm, is, iu, iw, ju, ix, iz, nu, iy;
00693 
00694     static doublereal parinc, parnew;
00695     static integer ira, irb, iwa;
00696     static doublereal one;
00697     static integer iwd, iza;
00698     static doublereal res;
00699     static integer ipp, iwr, iws;
00700     static doublereal sum;
00701     static integer iww, iwx, iwy;
00702     static doublereal two;
00703     static integer iwz;
00704 
00705 
00706 /*       WHETHER THE CONSTRAINT IS ACTIVE. */
00707 
00708 
00709 /*   AUTHOR:    K. SCHITTKOWSKI, */
00710 /*              MATHEMATISCHES INSTITUT, */
00711 /*              UNIVERSITAET BAYREUTH, */
00712 /*              8580 BAYREUTH, */
00713 /*              GERMANY, F.R. */
00714 
00715 /*   AUTHOR OF ORIGINAL VERSION: */
00716 /*              M.J.D. POWELL, DAMTP, */
00717 /*              UNIVERSITY OF CAMBRIDGE, SILVER STREET */
00718 /*              CAMBRIDGE, */
00719 /*              ENGLAND */
00720 
00721 
00722 /*   REFERENCE: M.J.D. POWELL: ZQPCVX, A FORTRAN SUBROUTINE FOR CONVEX */
00723 /*              PROGRAMMING, REPORT DAMTP/1983/NA17, UNIVERSITY OF */
00724 /*              CAMBRIDGE, ENGLAND, 1983. */
00725 
00726 
00727 /*   VERSION :  2.0 (MARCH, 1987) */
00728 
00729 
00730 /************************************************************************
00731  ***/
00732 
00733 
00734 /*   INTRINSIC FUNCTIONS:   DMAX1,DSQRT,DABS,DMIN1 */
00735 
00736 
00737 /*   INITIAL ADDRESSES */
00738 
00739     /* Parameter adjustments */
00740     --w;
00741     --iact;
00742     --x;
00743     --xu;
00744     --xl;
00745     g_dim1 = *nmax;
00746     g_offset = g_dim1 + 1;
00747     g -= g_offset;
00748     --grad;
00749     --b;
00750     a_dim1 = *mmax;
00751     a_offset = a_dim1 + 1;
00752     a -= a_offset;
00753 
00754     /* Function Body */
00755     iwz = *nmax;
00756     iwr = iwz + *nmax * *nmax;
00757     iww = iwr + *nmax * (*nmax + 3) / 2;
00758     iwd = iww + *nmax;
00759     iwx = iwd + *nmax;
00760     iwa = iwx + *nmax;
00761 
00762 /*     SET SOME CONSTANTS. */
00763 
00764     zero = 0.;
00765     one = 1.;
00766     two = 2.;
00767     onha = 1.5;
00768     vfact = 1.;
00769 
00770 /*     SET SOME PARAMETERS. */
00771 /*     NUMBER LESS THAN VSMALL ARE ASSUMED TO BE NEGLIGIBLE. */
00772 /*     THE MULTIPLE OF I THAT IS ADDED TO G IS AT MOST DIAGR TIMES */
00773 /*       THE LEAST MULTIPLE OF I THAT GIVES POSITIVE DEFINITENESS. */
00774 /*     X IS RE-INITIALISED IF ITS MAGNITUDE IS REDUCED BY THE */
00775 /*       FACTOR XMAGR. */
00776 /*     A CHECK IS MADE FOR AN INCREASE IN F EVERY IFINC ITERATIONS, */
00777 /*       AFTER KFINC ITERATIONS ARE COMPLETED. */
00778 
00779     diagr = two;
00780     xmagr = .01;
00781     ifinc = 3;
00782     kfinc = max(10,*n);
00783 
00784 /*     FIND THE RECIPROCALS OF THE LENGTHS OF THE CONSTRAINT NORMALS. */
00785 /*     RETURN IF A CONSTRAINT IS INFEASIBLE DUE TO A ZERO NORMAL. */
00786 
00787     *nact = 0;
00788     if (*m <= 0) {
00789         goto L45;
00790     }
00791     i__1 = *m;
00792     for (k = 1; k <= i__1; ++k) {
00793         sum = zero;
00794         i__2 = *n;
00795         for (i = 1; i <= i__2; ++i) {
00796 /* L10: */
00797 /* Computing 2nd power */
00798             d__1 = a[k + i * a_dim1];
00799             sum += d__1 * d__1;
00800         }
00801         if (sum > zero) {
00802             goto L20;
00803         }
00804         if (b[k] == zero) {
00805             goto L30;
00806         }
00807         *info = -k;
00808         if (k <= *meq) {
00809             goto L730;
00810         }
00811         if (b[k] <= 0.) {
00812             goto L30;
00813         } else {
00814             goto L730;
00815         }
00816     L20:
00817         sum = one / sqrt(sum);
00818     L30:
00819         ia = iwa + k;
00820 /* L40: */
00821         w[ia] = sum;
00822     }
00823  L45:
00824     i__1 = *n;
00825     for (k = 1; k <= i__1; ++k) {
00826         ia = iwa + *m + k;
00827 /* L50: */
00828         w[ia] = one;
00829     }
00830 
00831 /*     IF NECESSARY INCREASE THE DIAGONAL ELEMENTS OF G. */
00832 
00833     if (! (*lql)) {
00834         goto L165;
00835     }
00836     *diag = zero;
00837     i__1 = *n;
00838     for (i = 1; i <= i__1; ++i) {
00839         id = iwd + i;
00840         w[id] = g[i + i * g_dim1];
00841 /* Computing MAX */
00842         d__1 = *diag, d__2 = *vsmall - w[id];
00843         *diag = max(d__1,d__2);
00844         if (i == *n) {
00845             goto L60;
00846         }
00847         ii = i + 1;
00848         i__2 = *n;
00849         for (j = ii; j <= i__2; ++j) {
00850 /* Computing MIN */
00851             d__1 = w[id], d__2 = g[j + j * g_dim1];
00852             ga = -min(d__1,d__2);
00853             gb = (d__1 = w[id] - g[j + j * g_dim1], abs(d__1)) + (d__2 = g[i 
00854                                                                            + j * g_dim1], abs(d__2));
00855             if (gb > zero) {
00856 /* Computing 2nd power */
00857                 d__1 = g[i + j * g_dim1];
00858                 ga += d__1 * d__1 / gb;
00859             }
00860 /* L55: */
00861             *diag = max(*diag,ga);
00862         }
00863     L60:
00864         ;
00865     }
00866     if (*diag <= zero) {
00867         goto L90;
00868     }
00869  L70:
00870     *diag = diagr * *diag;
00871     i__1 = *n;
00872     for (i = 1; i <= i__1; ++i) {
00873         id = iwd + i;
00874 /* L80: */
00875         g[i + i * g_dim1] = *diag + w[id];
00876     }
00877 
00878 /*     FORM THE CHOLESKY FACTORISATION OF G. THE TRANSPOSE */
00879 /*     OF THE FACTOR WILL BE PLACED IN THE R-PARTITION OF W. */
00880 
00881  L90:
00882     ir = iwr;
00883     i__1 = *n;
00884     for (j = 1; j <= i__1; ++j) {
00885         ira = iwr;
00886         irb = ir + 1;
00887         i__2 = j;
00888         for (i = 1; i <= i__2; ++i) {
00889             temp = g[i + j * g_dim1];
00890             if (i == 1) {
00891                 goto L110;
00892             }
00893             i__3 = ir;
00894             for (k = irb; k <= i__3; ++k) {
00895                 ++ira;
00896 /* L100: */
00897                 temp -= w[k] * w[ira];
00898             }
00899         L110:
00900             ++ir;
00901             ++ira;
00902             if (i < j) {
00903                 w[ir] = temp / w[ira];
00904             }
00905 /* L120: */
00906         }
00907         if (temp < *vsmall) {
00908             goto L140;
00909         }
00910 /* L130: */
00911         w[ir] = sqrt(temp);
00912     }
00913     goto L170;
00914 
00915 /*     INCREASE FURTHER THE DIAGONAL ELEMENT OF G. */
00916 
00917  L140:
00918     w[j] = one;
00919     sumx = one;
00920     k = j;
00921  L150:
00922     sum = zero;
00923     ira = ir - 1;
00924     i__1 = j;
00925     for (i = k; i <= i__1; ++i) {
00926         sum -= w[ira] * w[i];
00927 /* L160: */
00928         ira += i;
00929     }
00930     ir -= k;
00931     --k;
00932     w[k] = sum / w[ir];
00933 /* Computing 2nd power */
00934     d__1 = w[k];
00935     sumx += d__1 * d__1;
00936     if (k >= 2) {
00937         goto L150;
00938     }
00939     *diag = *diag + *vsmall - temp / sumx;
00940     goto L70;
00941 
00942 /*     STORE THE CHOLESKY FACTORISATION IN THE R-PARTITION */
00943 /*     OF W. */
00944 
00945  L165:
00946     ir = iwr;
00947     i__1 = *n;
00948     for (i = 1; i <= i__1; ++i) {
00949         i__2 = i;
00950         for (j = 1; j <= i__2; ++j) {
00951             ++ir;
00952 /* L166: */
00953             w[ir] = g[j + i * g_dim1];
00954         }
00955     }
00956 
00957 /*     SET Z THE INVERSE OF THE MATRIX IN R. */
00958 
00959  L170:
00960     nm = *n - 1;
00961     i__2 = *n;
00962     for (i = 1; i <= i__2; ++i) {
00963         iz = iwz + i;
00964         if (i == 1) {
00965             goto L190;
00966         }
00967         i__1 = i;
00968         for (j = 2; j <= i__1; ++j) {
00969             w[iz] = zero;
00970 /* L180: */
00971             iz += *n;
00972         }
00973     L190:
00974         ir = iwr + (i + i * i) / 2;
00975         w[iz] = one / w[ir];
00976         if (i == *n) {
00977             goto L220;
00978         }
00979         iza = iz;
00980         i__1 = nm;
00981         for (j = i; j <= i__1; ++j) {
00982             ir += i;
00983             sum = zero;
00984             i__3 = iz;
00985             i__4 = *n;
00986             for (k = iza; i__4 < 0 ? k >= i__3 : k <= i__3; k += i__4) {
00987                 sum += w[k] * w[ir];
00988 /* L200: */
00989                 ++ir;
00990             }
00991             iz += *n;
00992 /* L210: */
00993             w[iz] = -sum / w[ir];
00994         }
00995     L220:
00996         ;
00997     }
00998 
00999 /*     SET THE INITIAL VALUES OF SOME VARIABLES. */
01000 /*     ITERC COUNTS THE NUMBER OF ITERATIONS. */
01001 /*     ITREF IS SET TO ONE WHEN ITERATIVE REFINEMENT IS REQUIRED. */
01002 /*     JFINC INDICATES WHEN TO TEST FOR AN INCREASE IN F. */
01003 
01004     iterc = 1;
01005     itref = 0;
01006     jfinc = -kfinc;
01007 
01008 /*     SET X TO ZERO AND SET THE CORRESPONDING RESIDUALS OF THE */
01009 /*     KUHN-TUCKER CONDITIONS. */
01010 
01011  L230:
01012     iflag = 1;
01013     iws = iww - *n;
01014     i__2 = *n;
01015     for (i = 1; i <= i__2; ++i) {
01016         x[i] = zero;
01017         iw = iww + i;
01018         w[iw] = grad[i];
01019         if (i > *nact) {
01020             goto L240;
01021         }
01022         w[i] = zero;
01023         is = iws + i;
01024         k = iact[i];
01025         if (k <= *m) {
01026             goto L235;
01027         }
01028         if (k > *mn) {
01029             goto L234;
01030         }
01031         k1 = k - *m;
01032         w[is] = xl[k1];
01033         goto L240;
01034     L234:
01035         k1 = k - *mn;
01036         w[is] = -xu[k1];
01037         goto L240;
01038     L235:
01039         w[is] = b[k];
01040     L240:
01041         ;
01042     }
01043     xmag = zero;
01044     vfact = 1.;
01045     if (*nact <= 0) {
01046         goto L340;
01047     } else {
01048         goto L280;
01049     }
01050 
01051 /*     SET THE RESIDUALS OF THE KUHN-TUCKER CONDITIONS FOR GENERAL X. */
01052 
01053  L250:
01054     iflag = 2;
01055     iws = iww - *n;
01056     i__2 = *n;
01057     for (i = 1; i <= i__2; ++i) {
01058         iw = iww + i;
01059         w[iw] = grad[i];
01060         if (*lql) {
01061             goto L259;
01062         }
01063         id = iwd + i;
01064         w[id] = zero;
01065         i__1 = *n;
01066         for (j = i; j <= i__1; ++j) {
01067 /* L251: */
01068             w[id] += g[i + j * g_dim1] * x[j];
01069         }
01070         i__1 = i;
01071         for (j = 1; j <= i__1; ++j) {
01072             id = iwd + j;
01073 /* L252: */
01074             w[iw] += g[j + i * g_dim1] * w[id];
01075         }
01076         goto L260;
01077     L259:
01078         i__1 = *n;
01079         for (j = 1; j <= i__1; ++j) {
01080 /* L261: */
01081             w[iw] += g[i + j * g_dim1] * x[j];
01082         }
01083     L260:
01084         ;
01085     }
01086     if (*nact == 0) {
01087         goto L340;
01088     }
01089     i__2 = *nact;
01090     for (k = 1; k <= i__2; ++k) {
01091         kk = iact[k];
01092         is = iws + k;
01093         if (kk > *m) {
01094             goto L265;
01095         }
01096         w[is] = b[kk];
01097         i__1 = *n;
01098         for (i = 1; i <= i__1; ++i) {
01099             iw = iww + i;
01100             w[iw] -= w[k] * a[kk + i * a_dim1];
01101 /* L264: */
01102             w[is] -= x[i] * a[kk + i * a_dim1];
01103         }
01104         goto L270;
01105     L265:
01106         if (kk > *mn) {
01107             goto L266;
01108         }
01109         k1 = kk - *m;
01110         iw = iww + k1;
01111         w[iw] -= w[k];
01112         w[is] = xl[k1] - x[k1];
01113         goto L270;
01114     L266:
01115         k1 = kk - *mn;
01116         iw = iww + k1;
01117         w[iw] += w[k];
01118         w[is] = -xu[k1] + x[k1];
01119     L270:
01120         ;
01121     }
01122 
01123 /*     PRE-MULTIPLY THE VECTOR IN THE S-PARTITION OF W BY THE */
01124 /*     INVERS OF R TRANSPOSE. */
01125 
01126  L280:
01127     ir = iwr;
01128     ip = iww + 1;
01129     ipp = iww + *n;
01130     il = iws + 1;
01131     iu = iws + *nact;
01132     i__2 = iu;
01133     for (i = il; i <= i__2; ++i) {
01134         sum = zero;
01135         if (i == il) {
01136             goto L300;
01137         }
01138         ju = i - 1;
01139         i__1 = ju;
01140         for (j = il; j <= i__1; ++j) {
01141             ++ir;
01142 /* L290: */
01143             sum += w[ir] * w[j];
01144         }
01145     L300:
01146         ++ir;
01147 /* L310: */
01148         w[i] = (w[i] - sum) / w[ir];
01149     }
01150 
01151 /*     SHIFT X TO SATISFY THE ACTIVE CONSTRAINTS AND MAKE THE */
01152 /*     CORRESPONDING CHANGE TO THE GRADIENT RESIDUALS. */
01153 
01154     i__2 = *n;
01155     for (i = 1; i <= i__2; ++i) {
01156         iz = iwz + i;
01157         sum = zero;
01158         i__1 = iu;
01159         for (j = il; j <= i__1; ++j) {
01160             sum += w[j] * w[iz];
01161 /* L320: */
01162             iz += *n;
01163         }
01164         x[i] += sum;
01165         if (*lql) {
01166             goto L329;
01167         }
01168         id = iwd + i;
01169         w[id] = zero;
01170         i__1 = *n;
01171         for (j = i; j <= i__1; ++j) {
01172 /* L321: */
01173             w[id] += g[i + j * g_dim1] * sum;
01174         }
01175         iw = iww + i;
01176         i__1 = i;
01177         for (j = 1; j <= i__1; ++j) {
01178             id = iwd + j;
01179 /* L322: */
01180             w[iw] += g[j + i * g_dim1] * w[id];
01181         }
01182         goto L330;
01183     L329:
01184         i__1 = *n;
01185         for (j = 1; j <= i__1; ++j) {
01186             iw = iww + j;
01187 /* L331: */
01188             w[iw] += sum * g[i + j * g_dim1];
01189         }
01190     L330:
01191         ;
01192     }
01193 
01194 /*     FORM THE SCALAR PRODUCT OF THE CURRENT GRADIENT RESIDUALS */
01195 /*     WITH EACH COLUMN OF Z. */
01196 
01197  L340:
01198     kflag = 1;
01199     goto L930;
01200  L350:
01201     if (*nact == *n) {
01202         goto L380;
01203     }
01204 
01205 /*     SHIFT X SO THAT IT SATISFIES THE REMAINING KUHN-TUCKER */
01206 /*     CONDITIONS. */
01207 
01208     il = iws + *nact + 1;
01209     iza = iwz + *nact * *n;
01210     i__2 = *n;
01211     for (i = 1; i <= i__2; ++i) {
01212         sum = zero;
01213         iz = iza + i;
01214         i__1 = iww;
01215         for (j = il; j <= i__1; ++j) {
01216             sum += w[iz] * w[j];
01217 /* L360: */
01218             iz += *n;
01219         }
01220 /* L370: */
01221         x[i] -= sum;
01222     }
01223     *info = 0;
01224     if (*nact == 0) {
01225         goto L410;
01226     }
01227 
01228 /*     UPDATE THE LAGRANGE MULTIPLIERS. */
01229 
01230  L380:
01231     lflag = 3;
01232     goto L740;
01233  L390:
01234     i__2 = *nact;
01235     for (k = 1; k <= i__2; ++k) {
01236         iw = iww + k;
01237 /* L400: */
01238         w[k] += w[iw];
01239     }
01240 
01241 /*     REVISE THE VALUES OF XMAG. */
01242 /*     BRANCH IF ITERATIVE REFINEMENT IS REQUIRED. */
01243 
01244  L410:
01245     jflag = 1;
01246     goto L910;
01247  L420:
01248     if (iflag == itref) {
01249         goto L250;
01250     }
01251 
01252 /*     DELETE A CONSTRAINT IF A LAGRANGE MULTIPLIER OF AN */
01253 /*     INEQUALITY CONSTRAINT IS NEGATIVE. */
01254 
01255     kdrop = 0;
01256     goto L440;
01257  L430:
01258     ++kdrop;
01259     if (w[kdrop] >= zero) {
01260         goto L440;
01261     }
01262     if (iact[kdrop] <= *meq) {
01263         goto L440;
01264     }
01265     nu = *nact;
01266     mflag = 1;
01267     goto L800;
01268  L440:
01269     if (kdrop < *nact) {
01270         goto L430;
01271     }
01272 
01273 /*     SEEK THE GREATEAST NORMALISED CONSTRAINT VIOLATION, DISREGARDING */
01274 
01275 /*     ANY THAT MAY BE DUE TO COMPUTER ROUNDING ERRORS. */
01276 
01277  L450:
01278     cvmax = zero;
01279     if (*m <= 0) {
01280         goto L481;
01281     }
01282     i__2 = *m;
01283     for (k = 1; k <= i__2; ++k) {
01284         ia = iwa + k;
01285         if (w[ia] <= zero) {
01286             goto L480;
01287         }
01288         sum = -b[k];
01289         i__1 = *n;
01290         for (i = 1; i <= i__1; ++i) {
01291 /* L460: */
01292             sum += x[i] * a[k + i * a_dim1];
01293         }
01294         sumx = -sum * w[ia];
01295         if (k <= *meq) {
01296             sumx = abs(sumx);
01297         }
01298         if (sumx <= cvmax) {
01299             goto L480;
01300         }
01301         temp = (d__1 = b[k], abs(d__1));
01302         i__1 = *n;
01303         for (i = 1; i <= i__1; ++i) {
01304 /* L470: */
01305             temp += (d__1 = x[i] * a[k + i * a_dim1], abs(d__1));
01306         }
01307         tempa = temp + abs(sum);
01308         if (tempa <= temp) {
01309             goto L480;
01310         }
01311         temp += onha * abs(sum);
01312         if (temp <= tempa) {
01313             goto L480;
01314         }
01315         cvmax = sumx;
01316         res = sum;
01317         knext = k;
01318     L480:
01319         ;
01320     }
01321  L481:
01322     i__2 = *n;
01323     for (k = 1; k <= i__2; ++k) {
01324         lower = TRUE_;
01325         ia = iwa + *m + k;
01326         if (w[ia] <= zero) {
01327             goto L485;
01328         }
01329         sum = xl[k] - x[k];
01330         if (sum < 0.) {
01331             goto L482;
01332         } else if (sum == 0) {
01333             goto L485;
01334         } else {
01335             goto L483;
01336         }
01337     L482:
01338         sum = x[k] - xu[k];
01339         lower = FALSE_;
01340     L483:
01341         if (sum <= cvmax) {
01342             goto L485;
01343         }
01344         cvmax = sum;
01345         res = -sum;
01346         knext = k + *m;
01347         if (lower) {
01348             goto L485;
01349         }
01350         knext = k + *mn;
01351     L485:
01352         ;
01353     }
01354 
01355 /*     TEST FOR CONVERGENCE */
01356 
01357     *info = 0;
01358     if (cvmax <= *vsmall) {
01359         goto L700;
01360     }
01361 
01362 /*     RETURN IF, DUE TO ROUNDING ERRORS, THE ACTUAL CHANGE IN */
01363 /*     X MAY NOT INCREASE THE OBJECTIVE FUNCTION */
01364 
01365     ++jfinc;
01366     if (jfinc == 0) {
01367         goto L510;
01368     }
01369     if (jfinc != ifinc) {
01370         goto L530;
01371     }
01372     fdiff = zero;
01373     fdiffa = zero;
01374     i__2 = *n;
01375     for (i = 1; i <= i__2; ++i) {
01376         sum = two * grad[i];
01377         sumx = abs(sum);
01378         if (*lql) {
01379             goto L489;
01380         }
01381         id = iwd + i;
01382         w[id] = zero;
01383         i__1 = *n;
01384         for (j = i; j <= i__1; ++j) {
01385             ix = iwx + j;
01386 /* L486: */
01387             w[id] += g[i + j * g_dim1] * (w[ix] + x[j]);
01388         }
01389         i__1 = i;
01390         for (j = 1; j <= i__1; ++j) {
01391             id = iwd + j;
01392             temp = g[j + i * g_dim1] * w[id];
01393             sum += temp;
01394 /* L487: */
01395             sumx += abs(temp);
01396         }
01397         goto L495;
01398     L489:
01399         i__1 = *n;
01400         for (j = 1; j <= i__1; ++j) {
01401             ix = iwx + j;
01402             temp = g[i + j * g_dim1] * (w[ix] + x[j]);
01403             sum += temp;
01404 /* L490: */
01405             sumx += abs(temp);
01406         }
01407     L495:
01408         ix = iwx + i;
01409         fdiff += sum * (x[i] - w[ix]);
01410 /* L500: */
01411         fdiffa += sumx * (d__1 = x[i] - w[ix], abs(d__1));
01412     }
01413     *info = 2;
01414     sum = fdiffa + fdiff;
01415     if (sum <= fdiffa) {
01416         goto L700;
01417     }
01418     temp = fdiffa + onha * fdiff;
01419     if (temp <= sum) {
01420         goto L700;
01421     }
01422     jfinc = 0;
01423     *info = 0;
01424  L510:
01425     i__2 = *n;
01426     for (i = 1; i <= i__2; ++i) {
01427         ix = iwx + i;
01428 /* L520: */
01429         w[ix] = x[i];
01430     }
01431 
01432 /*     FORM THE SCALAR PRODUCT OF THE NEW CONSTRAINT NORMAL WITH EACH */
01433 /*     COLUMN OF Z. PARNEW WILL BECOME THE LAGRANGE MULTIPLIER OF */
01434 /*     THE NEW CONSTRAINT. */
01435 
01436  L530:
01437     ++iterc;
01438     if (iterc <= *maxit) {
01439         goto L531;
01440     }
01441     *info = 1;
01442     goto L710;
01443  L531:
01444     iws = iwr + (*nact + *nact * *nact) / 2;
01445     if (knext > *m) {
01446         goto L541;
01447     }
01448     i__2 = *n;
01449     for (i = 1; i <= i__2; ++i) {
01450         iw = iww + i;
01451 /* L540: */
01452         w[iw] = a[knext + i * a_dim1];
01453     }
01454     goto L549;
01455  L541:
01456     i__2 = *n;
01457     for (i = 1; i <= i__2; ++i) {
01458         iw = iww + i;
01459 /* L542: */
01460         w[iw] = zero;
01461     }
01462     k1 = knext - *m;
01463     if (k1 > *n) {
01464         goto L545;
01465     }
01466     iw = iww + k1;
01467     w[iw] = one;
01468     iz = iwz + k1;
01469     i__2 = *n;
01470     for (i = 1; i <= i__2; ++i) {
01471         is = iws + i;
01472         w[is] = w[iz];
01473 /* L543: */
01474         iz += *n;
01475     }
01476     goto L550;
01477  L545:
01478     k1 = knext - *mn;
01479     iw = iww + k1;
01480     w[iw] = -one;
01481     iz = iwz + k1;
01482     i__2 = *n;
01483     for (i = 1; i <= i__2; ++i) {
01484         is = iws + i;
01485         w[is] = -w[iz];
01486 /* L546: */
01487         iz += *n;
01488     }
01489     goto L550;
01490  L549:
01491     kflag = 2;
01492     goto L930;
01493  L550:
01494     parnew = zero;
01495 
01496 /*     APPLY GIVENS ROTATIONS TO MAKE THE LAST (N-NACT-2) SCALAR */
01497 /*     PRODUCTS EQUAL TO ZERO. */
01498 
01499     if (*nact == *n) {
01500         goto L570;
01501     }
01502     nu = *n;
01503     nflag = 1;
01504     goto L860;
01505 
01506 /*     BRANCH IF THERE IS NO NEED TO DELETE A CONSTRAINT. */
01507 
01508  L560:
01509     is = iws + *nact;
01510     if (*nact == 0) {
01511         goto L640;
01512     }
01513     suma = zero;
01514     sumb = zero;
01515     sumc = zero;
01516     iz = iwz + *nact * *n;
01517     i__2 = *n;
01518     for (i = 1; i <= i__2; ++i) {
01519         ++iz;
01520         iw = iww + i;
01521         suma += w[iw] * w[iz];
01522         sumb += (d__1 = w[iw] * w[iz], abs(d__1));
01523 /* L563: */
01524 /* Computing 2nd power */
01525         d__1 = w[iz];
01526         sumc += d__1 * d__1;
01527     }
01528     temp = sumb + abs(suma) * .1;
01529     tempa = sumb + abs(suma) * .2;
01530     if (temp <= sumb) {
01531         goto L570;
01532     }
01533     if (tempa <= temp) {
01534         goto L570;
01535     }
01536     if (sumb > *vsmall) {
01537         goto L5;
01538     }
01539     goto L570;
01540  L5:
01541     sumc = sqrt(sumc);
01542     ia = iwa + knext;
01543     if (knext <= *m) {
01544         sumc /= w[ia];
01545     }
01546     temp = sumc + abs(suma) * .1;
01547     tempa = sumc + abs(suma) * .2;
01548     if (temp <= sumc) {
01549         goto L567;
01550     }
01551     if (tempa <= temp) {
01552         goto L567;
01553     }
01554     goto L640;
01555 
01556 /*     CALCULATE THE MULTIPLIERS FOR THE NEW CONSTRAINT NORMAL */
01557 /*     EXPRESSED IN TERMS OF THE ACTIVE CONSTRAINT NORMALS. */
01558 /*     THEN WORK OUT WHICH CONTRAINT TO DROP. */
01559 
01560  L567:
01561     lflag = 4;
01562     goto L740;
01563  L570:
01564     lflag = 1;
01565     goto L740;
01566 
01567 /*     COMPLETE THE TEST FOR LINEARLY DEPENDENT CONSTRAINTS. */
01568 
01569  L571:
01570     if (knext > *m) {
01571         goto L574;
01572     }
01573     i__2 = *n;
01574     for (i = 1; i <= i__2; ++i) {
01575         suma = a[knext + i * a_dim1];
01576         sumb = abs(suma);
01577         if (*nact == 0) {
01578             goto L581;
01579         }
01580         i__1 = *nact;
01581         for (k = 1; k <= i__1; ++k) {
01582             kk = iact[k];
01583             if (kk <= *m) {
01584                 goto L568;
01585             }
01586             kk -= *m;
01587             temp = zero;
01588             if (kk == i) {
01589                 temp = w[iww + kk];
01590             }
01591             kk -= *n;
01592             if (kk == i) {
01593                 temp = -w[iww + kk];
01594             }
01595             goto L569;
01596         L568:
01597             iw = iww + k;
01598             temp = w[iw] * a[kk + i * a_dim1];
01599         L569:
01600             suma -= temp;
01601 /* L572: */
01602             sumb += abs(temp);
01603         }
01604     L581:
01605         if (suma <= *vsmall) {
01606             goto L573;
01607         }
01608         temp = sumb + abs(suma) * .1;
01609         tempa = sumb + abs(suma) * .2;
01610         if (temp <= sumb) {
01611             goto L573;
01612         }
01613         if (tempa <= temp) {
01614             goto L573;
01615         }
01616         goto L630;
01617     L573:
01618         ;
01619     }
01620     lflag = 1;
01621     goto L775;
01622  L574:
01623     k1 = knext - *m;
01624     if (k1 > *n) {
01625         k1 -= *n;
01626     }
01627     i__2 = *n;
01628     for (i = 1; i <= i__2; ++i) {
01629         suma = zero;
01630         if (i != k1) {
01631             goto L575;
01632         }
01633         suma = one;
01634         if (knext > *mn) {
01635             suma = -one;
01636         }
01637     L575:
01638         sumb = abs(suma);
01639         if (*nact == 0) {
01640             goto L582;
01641         }
01642         i__1 = *nact;
01643         for (k = 1; k <= i__1; ++k) {
01644             kk = iact[k];
01645             if (kk <= *m) {
01646                 goto L579;
01647             }
01648             kk -= *m;
01649             temp = zero;
01650             if (kk == i) {
01651                 temp = w[iww + kk];
01652             }
01653             kk -= *n;
01654             if (kk == i) {
01655                 temp = -w[iww + kk];
01656             }
01657             goto L576;
01658         L579:
01659             iw = iww + k;
01660             temp = w[iw] * a[kk + i * a_dim1];
01661         L576:
01662             suma -= temp;
01663 /* L577: */
01664             sumb += abs(temp);
01665         }
01666     L582:
01667         temp = sumb + abs(suma) * .1;
01668         tempa = sumb + abs(suma) * .2;
01669         if (temp <= sumb) {
01670             goto L578;
01671         }
01672         if (tempa <= temp) {
01673             goto L578;
01674         }
01675         goto L630;
01676     L578:
01677         ;
01678     }
01679     lflag = 1;
01680     goto L775;
01681 
01682 /*     BRANCH IF THE CONTRAINTS ARE INCONSISTENT. */
01683 
01684  L580:
01685     *info = -knext;
01686     if (kdrop == 0) {
01687         goto L700;
01688     }
01689     parinc = ratio;
01690     parnew = parinc;
01691 
01692 /*     REVISE THE LAGRANGE MULTIPLIERS OF THE ACTIVE CONSTRAINTS. */
01693 
01694  L590:
01695     if (*nact == 0) {
01696         goto L601;
01697     }
01698     i__2 = *nact;
01699     for (k = 1; k <= i__2; ++k) {
01700         iw = iww + k;
01701         w[k] -= parinc * w[iw];
01702         if (iact[k] > *meq) {
01703 /* Computing MAX */
01704             d__1 = zero, d__2 = w[k];
01705             w[k] = max(d__1,d__2);
01706         }
01707 /* L600: */
01708     }
01709  L601:
01710     if (kdrop == 0) {
01711         goto L680;
01712     }
01713 
01714 /*     DELETE THE CONSTRAINT TO BE DROPPED. */
01715 /*     SHIFT THE VECTOR OF SCALAR PRODUCTS. */
01716 /*     THEN, IF APPROPRIATE, MAKE ONE MORE SCALAR PRODUCT ZERO. */
01717 
01718     nu = *nact + 1;
01719     mflag = 2;
01720     goto L800;
01721  L610:
01722     iws = iws - *nact - 1;
01723     nu = min(*n,nu);
01724     i__2 = nu;
01725     for (i = 1; i <= i__2; ++i) {
01726         is = iws + i;
01727         j = is + *nact;
01728 /* L620: */
01729         w[is] = w[j + 1];
01730     }
01731     nflag = 2;
01732     goto L860;
01733 
01734 /*     CALCULATE THE STEP TO THE VIOLATED CONSTRAINT. */
01735 
01736  L630:
01737     is = iws + *nact;
01738  L640:
01739     sumy = w[is + 1];
01740     step = -res / sumy;
01741     parinc = step / sumy;
01742     if (*nact == 0) {
01743         goto L660;
01744     }
01745 
01746 /*     CALCULATE THE CHANGES TO THE LAGRANGE MULTIPLIERS, AND REDUCE */
01747 /*     THE STEP ALONG THE NEW SEARCH DIRECTION IF NECESSARY. */
01748 
01749     lflag = 2;
01750     goto L740;
01751  L650:
01752     if (kdrop == 0) {
01753         goto L660;
01754     }
01755     temp = one - ratio / parinc;
01756     if (temp <= zero) {
01757         kdrop = 0;
01758     }
01759     if (kdrop == 0) {
01760         goto L660;
01761     }
01762     step = ratio * sumy;
01763     parinc = ratio;
01764     res = temp * res;
01765 
01766 /*     UPDATE X AND THE LAGRANGE MULTIPIERS. */
01767 /*     DROP A CONSTRAINT IF THE FULL STEP IS NOT TAKEN. */
01768 
01769  L660:
01770     iwy = iwz + *nact * *n;
01771     i__2 = *n;
01772     for (i = 1; i <= i__2; ++i) {
01773         iy = iwy + i;
01774 /* L670: */
01775         x[i] += step * w[iy];
01776     }
01777     parnew += parinc;
01778     if (*nact >= 1) {
01779         goto L590;
01780     }
01781 
01782 /*     ADD THE NEW CONSTRAINT TO THE ACTIVE SET. */
01783 
01784  L680:
01785     ++(*nact);
01786     w[*nact] = parnew;
01787     iact[*nact] = knext;
01788     ia = iwa + knext;
01789     if (knext > *mn) {
01790         ia -= *n;
01791     }
01792     w[ia] = -w[ia];
01793 
01794 /*     ESTIMATE THE MAGNITUDE OF X. THEN BEGIN A NEW ITERATION, */
01795 /*     RE-INITILISING X IF THIS MAGNITUDE IS SMALL. */
01796 
01797     jflag = 2;
01798     goto L910;
01799  L690:
01800     if (sum < xmagr * xmag) {
01801         goto L230;
01802     }
01803     if (itref <= 0) {
01804         goto L450;
01805     } else {
01806         goto L250;
01807     }
01808 
01809 /*     INITIATE ITERATIVE REFINEMENT IF IT HAS NOT YET BEEN USED, */
01810 /*     OR RETURN AFTER RESTORING THE DIAGONAL ELEMENTS OF G. */
01811 
01812  L700:
01813     if (iterc == 0) {
01814         goto L710;
01815     }
01816     ++itref;
01817     jfinc = -1;
01818     if (itref == 1) {
01819         goto L250;
01820     }
01821  L710:
01822     if (! (*lql)) {
01823         return 0;
01824     }
01825     i__2 = *n;
01826     for (i = 1; i <= i__2; ++i) {
01827         id = iwd + i;
01828 /* L720: */
01829         g[i + i * g_dim1] = w[id];
01830     }
01831  L730:
01832     return 0;
01833 
01834 
01835 /*     THE REMAINIG INSTRUCTIONS ARE USED AS SUBROUTINES. */
01836 
01837 
01838 /* ******************************************************************** */
01839 
01840 
01841 
01842 /*     CALCULATE THE LAGRANGE MULTIPLIERS BY PRE-MULTIPLYING THE */
01843 /*     VECTOR IN THE S-PARTITION OF W BY THE INVERSE OF R. */
01844 
01845  L740:
01846     ir = iwr + (*nact + *nact * *nact) / 2;
01847     i = *nact;
01848     sum = zero;
01849     goto L770;
01850  L750:
01851     ira = ir - 1;
01852     sum = zero;
01853     if (*nact == 0) {
01854         goto L761;
01855     }
01856     i__2 = *nact;
01857     for (j = i; j <= i__2; ++j) {
01858         iw = iww + j;
01859         sum += w[ira] * w[iw];
01860 /* L760: */
01861         ira += j;
01862     }
01863  L761:
01864     ir -= i;
01865     --i;
01866  L770:
01867     iw = iww + i;
01868     is = iws + i;
01869     w[iw] = (w[is] - sum) / w[ir];
01870     if (i > 1) {
01871         goto L750;
01872     }
01873     if (lflag == 3) {
01874         goto L390;
01875     }
01876     if (lflag == 4) {
01877         goto L571;
01878     }
01879 
01880 /*     CALCULATE THE NEXT CONSTRAINT TO DROP. */
01881 
01882  L775:
01883     ip = iww + 1;
01884     ipp = iww + *nact;
01885     kdrop = 0;
01886     if (*nact == 0) {
01887         goto L791;
01888     }
01889     i__2 = *nact;
01890     for (k = 1; k <= i__2; ++k) {
01891         if (iact[k] <= *meq) {
01892             goto L790;
01893         }
01894         iw = iww + k;
01895         if (res * w[iw] >= zero) {
01896             goto L790;
01897         }
01898         temp = w[k] / w[iw];
01899         if (kdrop == 0) {
01900             goto L780;
01901         }
01902         if (abs(temp) >= abs(ratio)) {
01903             goto L790;
01904         }
01905     L780:
01906         kdrop = k;
01907         ratio = temp;
01908     L790:
01909         ;
01910     }
01911  L791:
01912     switch ((int)lflag) {
01913     case 1:  goto L580;
01914     case 2:  goto L650;
01915     }
01916 
01917 
01918 /* ******************************************************************** */
01919 
01920 
01921 
01922 /*     DROP THE CONSTRAINT IN POSITION KDROP IN THE ACTIVE SET. */
01923 
01924  L800:
01925     ia = iwa + iact[kdrop];
01926     if (iact[kdrop] > *mn) {
01927         ia -= *n;
01928     }
01929     w[ia] = -w[ia];
01930     if (kdrop == *nact) {
01931         goto L850;
01932     }
01933 
01934 /*     SET SOME INDICES AND CALCULATE THE ELEMENTS OF THE NEXT */
01935 /*     GIVENS ROTATION. */
01936 
01937     iz = iwz + kdrop * *n;
01938     ir = iwr + (kdrop + kdrop * kdrop) / 2;
01939  L810:
01940     ira = ir;
01941     ir = ir + kdrop + 1;
01942 /* Computing MAX */
01943     d__3 = (d__1 = w[ir - 1], abs(d__1)), d__4 = (d__2 = w[ir], abs(d__2));
01944     temp = max(d__3,d__4);
01945 /* Computing 2nd power */
01946     d__1 = w[ir - 1] / temp;
01947 /* Computing 2nd power */
01948     d__2 = w[ir] / temp;
01949     sum = temp * sqrt(d__1 * d__1 + d__2 * d__2);
01950     ga = w[ir - 1] / sum;
01951     gb = w[ir] / sum;
01952 
01953 /*     EXCHANGE THE COLUMNS OF R. */
01954 
01955     i__2 = kdrop;
01956     for (i = 1; i <= i__2; ++i) {
01957         ++ira;
01958         j = ira - kdrop;
01959         temp = w[ira];
01960         w[ira] = w[j];
01961 /* L820: */
01962         w[j] = temp;
01963     }
01964     w[ir] = zero;
01965 
01966 /*     APPLY THE ROTATION TO THE ROWS OF R. */
01967 
01968     w[j] = sum;
01969     ++kdrop;
01970     i__2 = nu;
01971     for (i = kdrop; i <= i__2; ++i) {
01972         temp = ga * w[ira] + gb * w[ira + 1];
01973         w[ira + 1] = ga * w[ira + 1] - gb * w[ira];
01974         w[ira] = temp;
01975 /* L830: */
01976         ira += i;
01977     }
01978 
01979 /*     APPLY THE ROTATION TO THE COLUMNS OF Z. */
01980 
01981     i__2 = *n;
01982     for (i = 1; i <= i__2; ++i) {
01983         ++iz;
01984         j = iz - *n;
01985         temp = ga * w[j] + gb * w[iz];
01986         w[iz] = ga * w[iz] - gb * w[j];
01987 /* L840: */
01988         w[j] = temp;
01989     }
01990 
01991 /*     REVISE IACT AND THE LAGRANGE MULTIPLIERS. */
01992 
01993     iact[kdrop - 1] = iact[kdrop];
01994     w[kdrop - 1] = w[kdrop];
01995     if (kdrop < *nact) {
01996         goto L810;
01997     }
01998  L850:
01999     --(*nact);
02000     switch ((int)mflag) {
02001     case 1:  goto L250;
02002     case 2:  goto L610;
02003     }
02004 
02005 
02006 /* ******************************************************************** */
02007 
02008 
02009 
02010 /*     APPLY GIVENS ROTATION TO REDUCE SOME OF THE SCALAR */
02011 /*     PRODUCTS IN THE S-PARTITION OF W TO ZERO. */
02012 
02013  L860:
02014     iz = iwz + nu * *n;
02015  L870:
02016     iz -= *n;
02017  L880:
02018     is = iws + nu;
02019     --nu;
02020     if (nu == *nact) {
02021         goto L900;
02022     }
02023     if (w[is] == zero) {
02024         goto L870;
02025     }
02026 /* Computing MAX */
02027     d__3 = (d__1 = w[is - 1], abs(d__1)), d__4 = (d__2 = w[is], abs(d__2));
02028     temp = max(d__3,d__4);
02029 /* Computing 2nd power */
02030     d__1 = w[is - 1] / temp;
02031 /* Computing 2nd power */
02032     d__2 = w[is] / temp;
02033     sum = temp * sqrt(d__1 * d__1 + d__2 * d__2);
02034     ga = w[is - 1] / sum;
02035     gb = w[is] / sum;
02036     w[is - 1] = sum;
02037     i__2 = *n;
02038     for (i = 1; i <= i__2; ++i) {
02039         k = iz + *n;
02040         temp = ga * w[iz] + gb * w[k];
02041         w[k] = ga * w[k] - gb * w[iz];
02042         w[iz] = temp;
02043 /* L890: */
02044         --iz;
02045     }
02046     goto L880;
02047  L900:
02048     switch ((int)nflag) {
02049     case 1:  goto L560;
02050     case 2:  goto L630;
02051     }
02052 
02053 
02054 /* ******************************************************************** */
02055 
02056 
02057 
02058 /*     CALCULATE THE MAGNITUDE OF X AN REVISE XMAG. */
02059 
02060  L910:
02061     sum = zero;
02062     i__2 = *n;
02063     for (i = 1; i <= i__2; ++i) {
02064         sum += (d__1 = x[i], abs(d__1)) * vfact * ((d__2 = grad[i], abs(d__2))
02065                                                    + (d__3 = g[i + i * g_dim1] * x[i], abs(d__3)));
02066         if (*lql) {
02067             goto L920;
02068         }
02069         if (sum < 1e-30) {
02070             goto L920;
02071         }
02072         vfact *= 1e-10;
02073         sum *= 1e-10;
02074         xmag *= 1e-10;
02075     L920:
02076         ;
02077     }
02078 /* L925: */
02079     xmag = max(xmag,sum);
02080     switch ((int)jflag) {
02081     case 1:  goto L420;
02082     case 2:  goto L690;
02083     }
02084 
02085 
02086 /* ******************************************************************** */
02087 
02088 
02089 
02090 /*     PRE-MULTIPLY THE VECTOR IN THE W-PARTITION OF W BY Z TRANSPOSE. */
02091 
02092  L930:
02093     jl = iww + 1;
02094     iz = iwz;
02095     i__2 = *n;
02096     for (i = 1; i <= i__2; ++i) {
02097         is = iws + i;
02098         w[is] = zero;
02099         iwwn = iww + *n;
02100         i__1 = iwwn;
02101         for (j = jl; j <= i__1; ++j) {
02102             ++iz;
02103 /* L940: */
02104             w[is] += w[iz] * w[j];
02105         }
02106     }
02107     switch ((int)kflag) {
02108     case 1:  goto L350;
02109     case 2:  goto L550;
02110     }
02111     return 0;
02112 } /* ql0002_ */
02113 
02114 #ifdef uNdEfInEd
02115 comments from the converter:  (stderr from f2c)
02116     ql0001:
02117 ql0002:
02118 #endif
02119 
02120 
02121 /*
02122   Local Variables:
02123   mode:c++
02124   c-basic-offset:4
02125   c-file-style:"stroustrup"
02126   c-file-offsets:((innamespace . 0)(inline-open . 0))
02127   indent-tabs-mode:nil
02128   fill-column:79
02129   End:
02130 */
02131 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines