home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / fsqp / ql0002.f < prev    next >
Text File  |  1996-09-28  |  24KB  |  942 lines
  1. C
  2.       SUBROUTINE QL0002(N,M,MEQ,MMAX,MN,MNN,NMAX,LQL,A,B,GRAD,G,
  3.      1      XL,XU,X,NACT,IACT,MAXIT,VSMALL,INFO,DIAG,W,LW)
  4. C
  5. C**************************************************************************
  6. C
  7. C
  8. C   THIS SUBROUTINE SOLVES THE QUADRATIC PROGRAMMING PROBLEM 
  9. C
  10. C       MINIMIZE      GRAD'*X  +  0.5 * X*G*X
  11. C       SUBJECT TO    A(K)*X  =  B(K)   K=1,2,...,MEQ,
  12. C                     A(K)*X >=  B(K)   K=MEQ+1,...,M,
  13. C                     XL  <=  X  <=  XU
  14. C
  15. C   THE QUADRATIC PROGRAMMING METHOD PROCEEDS FROM AN INITIAL CHOLESKY-
  16. C   DECOMPOSITION OF THE OBJECTIVE FUNCTION MATRIX, TO CALCULATE THE
  17. C   UNIQUELY DETERMINED MINIMIZER OF THE UNCONSTRAINED PROBLEM. 
  18. C   SUCCESSIVELY ALL VIOLATED CONSTRAINTS ARE ADDED TO A WORKING SET 
  19. C   AND A MINIMIZER OF THE OBJECTIVE FUNCTION SUBJECT TO ALL CONSTRAINTS 
  20. C   IN THIS WORKING SET IS COMPUTED. IT IS POSSIBLE THAT CONSTRAINTS
  21. C   HAVE TO LEAVE THE WORKING SET.
  22. C
  23. C
  24. C   DESCRIPTION OF PARAMETERS:    
  25. C
  26. C     N        : IS THE NUMBER OF VARIABLES.
  27. C     M        : TOTAL NUMBER OF CONSTRAINTS.
  28. C     MEQ      : NUMBER OF EQUALITY CONTRAINTS.
  29. C     MMAX     : ROW DIMENSION OF A, DIMENSION OF B. MMAX MUST BE AT
  30. C                LEAST ONE AND GREATER OR EQUAL TO M.
  31. C     MN       : MUST BE EQUAL M + N.
  32. C     MNN      : MUST BE EQUAL M + N + N.
  33. C     NMAX     : ROW DIEMSION OF G. MUST BE AT LEAST N.
  34. C     LQL      : DETERMINES INITIAL DECOMPOSITION.
  35. C        LQL = .FALSE.  : THE UPPER TRIANGULAR PART OF THE MATRIX G
  36. C                         CONTAINS INITIALLY THE CHOLESKY-FACTOR OF A SUITABLE 
  37. C                         DECOMPOSITION.
  38. C        LQL = .TRUE.   : THE INITIAL CHOLESKY-FACTORISATION OF G IS TO BE
  39. C                         PERFORMED BY THE ALGORITHM.
  40. C     A(MMAX,NMAX) : A IS A MATRIX WHOSE COLUMNS ARE THE CONSTRAINTS NORMALS.
  41. C     B(MMAX)  : CONTAINS THE RIGHT HAND SIDES OF THE CONSTRAINTS.
  42. C     GRAD(N)  : CONTAINS THE OBJECTIVE FUNCTION VECTOR GRAD.
  43. C     G(NMAX,N): CONTAINS THE SYMMETRIC OBJECTIVE FUNCTION MATRIX.
  44. C     XL(N), XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR X.
  45. C     X(N)     : VECTOR OF VARIABLES.
  46. C     NACT     : FINAL NUMBER OF ACTIVE CONSTRAINTS.
  47. C     IACT(K) (K=1,2,...,NACT): INDICES OF THE FINAL ACTIVE CONSTRAINTS.
  48. C     INFO     : REASON FOR THE RETURN FROM THE SUBROUTINE.
  49. C         INFO = 0 : CALCULATION WAS TERMINATED SUCCESSFULLY.
  50. C         INFO = 1 : MAXIMUM NUMBER OF ITERATIONS ATTAINED.
  51. C         INFO = 2 : ACCURACY IS INSUFFICIENT TO MAINTAIN INCREASING
  52. C                    FUNCTION VALUES.
  53. C         INFO < 0 : THE CONSTRAINT WITH INDEX ABS(INFO) AND THE CON-
  54. C                    STRAINTS WHOSE INDICES ARE IACT(K), K=1,2,...,NACT,
  55. C                    ARE INCONSISTENT.
  56. C     MAXIT    : MAXIMUM NUMBER OF ITERATIONS.
  57. C     VSMALL   : REQUIRED ACCURACY TO BE ACHIEVED (E.G. IN THE ORDER OF THE 
  58. C                MACHINE PRECISION FOR SMALL AND WELL-CONDITIONED PROBLEMS).
  59. C     DIAG     : ON RETURN DIAG IS EQUAL TO THE MULTIPLE OF THE UNIT MATRIX
  60. C                THAT WAS ADDED TO G TO ACHIEVE POSITIVE DEFINITENESS.
  61. C     W(LW)    : THE ELEMENTS OF W(.) ARE USED FOR WORKING SPACE. THE LENGTH
  62. C                OF W MUST NOT BE LESS THAN (1.5*NMAX*NMAX + 10*NMAX + M).
  63. C                WHEN INFO = 0 ON RETURN, THE LAGRANGE MULTIPLIERS OF THE
  64. C                FINAL ACTIVE CONSTRAINTS ARE HELD IN W(K), K=1,2,...,NACT.
  65. C   THE VALUES OF N, M, MEQ, MMAX, MN, MNN AND NMAX AND THE ELEMENTS OF
  66. C   A, B, GRAD AND G ARE NOT ALTERED.
  67. C
  68. C   THE FOLLOWING INTEGERS ARE USED TO PARTITION W:
  69. C     THE FIRST N ELEMENTS OF W HOLD LAGRANGE MULTIPLIER ESTIMATES.
  70. C     W(IWZ+I+(N-1)*J) HOLDS THE MATRIX ELEMENT Z(I,J).
  71. C     W(IWR+I+0.5*J*(J-1)) HOLDS THE UPPER TRIANGULAR MATRIX
  72. C       ELEMENT R(I,J). THE SUBSEQUENT N COMPONENTS OF W MAY BE
  73. C       TREATED AS AN EXTRA COLUMN OF R(.,.).
  74. C     W(IWW-N+I) (I=1,2,...,N) ARE USED FOR TEMPORARY STORAGE.
  75. C     W(IWW+I) (I=1,2,...,N) ARE USED FOR TEMPORARY STORAGE.
  76. C     W(IWD+I) (I=1,2,...,N) HOLDS G(I,I) DURING THE CALCULATION.
  77. C     W(IWX+I) (I=1,2,...,N) HOLDS VARIABLES THAT WILL BE USED TO
  78. C       TEST THAT THE ITERATIONS INCREASE THE OBJECTIVE FUNCTION.
  79. C     W(IWA+K) (K=1,2,...,M) USUALLY HOLDS THE RECIPROCAL OF THE
  80. C       LENGTH OF THE K-TH CONSTRAINT, BUT ITS SIGN INDICATES
  81. C       WHETHER THE CONSTRAINT IS ACTIVE.
  82. C
  83. C   
  84. C   AUTHOR:    K. SCHITTKOWSKI,
  85. C              MATHEMATISCHES INSTITUT,
  86. C              UNIVERSITAET BAYREUTH,
  87. C              8580 BAYREUTH,
  88. C              GERMANY, F.R.
  89. C
  90. C   AUTHOR OF ORIGINAL VERSION:
  91. C              M.J.D. POWELL, DAMTP,
  92. C              UNIVERSITY OF CAMBRIDGE, SILVER STREET
  93. C              CAMBRIDGE,
  94. C              ENGLAND
  95. C
  96. C
  97. C   REFERENCE: M.J.D. POWELL: ZQPCVX, A FORTRAN SUBROUTINE FOR CONVEX
  98. C              PROGRAMMING, REPORT DAMTP/1983/NA17, UNIVERSITY OF
  99. C              CAMBRIDGE, ENGLAND, 1983.
  100. C
  101. C
  102. C   VERSION :  2.0 (MARCH, 1987)
  103. C
  104. C
  105. C*************************************************************************
  106. C
  107.       INTEGER MMAX,NMAX,N,LW
  108.       DIMENSION A(MMAX,NMAX),B(MMAX),GRAD(N),G(NMAX,N),X(N),IACT(N),
  109.      1      W(LW),XL(N),XU(N)
  110.       INTEGER M,MEQ,MN,MNN,NACT,IACT,INFO,MAXIT
  111.       DOUBLE PRECISION CVMAX,DIAG,DIAGR,FDIFF,FDIFFA,GA,GB,PARINC,PARNEW       
  112.      1      ,RATIO,RES,STEP,SUM,SUMX,SUMY,SUMA,SUMB,SUMC,TEMP,TEMPA,
  113.      2       VSMALL,XMAG,XMAGR,ZERO,ONE,TWO,ONHA,VFACT
  114.       DOUBLE PRECISION A,B,G,GRAD,W,X,XL,XU
  115. C
  116. C   INTRINSIC FUNCTIONS:   DMAX1,DSQRT,DABS,DMIN1
  117. C
  118.       INTEGER IWZ,IWR,IWW,IWD,IWA,IFINC,KFINC,K,I,IA,ID,II,IR,IRA,
  119.      1     IRB,J,NM,IZ,IZA,ITERC,ITREF,JFINC,IFLAG,IWS,IS,K1,IW,KK,IP,
  120.      2     IPP,IL,IU,JU,KFLAG,LFLAG,JFLAG,KDROP,NU,MFLAG,KNEXT,IX,IWX,
  121.      3     IWY,IY,JL
  122.       LOGICAL LQL,LOWER
  123. C
  124. C   INITIAL ADDRESSES
  125. C
  126.       IWZ=NMAX
  127.       IWR=IWZ+NMAX*NMAX
  128.       IWW=IWR+(NMAX*(NMAX+3))/2
  129.       IWD=IWW+NMAX
  130.       IWX=IWD+NMAX
  131.       IWA=IWX+NMAX
  132. C
  133. C     SET SOME CONSTANTS.
  134. C
  135.       ZERO=0.D+0
  136.       ONE=1.D+0
  137.       TWO=2.D+0
  138.       ONHA=1.5D+0
  139.       VFACT=1.D+0
  140. C
  141. C     SET SOME PARAMETERS.
  142. C     NUMBER LESS THAN VSMALL ARE ASSUMED TO BE NEGLIGIBLE.
  143. C     THE MULTIPLE OF I THAT IS ADDED TO G IS AT MOST DIAGR TIMES
  144. C       THE LEAST MULTIPLE OF I THAT GIVES POSITIVE DEFINITENESS.
  145. C     X IS RE-INITIALISED IF ITS MAGNITUDE IS REDUCED BY THE
  146. C       FACTOR XMAGR.
  147. C     A CHECK IS MADE FOR AN INCREASE IN F EVERY IFINC ITERATIONS,
  148. C       AFTER KFINC ITERATIONS ARE COMPLETED.
  149. C
  150.       DIAGR=TWO
  151.       XMAGR=1.0D-2
  152.       IFINC=3
  153.       KFINC=MAX0(10,N)
  154. C
  155. C     FIND THE RECIPROCALS OF THE LENGTHS OF THE CONSTRAINT NORMALS.
  156. C     RETURN IF A CONSTRAINT IS INFEASIBLE DUE TO A ZERO NORMAL.
  157. C
  158.       NACT=0
  159.       IF (M .LE. 0) GOTO 45
  160.       DO 40 K=1,M
  161.       SUM=ZERO
  162.       DO 10 I=1,N
  163.    10 SUM=SUM+A(K,I)**2
  164.       IF (SUM .GT. ZERO) GOTO 20
  165.       IF (B(K) .EQ. ZERO) GOTO 30
  166.       INFO=-K
  167.       IF (K .LE. MEQ) GOTO 730
  168.       IF (B(K)) 30,30,730
  169.    20 SUM=ONE/DSQRT(SUM)
  170.    30 IA=IWA+K
  171.    40 W(IA)=SUM
  172.    45 DO 50 K=1,N
  173.       IA=IWA+M+K
  174.    50 W(IA)=ONE
  175. C
  176. C     IF NECESSARY INCREASE THE DIAGONAL ELEMENTS OF G.
  177. C
  178.       IF (.NOT. LQL) GOTO 165
  179.       DIAG=ZERO
  180.       DO 60 I=1,N
  181.       ID=IWD+I
  182.       W(ID)=G(I,I)
  183.       DIAG=DMAX1(DIAG,VSMALL-W(ID))
  184.       IF (I .EQ. N) GOTO 60
  185.       II=I+1
  186.       DO 55 J=II,N
  187.       GA=-DMIN1(W(ID),G(J,J))
  188.       GB=DABS(W(ID)-G(J,J))+DABS(G(I,J))
  189.       IF (GB .GT. ZERO) GA=GA+G(I,J)**2/GB
  190.    55 DIAG=DMAX1(DIAG,GA)
  191.    60 CONTINUE
  192.       IF (DIAG .LE. ZERO) GOTO 90
  193.    70 DIAG=DIAGR*DIAG
  194.       DO 80 I=1,N
  195.       ID=IWD+I
  196.    80 G(I,I)=DIAG+W(ID)
  197. C
  198. C     FORM THE CHOLESKY FACTORISATION OF G. THE TRANSPOSE
  199. C     OF THE FACTOR WILL BE PLACED IN THE R-PARTITION OF W.
  200. C
  201.    90 IR=IWR
  202.       DO 130 J=1,N
  203.       IRA=IWR
  204.       IRB=IR+1
  205.       DO 120 I=1,J
  206.       TEMP=G(I,J)
  207.       IF (I .EQ. 1) GOTO 110
  208.       DO 100 K=IRB,IR
  209.       IRA=IRA+1
  210.   100 TEMP=TEMP-W(K)*W(IRA)
  211.   110 IR=IR+1
  212.       IRA=IRA+1
  213.       IF (I .LT. J) W(IR)=TEMP/W(IRA)
  214.   120 CONTINUE
  215.       IF (TEMP .LT. VSMALL) GOTO 140
  216.   130 W(IR)=DSQRT(TEMP)
  217.       GOTO 170
  218. C
  219. C     INCREASE FURTHER THE DIAGONAL ELEMENT OF G.
  220. C
  221.   140 W(J)=ONE
  222.       SUMX=ONE
  223.       K=J
  224.   150 SUM=ZERO
  225.       IRA=IR-1
  226.       DO 160 I=K,J
  227.       SUM=SUM-W(IRA)*W(I)
  228.   160 IRA=IRA+I
  229.       IR=IR-K
  230.       K=K-1
  231.       W(K)=SUM/W(IR)
  232.       SUMX=SUMX+W(K)**2
  233.       IF (K .GE. 2) GOTO 150
  234.       DIAG=DIAG+VSMALL-TEMP/SUMX
  235.       GOTO 70
  236. C
  237. C     STORE THE CHOLESKY FACTORISATION IN THE R-PARTITION
  238. C     OF W.
  239. C
  240.   165 IR=IWR
  241.       DO 166 I=1,N
  242.       DO 166 J=1,I
  243.       IR=IR+1
  244.   166 W(IR)=G(J,I)
  245. C
  246. C     SET Z THE INVERSE OF THE MATRIX IN R.
  247. C
  248.   170 NM=N-1
  249.       DO 220 I=1,N
  250.       IZ=IWZ+I
  251.       IF (I .EQ. 1) GOTO 190
  252.       DO 180 J=2,I
  253.       W(IZ)=ZERO
  254.   180 IZ=IZ+N
  255.   190 IR=IWR+(I+I*I)/2
  256.       W(IZ)=ONE/W(IR)
  257.       IF (I .EQ. N) GOTO 220
  258.       IZA=IZ
  259.       DO 210 J=I,NM
  260.       IR=IR+I
  261.       SUM=ZERO
  262.       DO 200 K=IZA,IZ,N
  263.       SUM=SUM+W(K)*W(IR)
  264.   200 IR=IR+1
  265.       IZ=IZ+N
  266.   210 W(IZ)=-SUM/W(IR)
  267.   220 CONTINUE
  268. C
  269. C     SET THE INITIAL VALUES OF SOME VARIABLES.
  270. C     ITERC COUNTS THE NUMBER OF ITERATIONS.
  271. C     ITREF IS SET TO ONE WHEN ITERATIVE REFINEMENT IS REQUIRED.
  272. C     JFINC INDICATES WHEN TO TEST FOR AN INCREASE IN F.
  273. C
  274.       ITERC=1
  275.       ITREF=0
  276.       JFINC=-KFINC
  277. C
  278. C     SET X TO ZERO AND SET THE CORRESPONDING RESIDUALS OF THE
  279. C     KUHN-TUCKER CONDITIONS.
  280. C
  281.   230 IFLAG=1
  282.       IWS=IWW-N
  283.       DO 240 I=1,N
  284.       X(I)=ZERO
  285.       IW=IWW+I
  286.       W(IW)=GRAD(I)
  287.       IF (I .GT. NACT) GOTO 240
  288.       W(I)=ZERO
  289.       IS=IWS+I
  290.       K=IACT(I)
  291.       IF (K .LE. M) GOTO 235
  292.       IF (K .GT. MN) GOTO 234
  293.       K1=K-M
  294.       W(IS)=XL(K1)
  295.       GOTO 240
  296.   234 K1=K-MN
  297.       W(IS)=-XU(K1)
  298.       GOTO 240
  299.   235 W(IS)=B(K)
  300.   240 CONTINUE
  301.       XMAG=ZERO
  302.       VFACT=1.D+0
  303.       IF (NACT) 340,340,280
  304. C
  305. C     SET THE RESIDUALS OF THE KUHN-TUCKER CONDITIONS FOR GENERAL X.
  306. C
  307.   250 IFLAG=2
  308.       IWS=IWW-N
  309.       DO 260 I=1,N
  310.       IW=IWW+I
  311.       W(IW)=GRAD(I)
  312.       IF (LQL) GOTO 259
  313.       ID=IWD+I
  314.       W(ID)=ZERO
  315.       DO 251 J=I,N
  316.   251 W(ID)=W(ID)+G(I,J)*X(J)
  317.       DO 252 J=1,I
  318.       ID=IWD+J
  319.   252 W(IW)=W(IW)+G(J,I)*W(ID)
  320.       GOTO 260
  321.   259 DO 261 J=1,N
  322.   261 W(IW)=W(IW)+G(I,J)*X(J)
  323.   260 CONTINUE
  324.       IF (NACT .EQ. 0) GOTO 340
  325.       DO 270 K=1,NACT
  326.       KK=IACT(K)
  327.       IS=IWS+K
  328.       IF (KK .GT. M) GOTO 265
  329.       W(IS)=B(KK)
  330.       DO 264 I=1,N
  331.       IW=IWW+I
  332.       W(IW)=W(IW)-W(K)*A(KK,I)
  333.   264 W(IS)=W(IS)-X(I)*A(KK,I)
  334.       GOTO 270
  335.   265 IF (KK .GT. MN) GOTO 266
  336.       K1=KK-M
  337.       IW=IWW+K1
  338.       W(IW)=W(IW)-W(K)
  339.       W(IS)=XL(K1)-X(K1)
  340.       GOTO 270
  341.   266 K1=KK-MN
  342.       IW=IWW+K1
  343.       W(IW)=W(IW)+W(K)
  344.       W(IS)=-XU(K1)+X(K1)
  345.   270 CONTINUE
  346. C
  347. C     PRE-MULTIPLY THE VECTOR IN THE S-PARTITION OF W BY THE
  348. C     INVERS OF R TRANSPOSE.
  349. C
  350.   280 IR=IWR
  351.       IP=IWW+1
  352.       IPP=IWW+N
  353.       IL=IWS+1
  354.       IU=IWS+NACT
  355.       DO 310 I=IL,IU
  356.       SUM=ZERO
  357.       IF (I .EQ. IL) GOTO 300
  358.       JU=I-1
  359.       DO 290 J=IL,JU
  360.       IR=IR+1
  361.   290 SUM=SUM+W(IR)*W(J)
  362.   300 IR=IR+1
  363.   310 W(I)=(W(I)-SUM)/W(IR)
  364. C
  365. C     SHIFT X TO SATISFY THE ACTIVE CONSTRAINTS AND MAKE THE
  366. C     CORRESPONDING CHANGE TO THE GRADIENT RESIDUALS.
  367. C
  368.       DO 330 I=1,N
  369.       IZ=IWZ+I
  370.       SUM=ZERO
  371.       DO 320 J=IL,IU
  372.       SUM=SUM+W(J)*W(IZ)
  373.   320 IZ=IZ+N
  374.       X(I)=X(I)+SUM
  375.       IF (LQL) GOTO 329
  376.       ID=IWD+I
  377.       W(ID)=ZERO
  378.       DO 321 J=I,N
  379.   321 W(ID)=W(ID)+G(I,J)*SUM
  380.       IW=IWW+I
  381.       DO 322 J=1,I
  382.       ID=IWD+J
  383.   322 W(IW)=W(IW)+G(J,I)*W(ID)
  384.       GOTO 330
  385.   329 DO 331 J=1,N
  386.       IW=IWW+J
  387.   331 W(IW)=W(IW)+SUM*G(I,J)
  388.   330 CONTINUE
  389. C
  390. C     FORM THE SCALAR PRODUCT OF THE CURRENT GRADIENT RESIDUALS
  391. C     WITH EACH COLUMN OF Z.
  392. C
  393.   340 KFLAG=1
  394.       GOTO 930
  395.   350 IF (NACT .EQ. N) GOTO 380
  396. C
  397. C     SHIFT X SO THAT IT SATISFIES THE REMAINING KUHN-TUCKER
  398. C     CONDITIONS.
  399. C
  400.       IL=IWS+NACT+1
  401.       IZA=IWZ+NACT*N
  402.       DO 370 I=1,N
  403.       SUM=ZERO
  404.       IZ=IZA+I
  405.       DO 360 J=IL,IWW
  406.       SUM=SUM+W(IZ)*W(J)
  407.   360 IZ=IZ+N
  408.   370 X(I)=X(I)-SUM
  409.       INFO=0
  410.       IF (NACT .EQ. 0) GOTO 410
  411. C
  412. C     UPDATE THE LAGRANGE MULTIPLIERS.
  413. C
  414.   380 LFLAG=3
  415.       GOTO 740
  416.   390 DO 400 K=1,NACT
  417.       IW=IWW+K
  418.   400 W(K)=W(K)+W(IW)
  419. C
  420. C     REVISE THE VALUES OF XMAG.
  421. C     BRANCH IF ITERATIVE REFINEMENT IS REQUIRED.
  422. C
  423.   410 JFLAG=1
  424.       GOTO 910
  425.   420 IF (IFLAG .EQ. ITREF) GOTO 250
  426. C
  427. C     DELETE A CONSTRAINT IF A LAGRANGE MULTIPLIER OF AN
  428. C     INEQUALITY CONSTRAINT IS NEGATIVE.
  429. C
  430.       KDROP=0
  431.       GOTO 440
  432.   430 KDROP=KDROP+1
  433.       IF (W(KDROP) .GE. ZERO) GOTO 440
  434.       IF (IACT(KDROP) .LE. MEQ) GOTO 440
  435.       NU=NACT
  436.       MFLAG=1
  437.       GOTO 800
  438.   440 IF (KDROP .LT. NACT) GOTO 430
  439. C
  440. C     SEEK THE GREATEAST NORMALISED CONSTRAINT VIOLATION, DISREGARDING
  441. C     ANY THAT MAY BE DUE TO COMPUTER ROUNDING ERRORS.
  442. C
  443.   450 CVMAX=ZERO
  444.       IF (M .LE. 0) GOTO 481
  445.       DO 480 K=1,M
  446.       IA=IWA+K
  447.       IF (W(IA) .LE. ZERO) GOTO 480
  448.       SUM=-B(K)
  449.       DO 460 I=1,N
  450.   460 SUM=SUM+X(I)*A(K,I)
  451.       SUMX=-SUM*W(IA)
  452.       IF (K .LE. MEQ) SUMX=DABS(SUMX)
  453.       IF (SUMX .LE. CVMAX) GOTO 480
  454.       TEMP=DABS(B(K))
  455.       DO 470 I=1,N
  456.   470 TEMP=TEMP+DABS(X(I)*A(K,I))
  457.       TEMPA=TEMP+DABS(SUM)
  458.       IF (TEMPA .LE. TEMP) GOTO 480
  459.       TEMP=TEMP+ONHA*DABS(SUM)
  460.       IF (TEMP .LE. TEMPA) GOTO 480
  461.       CVMAX=SUMX
  462.       RES=SUM
  463.       KNEXT=K
  464.   480 CONTINUE
  465.   481 DO 485 K=1,N
  466.       LOWER=.TRUE.
  467.       IA=IWA+M+K
  468.       IF (W(IA) .LE. ZERO) GOTO 485
  469.       SUM=XL(K)-X(K)
  470.       IF (SUM) 482,485,483
  471.   482 SUM=X(K)-XU(K)
  472.       LOWER=.FALSE.
  473.   483 IF (SUM .LE. CVMAX) GOTO 485
  474.       CVMAX=SUM
  475.       RES=-SUM
  476.       KNEXT=K+M
  477.       IF (LOWER) GOTO 485
  478.       KNEXT=K+MN
  479.   485 CONTINUE
  480. C
  481. C     TEST FOR CONVERGENCE
  482. C
  483.       INFO=0
  484.       IF (CVMAX .LE. VSMALL) GOTO 700
  485. C
  486. C     RETURN IF, DUE TO ROUNDING ERRORS, THE ACTUAL CHANGE IN
  487. C     X MAY NOT INCREASE THE OBJECTIVE FUNCTION
  488. C
  489.       JFINC=JFINC+1
  490.       IF (JFINC .EQ. 0) GOTO 510
  491.       IF (JFINC .NE. IFINC) GOTO 530
  492.       FDIFF=ZERO
  493.       FDIFFA=ZERO
  494.       DO 500 I=1,N
  495.       SUM=TWO*GRAD(I)
  496.       SUMX=DABS(SUM)
  497.       IF (LQL) GOTO 489
  498.       ID=IWD+I
  499.       W(ID)=ZERO
  500.       DO 486 J=I,N
  501.       IX=IWX+J
  502.   486 W(ID)=W(ID)+G(I,J)*(W(IX)+X(J))
  503.       DO 487 J=1,I
  504.       ID=IWD+J
  505.       TEMP=G(J,I)*W(ID)
  506.       SUM=SUM+TEMP
  507.   487 SUMX=SUMX+DABS(TEMP)
  508.       GOTO 495
  509.   489 DO 490 J=1,N
  510.       IX=IWX+J
  511.       TEMP=G(I,J)*(W(IX)+X(J))
  512.       SUM=SUM+TEMP
  513.   490 SUMX=SUMX+DABS(TEMP)
  514.   495 IX=IWX+I
  515.       FDIFF=FDIFF+SUM*(X(I)-W(IX))
  516.   500 FDIFFA=FDIFFA+SUMX*DABS(X(I)-W(IX))
  517.       INFO=2
  518.       SUM=FDIFFA+FDIFF
  519.       IF (SUM .LE. FDIFFA) GOTO 700
  520.       TEMP=FDIFFA+ONHA*FDIFF
  521.       IF (TEMP .LE. SUM) GOTO 700
  522.       JFINC=0
  523.       INFO=0
  524.   510 DO 520 I=1,N
  525.       IX=IWX+I
  526.   520 W(IX)=X(I)
  527. C
  528. C     FORM THE SCALAR PRODUCT OF THE NEW CONSTRAINT NORMAL WITH EACH
  529. C     COLUMN OF Z. PARNEW WILL BECOME THE LAGRANGE MULTIPLIER OF
  530. C     THE NEW CONSTRAINT.
  531. C
  532.   530 ITERC=ITERC+1
  533.       IF (ITERC.LE.MAXIT) GOTO 531
  534.       INFO=1
  535.       GOTO 710
  536.   531 CONTINUE
  537.       IWS=IWR+(NACT+NACT*NACT)/2
  538.       IF (KNEXT .GT. M) GOTO 541
  539.       DO 540 I=1,N
  540.       IW=IWW+I
  541.   540 W(IW)=A(KNEXT,I)
  542.       GOTO 549
  543.   541 DO 542 I=1,N
  544.       IW=IWW+I
  545.   542 W(IW)=ZERO
  546.       K1=KNEXT-M
  547.       IF (K1 .GT. N) GOTO 545
  548.       IW=IWW+K1
  549.       W(IW)=ONE
  550.       IZ=IWZ+K1
  551.       DO 543 I=1,N
  552.       IS=IWS+I
  553.       W(IS)=W(IZ)
  554.   543 IZ=IZ+N
  555.       GOTO 550
  556.   545 K1=KNEXT-MN
  557.       IW=IWW+K1
  558.       W(IW)=-ONE
  559.       IZ=IWZ+K1
  560.       DO 546 I=1,N
  561.       IS=IWS+I
  562.       W(IS)=-W(IZ)
  563.   546 IZ=IZ+N
  564.       GOTO 550
  565.   549 KFLAG=2
  566.       GOTO 930
  567.   550 PARNEW=ZERO
  568. C
  569. C     APPLY GIVENS ROTATIONS TO MAKE THE LAST (N-NACT-2) SCALAR
  570. C     PRODUCTS EQUAL TO ZERO.
  571. C
  572.       IF (NACT .EQ. N) GOTO 570
  573.       NU=N
  574.       NFLAG=1
  575.       GOTO 860
  576. C
  577. C     BRANCH IF THERE IS NO NEED TO DELETE A CONSTRAINT.
  578. C
  579.   560 IS=IWS+NACT
  580.       IF (NACT .EQ. 0) GOTO 640
  581.       SUMA=ZERO
  582.       SUMB=ZERO
  583.       SUMC=ZERO
  584.       IZ=IWZ+NACT*N
  585.       DO 563 I=1,N
  586.       IZ=IZ+1
  587.       IW=IWW+I
  588.       SUMA=SUMA+W(IW)*W(IZ)
  589.       SUMB=SUMB+DABS(W(IW)*W(IZ))
  590.   563 SUMC=SUMC+W(IZ)**2
  591.       TEMP=SUMB+.1D+0*DABS(SUMA)
  592.       TEMPA=SUMB+.2D+0*DABS(SUMA)
  593.       IF (TEMP .LE. SUMB) GOTO 570
  594.       IF (TEMPA .LE. TEMP) GOTO 570
  595.       IF (SUMB .GT. VSMALL) GOTO 5
  596.       GOTO 570
  597.     5 SUMC=DSQRT(SUMC)
  598.       IA=IWA+KNEXT
  599.       IF (KNEXT .LE. M) SUMC=SUMC/W(IA)
  600.       TEMP=SUMC+.1D+0*DABS(SUMA)
  601.       TEMPA=SUMC+.2D+0*DABS(SUMA)
  602.       IF (TEMP .LE. SUMC) GOTO 567
  603.       IF (TEMPA .LE. TEMP) GOTO 567
  604.       GOTO 640
  605. C
  606. C     CALCULATE THE MULTIPLIERS FOR THE NEW CONSTRAINT NORMAL
  607. C     EXPRESSED IN TERMS OF THE ACTIVE CONSTRAINT NORMALS.
  608. C     THEN WORK OUT WHICH CONTRAINT TO DROP.
  609. C
  610.   567 LFLAG=4
  611.       GOTO 740
  612.   570 LFLAG=1
  613.       GOTO 740
  614. C
  615. C     COMPLETE THE TEST FOR LINEARLY DEPENDENT CONSTRAINTS.
  616. C
  617.   571 IF (KNEXT .GT. M) GOTO 574
  618.       DO 573 I=1,N
  619.       SUMA=A(KNEXT,I)
  620.       SUMB=DABS(SUMA)
  621.       IF (NACT.EQ.0) GOTO 581
  622.       DO 572 K=1,NACT
  623.       KK=IACT(K)
  624.       IF (KK.LE.M) GOTO 568
  625.       KK=KK-M
  626.       TEMP=ZERO
  627.       IF (KK.EQ.I) TEMP=W(IWW+KK)
  628.       KK=KK-N
  629.       IF (KK.EQ.I) TEMP=-W(IWW+KK)
  630.       GOTO 569
  631.   568 CONTINUE
  632.       IW=IWW+K
  633.       TEMP=W(IW)*A(KK,I)
  634.   569 CONTINUE
  635.       SUMA=SUMA-TEMP
  636.   572 SUMB=SUMB+DABS(TEMP)
  637.   581 IF (SUMA .LE. VSMALL) GOTO 573
  638.       TEMP=SUMB+.1D+0*DABS(SUMA)
  639.       TEMPA=SUMB+.2D+0*DABS(SUMA)
  640.       IF (TEMP .LE. SUMB) GOTO 573
  641.       IF (TEMPA .LE. TEMP) GOTO 573
  642.       GOTO 630
  643.   573 CONTINUE
  644.       LFLAG=1
  645.       GOTO 775
  646.   574 K1=KNEXT-M
  647.       IF (K1 .GT. N) K1=K1-N
  648.       DO 578 I=1,N
  649.       SUMA=ZERO
  650.       IF (I .NE. K1) GOTO 575
  651.       SUMA=ONE
  652.       IF (KNEXT .GT. MN) SUMA=-ONE
  653.   575 SUMB=DABS(SUMA)
  654.       IF (NACT.EQ.0) GOTO 582
  655.       DO 577 K=1,NACT
  656.       KK=IACT(K)
  657.       IF (KK .LE. M) GOTO 579
  658.       KK=KK-M
  659.       TEMP=ZERO
  660.       IF (KK.EQ.I) TEMP=W(IWW+KK)
  661.       KK=KK-N
  662.       IF (KK.EQ.I) TEMP=-W(IWW+KK)
  663.       GOTO 576
  664.   579 IW=IWW+K
  665.       TEMP=W(IW)*A(KK,I)
  666.   576 SUMA=SUMA-TEMP
  667.   577 SUMB=SUMB+DABS(TEMP)
  668.   582 TEMP=SUMB+.1D+0*DABS(SUMA)
  669.       TEMPA=SUMB+.2D+0*DABS(SUMA)
  670.       IF (TEMP .LE. SUMB) GOTO 578
  671.       IF (TEMPA .LE. TEMP) GOTO 578
  672.       GOTO 630
  673.   578 CONTINUE
  674.       LFLAG=1
  675.       GOTO 775
  676. C
  677. C     BRANCH IF THE CONTRAINTS ARE INCONSISTENT.
  678. C
  679.   580 INFO=-KNEXT
  680.       IF (KDROP .EQ. 0) GOTO 700
  681.       PARINC=RATIO
  682.       PARNEW=PARINC
  683. C
  684. C     REVISE THE LAGRANGE MULTIPLIERS OF THE ACTIVE CONSTRAINTS.
  685. C
  686.   590 IF (NACT.EQ.0) GOTO 601
  687.       DO 600 K=1,NACT
  688.       IW=IWW+K
  689.       W(K)=W(K)-PARINC*W(IW)
  690.       IF (IACT(K) .GT. MEQ) W(K)=DMAX1(ZERO,W(K))
  691.   600 CONTINUE
  692.   601 IF (KDROP .EQ. 0) GOTO 680
  693. C
  694. C     DELETE THE CONSTRAINT TO BE DROPPED.
  695. C     SHIFT THE VECTOR OF SCALAR PRODUCTS.
  696. C     THEN, IF APPROPRIATE, MAKE ONE MORE SCALAR PRODUCT ZERO.
  697. C
  698.       NU=NACT+1
  699.       MFLAG=2
  700.       GOTO 800
  701.   610 IWS=IWS-NACT-1
  702.       NU=MIN0(N,NU)
  703.       DO 620 I=1,NU
  704.       IS=IWS+I
  705.       J=IS+NACT
  706.   620 W(IS)=W(J+1)
  707.       NFLAG=2
  708.       GOTO 860
  709. C
  710. C     CALCULATE THE STEP TO THE VIOLATED CONSTRAINT.
  711. C
  712.   630 IS=IWS+NACT
  713.   640 SUMY=W(IS+1)
  714.       STEP=-RES/SUMY
  715.       PARINC=STEP/SUMY
  716.       IF (NACT .EQ. 0) GOTO 660
  717. C
  718. C     CALCULATE THE CHANGES TO THE LAGRANGE MULTIPLIERS, AND REDUCE
  719. C     THE STEP ALONG THE NEW SEARCH DIRECTION IF NECESSARY.
  720. C
  721.       LFLAG=2
  722.       GOTO 740
  723.   650 IF (KDROP .EQ. 0) GOTO 660
  724.       TEMP=ONE-RATIO/PARINC
  725.       IF (TEMP .LE. ZERO) KDROP=0
  726.       IF (KDROP .EQ. 0) GOTO 660
  727.       STEP=RATIO*SUMY
  728.       PARINC=RATIO
  729.       RES=TEMP*RES
  730. C
  731. C     UPDATE X AND THE LAGRANGE MULTIPIERS.
  732. C     DROP A CONSTRAINT IF THE FULL STEP IS NOT TAKEN.
  733. C
  734.   660 IWY=IWZ+NACT*N
  735.       DO 670 I=1,N
  736.       IY=IWY+I
  737.   670 X(I)=X(I)+STEP*W(IY)
  738.       PARNEW=PARNEW+PARINC
  739.       IF (NACT .GE. 1) GOTO 590
  740. C
  741. C     ADD THE NEW CONSTRAINT TO THE ACTIVE SET.
  742. C
  743.   680 NACT=NACT+1
  744.       W(NACT)=PARNEW
  745.       IACT(NACT)=KNEXT
  746.       IA=IWA+KNEXT
  747.       IF (KNEXT .GT. MN) IA=IA-N
  748.       W(IA)=-W(IA)
  749. C
  750. C     ESTIMATE THE MAGNITUDE OF X. THEN BEGIN A NEW ITERATION,
  751. C     RE-INITILISING X IF THIS MAGNITUDE IS SMALL.
  752. C
  753.       JFLAG=2
  754.       GOTO 910
  755.   690 IF (SUM .LT. (XMAGR*XMAG)) GOTO 230
  756.       IF (ITREF) 450,450,250
  757. C
  758. C     INITIATE ITERATIVE REFINEMENT IF IT HAS NOT YET BEEN USED,
  759. C     OR RETURN AFTER RESTORING THE DIAGONAL ELEMENTS OF G.
  760. C
  761.   700 IF (ITERC .EQ. 0) GOTO 710
  762.       ITREF=ITREF+1
  763.       JFINC=-1
  764.       IF (ITREF .EQ. 1) GOTO 250
  765.   710 IF (.NOT. LQL) RETURN
  766.       DO 720 I=1,N
  767.       ID=IWD+I
  768.   720 G(I,I)=W(ID)
  769.   730 RETURN
  770. C
  771. C
  772. C     THE REMAINIG INSTRUCTIONS ARE USED AS SUBROUTINES.
  773. C
  774. C
  775. C********************************************************************
  776. C
  777. C
  778. C     CALCULATE THE LAGRANGE MULTIPLIERS BY PRE-MULTIPLYING THE
  779. C     VECTOR IN THE S-PARTITION OF W BY THE INVERSE OF R.
  780. C
  781.   740 IR=IWR+(NACT+NACT*NACT)/2
  782.       I=NACT
  783.       SUM=ZERO
  784.       GOTO 770
  785.   750 IRA=IR-1
  786.       SUM=ZERO
  787.       IF (NACT.EQ.0) GOTO 761
  788.       DO 760 J=I,NACT
  789.       IW=IWW+J
  790.       SUM=SUM+W(IRA)*W(IW)
  791.   760 IRA=IRA+J
  792.   761 IR=IR-I
  793.       I=I-1
  794.   770 IW=IWW+I
  795.       IS=IWS+I
  796.       W(IW)=(W(IS)-SUM)/W(IR)
  797.       IF (I .GT. 1) GOTO 750
  798.       IF (LFLAG .EQ. 3) GOTO 390
  799.       IF (LFLAG .EQ. 4) GOTO 571
  800. C
  801. C     CALCULATE THE NEXT CONSTRAINT TO DROP.
  802. C
  803.   775 IP=IWW+1
  804.       IPP=IWW+NACT
  805.       KDROP=0
  806.       IF (NACT.EQ.0) GOTO 791
  807.       DO 790 K=1,NACT
  808.       IF (IACT(K) .LE. MEQ) GOTO 790
  809.       IW=IWW+K
  810.       IF ((RES*W(IW)) .GE. ZERO) GOTO 790
  811.       TEMP=W(K)/W(IW)
  812.       IF (KDROP .EQ. 0) GOTO 780
  813.       IF (DABS(TEMP) .GE. DABS(RATIO)) GOTO 790
  814.   780 KDROP=K
  815.       RATIO=TEMP
  816.   790 CONTINUE
  817.   791 GOTO (580,650), LFLAG
  818. C
  819. C
  820. C********************************************************************
  821. C
  822. C
  823. C     DROP THE CONSTRAINT IN POSITION KDROP IN THE ACTIVE SET.
  824. C
  825.   800 IA=IWA+IACT(KDROP)
  826.       IF (IACT(KDROP) .GT. MN) IA=IA-N
  827.       W(IA)=-W(IA)
  828.       IF (KDROP .EQ. NACT) GOTO 850
  829. C
  830. C     SET SOME INDICES AND CALCULATE THE ELEMENTS OF THE NEXT
  831. C     GIVENS ROTATION.
  832. C
  833.       IZ=IWZ+KDROP*N
  834.       IR=IWR+(KDROP+KDROP*KDROP)/2
  835.   810 IRA=IR
  836.       IR=IR+KDROP+1
  837.       TEMP=DMAX1(DABS(W(IR-1)),DABS(W(IR)))
  838.       SUM=TEMP*DSQRT((W(IR-1)/TEMP)**2+(W(IR)/TEMP)**2)
  839.       GA=W(IR-1)/SUM
  840.       GB=W(IR)/SUM
  841. C
  842. C     EXCHANGE THE COLUMNS OF R.
  843. C
  844.       DO 820 I=1,KDROP
  845.       IRA=IRA+1
  846.       J=IRA-KDROP
  847.       TEMP=W(IRA)
  848.       W(IRA)=W(J)
  849.   820 W(J)=TEMP
  850.       W(IR)=ZERO
  851. C
  852. C     APPLY THE ROTATION TO THE ROWS OF R.
  853. C
  854.       W(J)=SUM
  855.       KDROP=KDROP+1
  856.       DO 830 I=KDROP,NU
  857.       TEMP=GA*W(IRA)+GB*W(IRA+1)
  858.       W(IRA+1)=GA*W(IRA+1)-GB*W(IRA)
  859.       W(IRA)=TEMP
  860.   830 IRA=IRA+I
  861. C
  862. C     APPLY THE ROTATION TO THE COLUMNS OF Z.
  863. C
  864.       DO 840 I=1,N
  865.       IZ=IZ+1
  866.       J=IZ-N
  867.       TEMP=GA*W(J)+GB*W(IZ)
  868.       W(IZ)=GA*W(IZ)-GB*W(J)
  869.   840 W(J)=TEMP
  870. C
  871. C     REVISE IACT AND THE LAGRANGE MULTIPLIERS.
  872. C
  873.       IACT(KDROP-1)=IACT(KDROP)
  874.       W(KDROP-1)=W(KDROP)
  875.       IF (KDROP .LT. NACT) GOTO 810
  876.   850 NACT=NACT-1
  877.       GOTO (250,610), MFLAG
  878. C
  879. C
  880. C********************************************************************
  881. C
  882. C
  883. C     APPLY GIVENS ROTATION TO REDUCE SOME OF THE SCALAR
  884. C     PRODUCTS IN THE S-PARTITION OF W TO ZERO.
  885. C
  886.   860 IZ=IWZ+NU*N
  887.   870 IZ=IZ-N
  888.   880 IS=IWS+NU
  889.       NU=NU-1
  890.       IF (NU .EQ. NACT) GOTO 900
  891.       IF (W(IS) .EQ. ZERO) GOTO 870
  892.       TEMP=DMAX1(DABS(W(IS-1)),DABS(W(IS)))
  893.       SUM=TEMP*DSQRT((W(IS-1)/TEMP)**2+(W(IS)/TEMP)**2)
  894.       GA=W(IS-1)/SUM
  895.       GB=W(IS)/SUM
  896.       W(IS-1)=SUM
  897.       DO 890 I=1,N
  898.       K=IZ+N
  899.       TEMP=GA*W(IZ)+GB*W(K)
  900.       W(K)=GA*W(K)-GB*W(IZ)
  901.       W(IZ)=TEMP
  902.   890 IZ=IZ-1
  903.       GOTO 880
  904.   900 GOTO (560,630), NFLAG
  905. C
  906. C
  907. C********************************************************************
  908. C
  909. C
  910. C     CALCULATE THE MAGNITUDE OF X AN REVISE XMAG.
  911. C
  912.   910 SUM=ZERO
  913.       DO 920 I=1,N
  914.       SUM=SUM+DABS(X(I))*VFACT*(DABS(GRAD(I))+DABS(G(I,I)*X(I)))
  915.       IF (LQL) GOTO 920
  916.       IF (SUM .LT. 1.D-30) GOTO 920
  917.       VFACT=1.D-10*VFACT
  918.       SUM=1.D-10*SUM
  919.       XMAG=1.D-10*XMAG
  920.   920 CONTINUE
  921.   925 XMAG=DMAX1(XMAG,SUM)
  922.       GOTO (420,690), JFLAG
  923. C
  924. C
  925. C********************************************************************
  926. C
  927. C
  928. C     PRE-MULTIPLY THE VECTOR IN THE W-PARTITION OF W BY Z TRANSPOSE.
  929. C
  930.   930 JL=IWW+1
  931.       IZ=IWZ
  932.       DO 940 I=1,N
  933.       IS=IWS+I
  934.       W(IS)=ZERO
  935.       IWWN=IWW+N
  936.       DO 940 J=JL,IWWN
  937.       IZ=IZ+1
  938.   940 W(IS)=W(IS)+W(IZ)*W(J)
  939.       GOTO (350,550), KFLAG
  940.       RETURN
  941.       END
  942.