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 / ql0001.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  227 lines
  1.       SUBROUTINE QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,
  2.      1           X,U,IOUT,IFAIL,IPRINT,WAR,LWAR,IWAR,LIWAR)
  3. c
  4. cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  5. c
  6. c                     !!!! NOTICE !!!!
  7. c
  8. c 1. The routines contained in this file are due to Prof. K.Schittkowski
  9. c    of the University of Bayreuth, Germany (modification of routines
  10. c    due to Prof. MJD Powell at the University of Cambridge).  They can
  11. c    be freely distributed.
  12. c
  13. c 2. A minor modification was performed at the University of Maryland. 
  14. c    It is marked in the code by "c umd".
  15. c
  16. c                                      A.L. Tits and J.L. Zhou
  17. c                                      University of Maryland
  18. C
  19. C***********************************************************************
  20. C
  21. C
  22. C             SOLUTION OF QUADRATIC PROGRAMMING PROBLEMS
  23. C
  24. C
  25. C
  26. C   QL0001 SOLVES THE QUADRATIC PROGRAMMING PROBLEM
  27. C
  28. C   MINIMIZE        .5*X'*C*X + D'*X
  29. C   SUBJECT TO      A(J)*X  +  B(J)   =  0  ,  J=1,...,ME
  30. C                   A(J)*X  +  B(J)  >=  0  ,  J=ME+1,...,M
  31. C                   XL  <=  X  <=  XU
  32. C   
  33. C   HERE C MUST BE AN N BY N SYMMETRIC AND POSITIVE MATRIX, D AN N-DIMENSIONAL
  34. C   VECTOR, A AN M BY N MATRIX AND B AN M-DIMENSIONAL VECTOR. THE ABOVE 
  35. C   SITUATION IS INDICATED BY IWAR(1)=1. ALTERNATIVELY, I.E. IF IWAR(1)=0,
  36. C   THE OBJECTIVE FUNCTION MATRIX CAN ALSO BE PROVIDED IN FACTORIZED FORM.
  37. C   IN THIS CASE, C IS AN UPPER TRIANGULAR MATRIX.
  38. C
  39. C   THE SUBROUTINE REORGANIZES SOME DATA SO THAT THE PROBLEM CAN BE SOLVED
  40. C   BY A MODIFICATION OF AN ALGORITHM PROPOSED BY POWELL (1983).
  41. C
  42. C
  43. C   USAGE:
  44. C
  45. C      QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,X,U,IOUT,IFAIL,IPRINT,
  46. C             WAR,LWAR,IWAR,LIWAR)
  47. C
  48. C
  49. C   DEFINITION OF THE PARAMETERS:
  50. C
  51. C   M :        TOTAL NUMBER OF CONSTRAINTS.
  52. C   ME :       NUMBER OF EQUALITY CONSTRAINTS.
  53. C   MMAX :     ROW DIMENSION OF A. MMAX MUST BE AT LEAST ONE AND GREATER
  54. C              THAN M.
  55. C   N :        NUMBER OF VARIABLES.
  56. C   NMAX :     ROW DIMENSION OF C. NMAX MUST BE GREATER OR EQUAL TO N.
  57. C   MNN :      MUST BE EQUAL TO M + N + N.
  58. C   C(NMAX,NMAX): OBJECTIVE FUNCTION MATRIX WHICH SHOULD BE SYMMETRIC AND
  59. C              POSITIVE DEFINITE. IF IWAR(1) = 0, C IS SUPPOSED TO BE THE
  60. C              CHOLESKEY-FACTOR OF ANOTHER MATRIX, I.E. C IS UPPER
  61. C              TRIANGULAR.
  62. C   D(NMAX) :  CONTAINS THE CONSTANT VECTOR OF THE OBJECTIVE FUNCTION.
  63. C   A(MMAX,NMAX): CONTAINS THE DATA MATRIX OF THE LINEAR CONSTRAINTS.
  64. C   B(MMAX) :  CONTAINS THE CONSTANT DATA OF THE LINEAR CONSTRAINTS.
  65. C   XL(N),XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR THE VARIABLES.
  66. C   X(N) :     ON RETURN, X CONTAINS THE OPTIMAL SOLUTION VECTOR.
  67. C   U(MNN) :   ON RETURN, U CONTAINS THE LAGRANGE MULTIPLIERS. THE FIRST
  68. C              M POSITIONS ARE RESERVED FOR THE MULTIPLIERS OF THE M
  69. C              LINEAR CONSTRAINTS AND THE SUBSEQUENT ONES FOR THE 
  70. C              MULTIPLIERS OF THE LOWER AND UPPER BOUNDS. ON SUCCESSFUL
  71. C              TERMINATION, ALL VALUES OF U WITH RESPECT TO INEQUALITIES 
  72. C              AND BOUNDS SHOULD BE GREATER OR EQUAL TO ZERO.
  73. C   IOUT :     INTEGER INDICATING THE DESIRED OUTPUT UNIT NUMBER, I.E.
  74. C              ALL WRITE-STATEMENTS START WITH 'WRITE(IOUT,... '.
  75. C   IFAIL :    SHOWS THE TERMINATION REASON.
  76. C      IFAIL = 0 :   SUCCESSFUL RETURN.
  77. C      IFAIL = 1 :   TOO MANY ITERATIONS (MORE THAN 40*(N+M)).
  78. C      IFAIL = 2 :   ACCURACY INSUFFICIENT TO SATISFY CONVERGENCE
  79. C                    CRITERION.
  80. C      IFAIL = 5 :   LENGTH OF A WORKING ARRAY IS TOO SHORT.
  81. C      IFAIL > 10 :  THE CONSTRAINTS ARE INCONSISTENT.
  82. C   IPRINT :   OUTPUT CONTROL.
  83. C      IPRINT = 0 :  NO OUTPUT OF QL0001.
  84. C      IPRINT > 0 :  BRIEF OUTPUT IN ERROR CASES.
  85. C   WAR(LWAR) : REAL WORKING ARRAY. THE LENGTH LWAR SHOULD BE GRATER THAN
  86. C               NMAX*(3*NMAX+15)/2 + 2*M.
  87. C   IWAR(LIWAR): INTEGER WORKING ARRAY. THE LENGTH LIWAR SHOULD BE AT
  88. C              LEAST N.
  89. C              IF IWAR(1)=0 INITIALLY, THEN THE CHOLESKY DECOMPOSITION
  90. C              WHICH IS REQUIRED BY THE DUAL ALGORITHM TO GET THE FIRST
  91. C              UNCONSTRAINED MINIMUM OF THE OBJECTIVE FUNCTION, IS
  92. C              PERFORMED INTERNALLY. OTHERWISE, I.E. IF IWAR(1)=1, THEN
  93. C              IT IS ASSUMED THAT THE USER PROVIDES THE INITIAL FAC-
  94. C              TORIZATION BY HIMSELF AND STORES IT IN THE UPPER TRIAN-
  95. C              GULAR PART OF THE ARRAY C.
  96. C
  97. C   A NAMED COMMON-BLOCK  /CMACHE/EPS   MUST BE PROVIDED BY THE USER,
  98. C   WHERE EPS DEFINES A GUESS FOR THE UNDERLYING MACHINE PRECISION.
  99. C
  100. C
  101. C   AUTHOR:    K. SCHITTKOWSKI,
  102. C              MATHEMATISCHES INSTITUT,
  103. C              UNIVERSITAET BAYREUTH,
  104. C              8580 BAYREUTH,
  105. C              GERMANY, F.R.
  106. C
  107. C
  108. C   VERSION:   1.4  (MARCH, 1987)  
  109. C
  110. C
  111. C*********************************************************************
  112. C
  113. C
  114.       INTEGER NMAX,MMAX,N,MNN,LWAR,LIWAR
  115.       DIMENSION C(NMAX,NMAX),D(NMAX),A(MMAX,NMAX),B(MMAX),
  116.      1      XL(N),XU(N),X(N),U(MNN),WAR(LWAR),IWAR(LIWAR)
  117.       DOUBLE PRECISION C,D,A,B,X,XL,XU,U,WAR,DIAG,ZERO,
  118.      1      EPS,QPEPS,TEN
  119.       INTEGER M,ME,IOUT,IFAIL,IPRINT,IWAR,INW1,INW2,IN,J,LW,MN,I,         
  120.      1      IDIAG,INFO,NACT,MAXIT
  121.       LOGICAL LQL
  122. C
  123. C     INTRINSIC FUNCTIONS:  DSQRT
  124. C
  125.       COMMON /CMACHE/EPS
  126. C
  127. C     CONSTANT DATA
  128. C
  129. c#################################################################
  130. c
  131.  
  132.       if(c(nmax,nmax).eq.0.d0) c(nmax,nmax)=eps
  133. c
  134. c umd
  135. c  This prevents a subsequent more major modification of the Hessian
  136. c  matrix in the important case when a minmax problem (yielding a 
  137. c  singular Hessian matrix) is being solved.
  138. c                                 ----UMCP, April 1991, Jian L. Zhou
  139. c#################################################################
  140. c
  141.       LQL=.FALSE.
  142.       IF (IWAR(1).EQ.1) LQL=.TRUE.
  143.       ZERO=0.0D+0
  144.       TEN=1.D+1
  145.       MAXIT=40*(M+N)
  146.       QPEPS=EPS
  147.       INW1=1
  148.       INW2=INW1+M
  149. C
  150. C     PREPARE PROBLEM DATA FOR EXECUTION
  151. C
  152.       IF (M.LE.0) GOTO 20
  153.       IN=INW1
  154.       DO 10 J=1,M
  155.       WAR(IN)=-B(J)
  156.    10 IN=IN+1
  157.    20 LW=NMAX*(3*NMAX+15)/2 + M
  158.       IF (INW2+LW-1 .GT. LWAR) GOTO 80
  159.       IF (LIWAR.LT.N) GOTO 81
  160.       IF (MNN.LT.M+N+N) GOTO 82
  161.       MN=M+N
  162. C
  163. C     CALL OF QL0002
  164. C
  165.       CALL QL0002(N,M,ME,MMAX,MN,MNN,NMAX,LQL,A,WAR(INW1),
  166.      1    D,C,XL,XU,X,NACT,IWAR,MAXIT,QPEPS,INFO,DIAG,
  167.      2    WAR(INW2),LW)
  168. C
  169. C     TEST OF MATRIX CORRECTIONS
  170. C
  171.       IFAIL=0
  172.       IF (INFO.EQ.1) GOTO 40
  173.       IF (INFO.EQ.2) GOTO 90
  174.       IDIAG=0
  175.       IF ((DIAG.GT.ZERO).AND.(DIAG.LT.1000.0)) IDIAG=DIAG
  176.       IF ((IPRINT.GT.0).AND.(IDIAG.GT.0))
  177.      1   WRITE(IOUT,1000) IDIAG
  178.       IF (INFO .LT. 0) GOTO  70
  179. C
  180. C     REORDER MULTIPLIER
  181. C
  182.       DO  50 J=1,MNN
  183.    50 U(J)=ZERO
  184.       IN=INW2-1
  185.       IF (NACT.EQ.0) GOTO 30
  186.       DO  60 I=1,NACT
  187.       J=IWAR(I)
  188.       U(J)=WAR(IN+I)
  189.    60 CONTINUE
  190.    30 CONTINUE
  191.       RETURN
  192. C
  193. C     ERROR MESSAGES
  194. C
  195.    70 IFAIL=-INFO+10
  196.       IF ((IPRINT.GT.0).AND.(NACT.GT.0))
  197.      1     WRITE(IOUT,1100) -INFO,(IWAR(I),I=1,NACT)
  198.       RETURN
  199.    80 IFAIL=5
  200.       IF (IPRINT .GT. 0) WRITE(IOUT,1200)
  201.       RETURN
  202.    81 IFAIL=5
  203.       IF (IPRINT .GT. 0) WRITE(IOUT,1210)
  204.       RETURN
  205.    82 IFAIL=5
  206.       IF (IPRINT .GT. 0) WRITE(IOUT,1220)
  207.       RETURN
  208.    40 IFAIL=1
  209.       IF (IPRINT.GT.0) WRITE(IOUT,1300) MAXIT
  210.       RETURN
  211.    90 IFAIL=2
  212.       IF (IPRINT.GT.0) WRITE(IOUT,1400) 
  213.       RETURN
  214. C
  215. C     FORMAT-INSTRUCTIONS
  216. C
  217.  1000 FORMAT(/8X,28H***QL: MATRIX G WAS ENLARGED,I3,
  218.      1        20H-TIMES BY UNITMATRIX)
  219.  1100 FORMAT(/8X,18H***QL: CONSTRAINT ,I5,
  220.      1        19H NOT CONSISTENT TO ,/,(10X,10I5))
  221.  1200 FORMAT(/8X,21H***QL: LWAR TOO SMALL)
  222.  1210 FORMAT(/8X,22H***QL: LIWAR TOO SMALL)
  223.  1220 FORMAT(/8X,20H***QL: MNN TOO SMALL)
  224.  1300 FORMAT(/8X,37H***QL: TOO MANY ITERATIONS (MORE THAN,I6,1H))
  225.  1400 FORMAT(/8X,50H***QL: ACCURACY INSUFFICIENT TO ATTAIN CONVERGENCE) 
  226.       END
  227.