      SUBROUTINE QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,
     1           X,U,IOUT,IFAIL,IPRINT,WAR,LWAR,IWAR,LIWAR)
c
cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
c                     !!!! NOTICE !!!!
c
c 1. The routines contained in this file are due to Prof. K.Schittkowski
c    of the University of Bayreuth, Germany (modification of routines
c    due to Prof. MJD Powell at the University of Cambridge).  They can
c    be freely distributed.
c
c 2. A minor modification was performed at the University of Maryland. 
c    It is marked in the code by "c umd".
c
c                                      A.L. Tits and J.L. Zhou
c                                      University of Maryland
C
C***********************************************************************
C
C
C             SOLUTION OF QUADRATIC PROGRAMMING PROBLEMS
C
C
C
C   QL0001 SOLVES THE QUADRATIC PROGRAMMING PROBLEM
C
C   MINIMIZE        .5*X'*C*X + D'*X
C   SUBJECT TO      A(J)*X  +  B(J)   =  0  ,  J=1,...,ME
C                   A(J)*X  +  B(J)  >=  0  ,  J=ME+1,...,M
C                   XL  <=  X  <=  XU
C   
C   HERE C MUST BE AN N BY N SYMMETRIC AND POSITIVE MATRIX, D AN N-DIMENSIONAL
C   VECTOR, A AN M BY N MATRIX AND B AN M-DIMENSIONAL VECTOR. THE ABOVE 
C   SITUATION IS INDICATED BY IWAR(1)=1. ALTERNATIVELY, I.E. IF IWAR(1)=0,
C   THE OBJECTIVE FUNCTION MATRIX CAN ALSO BE PROVIDED IN FACTORIZED FORM.
C   IN THIS CASE, C IS AN UPPER TRIANGULAR MATRIX.
C
C   THE SUBROUTINE REORGANIZES SOME DATA SO THAT THE PROBLEM CAN BE SOLVED
C   BY A MODIFICATION OF AN ALGORITHM PROPOSED BY POWELL (1983).
C
C
C   USAGE:
C
C      QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,X,U,IOUT,IFAIL,IPRINT,
C             WAR,LWAR,IWAR,LIWAR)
C
C
C   DEFINITION OF THE PARAMETERS:
C
C   M :        TOTAL NUMBER OF CONSTRAINTS.
C   ME :       NUMBER OF EQUALITY CONSTRAINTS.
C   MMAX :     ROW DIMENSION OF A. MMAX MUST BE AT LEAST ONE AND GREATER
C              THAN M.
C   N :        NUMBER OF VARIABLES.
C   NMAX :     ROW DIMENSION OF C. NMAX MUST BE GREATER OR EQUAL TO N.
C   MNN :      MUST BE EQUAL TO M + N + N.
C   C(NMAX,NMAX): OBJECTIVE FUNCTION MATRIX WHICH SHOULD BE SYMMETRIC AND
C              POSITIVE DEFINITE. IF IWAR(1) = 0, C IS SUPPOSED TO BE THE
C              CHOLESKEY-FACTOR OF ANOTHER MATRIX, I.E. C IS UPPER
C              TRIANGULAR.
C   D(NMAX) :  CONTAINS THE CONSTANT VECTOR OF THE OBJECTIVE FUNCTION.
C   A(MMAX,NMAX): CONTAINS THE DATA MATRIX OF THE LINEAR CONSTRAINTS.
C   B(MMAX) :  CONTAINS THE CONSTANT DATA OF THE LINEAR CONSTRAINTS.
C   XL(N),XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR THE VARIABLES.
C   X(N) :     ON RETURN, X CONTAINS THE OPTIMAL SOLUTION VECTOR.
C   U(MNN) :   ON RETURN, U CONTAINS THE LAGRANGE MULTIPLIERS. THE FIRST
C              M POSITIONS ARE RESERVED FOR THE MULTIPLIERS OF THE M
C              LINEAR CONSTRAINTS AND THE SUBSEQUENT ONES FOR THE 
C              MULTIPLIERS OF THE LOWER AND UPPER BOUNDS. ON SUCCESSFUL
C              TERMINATION, ALL VALUES OF U WITH RESPECT TO INEQUALITIES 
C              AND BOUNDS SHOULD BE GREATER OR EQUAL TO ZERO.
C   IOUT :     INTEGER INDICATING THE DESIRED OUTPUT UNIT NUMBER, I.E.
C              ALL WRITE-STATEMENTS START WITH 'WRITE(IOUT,... '.
C   IFAIL :    SHOWS THE TERMINATION REASON.
C      IFAIL = 0 :   SUCCESSFUL RETURN.
C      IFAIL = 1 :   TOO MANY ITERATIONS (MORE THAN 40*(N+M)).
C      IFAIL = 2 :   ACCURACY INSUFFICIENT TO SATISFY CONVERGENCE
C                    CRITERION.
C      IFAIL = 5 :   LENGTH OF A WORKING ARRAY IS TOO SHORT.
C      IFAIL > 10 :  THE CONSTRAINTS ARE INCONSISTENT.
C   IPRINT :   OUTPUT CONTROL.
C      IPRINT = 0 :  NO OUTPUT OF QL0001.
C      IPRINT > 0 :  BRIEF OUTPUT IN ERROR CASES.
C   WAR(LWAR) : REAL WORKING ARRAY. THE LENGTH LWAR SHOULD BE GRATER THAN
C               NMAX*(3*NMAX+15)/2 + 2*M.
C   IWAR(LIWAR): INTEGER WORKING ARRAY. THE LENGTH LIWAR SHOULD BE AT
C              LEAST N.
C              IF IWAR(1)=0 INITIALLY, THEN THE CHOLESKY DECOMPOSITION
C              WHICH IS REQUIRED BY THE DUAL ALGORITHM TO GET THE FIRST
C              UNCONSTRAINED MINIMUM OF THE OBJECTIVE FUNCTION, IS
C              PERFORMED INTERNALLY. OTHERWISE, I.E. IF IWAR(1)=1, THEN
C              IT IS ASSUMED THAT THE USER PROVIDES THE INITIAL FAC-
C              TORIZATION BY HIMSELF AND STORES IT IN THE UPPER TRIAN-
C              GULAR PART OF THE ARRAY C.
C
C   A NAMED COMMON-BLOCK  /CMACHE/EPS   MUST BE PROVIDED BY THE USER,
C   WHERE EPS DEFINES A GUESS FOR THE UNDERLYING MACHINE PRECISION.
C
C
C   AUTHOR:    K. SCHITTKOWSKI,
C              MATHEMATISCHES INSTITUT,
C              UNIVERSITAET BAYREUTH,
C              8580 BAYREUTH,
C              GERMANY, F.R.
C
C
C   VERSION:   1.4  (MARCH, 1987)  
C
C
C*********************************************************************
C
C
      INTEGER NMAX,MMAX,N,MNN,LWAR,LIWAR
      DIMENSION C(NMAX,NMAX),D(NMAX),A(MMAX,NMAX),B(MMAX),
     1      XL(N),XU(N),X(N),U(MNN),WAR(LWAR),IWAR(LIWAR)
      DOUBLE PRECISION C,D,A,B,X,XL,XU,U,WAR,DIAG,ZERO,
     1      EPS,QPEPS,TEN
      INTEGER M,ME,IOUT,IFAIL,IPRINT,IWAR,INW1,INW2,IN,J,LW,MN,I,         
     1      IDIAG,INFO,NACT,MAXIT
      LOGICAL LQL
