PLearn 0.1
|
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 :