C
C     INTRINSIC FUNCTIONS:  DSQRT
C
      COMMON /CMACHE/EPS
C
C     CONSTANT DATA
C
c#################################################################
c

      if(c(nmax,nmax).eq.0.d0) c(nmax,nmax)=eps
c
c umd
c  This prevents a subsequent more major modification of the Hessian
c  matrix in the important case when a minmax problem (yielding a 
c  singular Hessian matrix) is being solved.
c                                 ----UMCP, April 1991, Jian L. Zhou
c#################################################################
c
      LQL=.FALSE.
      IF (IWAR(1).EQ.1) LQL=.TRUE.
      ZERO=0.0D+0
      TEN=1.D+1
      MAXIT=40*(M+N)
      QPEPS=EPS
      INW1=1
      INW2=INW1+M
C
C     PREPARE PROBLEM DATA FOR EXECUTION
C
      IF (M.LE.0) GOTO 20
      IN=INW1
      DO 10 J=1,M
      WAR(IN)=-B(J)
   10 IN=IN+1
   20 LW=NMAX*(3*NMAX+15)/2 + M
      IF (INW2+LW-1 .GT. LWAR) GOTO 80
      IF (LIWAR.LT.N) GOTO 81
      IF (MNN.LT.M+N+N) GOTO 82
      MN=M+N
C
C     CALL OF QL0002
C
      CALL QL0002(N,M,ME,MMAX,MN,MNN,NMAX,LQL,A,WAR(INW1),
     1    D,C,XL,XU,X,NACT,IWAR,MAXIT,QPEPS,INFO,DIAG,
     2    WAR(INW2),LW)
C
C     TEST OF MATRIX CORRECTIONS
C
      IFAIL=0
      IF (INFO.EQ.1) GOTO 40
      IF (INFO.EQ.2) GOTO 90
      IDIAG=0
      IF ((DIAG.GT.ZERO).AND.(DIAG.LT.1000.0)) IDIAG=DIAG
      IF ((IPRINT.GT.0).AND.(IDIAG.GT.0))
     1   WRITE(IOUT,1000) IDIAG
      IF (INFO .LT. 0) GOTO  70
C
C     REORDER MULTIPLIER
C
      DO  50 J=1,MNN
   50 U(J)=ZERO
      IN=INW2-1
      IF (NACT.EQ.0) GOTO 30
      DO  60 I=1,NACT
      J=IWAR(I)
      U(J)=WAR(IN+I)
   60 CONTINUE
   30 CONTINUE
      RETURN
C
C     ERROR MESSAGES
C
   70 IFAIL=-INFO+10
      IF ((IPRINT.GT.0).AND.(NACT.GT.0))
     1     WRITE(IOUT,1100) -INFO,(IWAR(I),I=1,NACT)
      RETURN
   80 IFAIL=5
      IF (IPRINT .GT. 0) WRITE(IOUT,1200)
      RETURN
   81 IFAIL=5
      IF (IPRINT .GT. 0) WRITE(IOUT,1210)
      RETURN
   82 IFAIL=5
      IF (IPRINT .GT. 0) WRITE(IOUT,1220)
      RETURN
   40 IFAIL=1
      IF (IPRINT.GT.0) WRITE(IOUT,1300) MAXIT
      RETURN
   90 IFAIL=2
      IF (IPRINT.GT.0) WRITE(IOUT,1400) 
      RETURN
C
C     FORMAT-INSTRUCTIONS
C
 1000 FORMAT(/8X,28H***QL: MATRIX G WAS ENLARGED,I3,
     1        20H-TIMES BY UNITMATRIX)
 1100 FORMAT(/8X,18H***QL: CONSTRAINT ,I5,
     1        19H NOT CONSISTENT TO ,/,(10X,10I5))
 1200 FORMAT(/8X,21H***QL: LWAR TOO SMALL)
 1210 FORMAT(/8X,22H***QL: LIWAR TOO SMALL)
 1220 FORMAT(/8X,20H***QL: MNN TOO SMALL)
 1300 FORMAT(/8X,37H***QL: TOO MANY ITERATIONS (MORE THAN,I6,1H))
 1400 FORMAT(/8X,50H***QL: ACCURACY INSUFFICIENT TO ATTAIN CONVERGENCE) 
      END
