      SUBROUTINE SLSQPB (M, MEQ, LA, N, X, XL, XU, F, C, G, A, ACC,
     *                   ITER, MODE, R, L, X0, MU, S, U, V, W, IW,
     *                   ALPHA, F0, GS, H1, H2, H3, H4, T, T0, TOL,
     *                   IEXACT, INCONS, IRESET, ITERMX, LINE, N1, 
     *                   N2, N3)

C   NONLINEAR PROGRAMMING BY SOLVING SEQUENTIALLY QUADRATIC PROGRAMS

C        -  L1 - LINE SEARCH,  POSITIVE DEFINITE  BFGS UPDATE  -

C                      BODY SUBROUTINE FOR SLSQP

      INTEGER          IW(*), I, IEXACT, INCONS, IRESET, ITER, ITERMX,
     *                 K, J, LA, LINE, M, MEQ, MODE, N, N1, N2, N3

      DOUBLE PRECISION A(LA,N+1), C(LA), G(N+1), L((N+1)*(N+2)/2),
     *                 MU(LA), R(M+N+N+2), S(N+1), U(N+1), V(N+1), W(*),
     *                 X(N), XL(N), XU(N), X0(N),
     *                 DDOT, DNRM2, LINMIN,
     *                 ACC, ALFMIN, ALPHA, F, F0, GS, H1, H2, H3, H4,
     *                 HUN, ONE, T, T0, TEN, TOL, TWO, ZERO

c     dim(W) =         N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1)  for LSQ
c                     +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ
c                     +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1       for LSEI
c                      with MINEQ = M - MEQ + 2*N1  &  N1 = N+1

      DATA             ZERO /0.0D0/, ONE /1.0D0/, ALFMIN /1.0D-1/,
     *                 HUN /1.0D+2/, TEN /1.0D+1/, TWO /2.0D0/

      IF (MODE) 260, 100, 220

  100 ITERMX = ITER
      IF (ACC.GE.ZERO) THEN
          IEXACT = 0
      ELSE
          IEXACT = 1
      ENDIF
      ACC = ABS(ACC)
      TOL = TEN*ACC
      ITER = 0
      IRESET = 0
      N1 = N + 1
      N2 = N1*N/2
      N3 = N2 + 1
      S(1) = ZERO
      MU(1) = ZERO
      CALL DCOPY(N, S(1),  0, S,  1)
      CALL DCOPY(M, MU(1), 0, MU, 1)

C   RESET BFGS MATRIX

  110 IRESET = IRESET + 1
      IF (IRESET.GT.5) GO TO 255
      L(1) = ZERO
      CALL DCOPY(N2, L(1), 0, L, 1)
      J = 1
      DO 120 I=1,N
         L(J) = ONE
         J = J + N1 - I
  120 CONTINUE

C   MAIN ITERATION : SEARCH DIRECTION, STEPLENGTH, LDL'-UPDATE

  130 ITER = ITER + 1
      MODE = 9
      IF (ITER.GT.ITERMX) GO TO 330

C   SEARCH DIRECTION AS SOLUTION OF QP - SUBPROBLEM

      CALL DCOPY(N, XL, 1, U, 1)
      CALL DCOPY(N, XU, 1, V, 1)
      CALL DAXPY(N, -ONE, X, 1, U, 1)
      CALL DAXPY(N, -ONE, X, 1, V, 1)
      H4 = ONE
      CALL LSQ (M, MEQ, N , N3, LA, L, G, A, C, U, V, S, R, W, IW, MODE)

C   AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION

      IF (MODE.EQ.6) THEN
          IF (N.EQ.MEQ) THEN
              MODE = 4
          ENDIF
      ENDIF
      IF (MODE.EQ.4) THEN
          DO 140 J=1,M
             IF (J.LE.MEQ) THEN
                 A(J,N1) = -C(J)
             ELSE
                 A(J,N1) = MAX(-C(J),ZERO)
             ENDIF
  140     CONTINUE
          S(1) = ZERO
          CALL DCOPY(N, S(1), 0, S, 1)
          H3 = ZERO
          G(N1) = ZERO
          L(N3) = HUN
          S(N1) = ONE
          U(N1) = ZERO
          V(N1) = ONE
          INCONS = 0
  150     CALL LSQ (M, MEQ, N1, N3, LA, L, G, A, C, U, V, S, R,
     *              W, IW, MODE)
          H4 = ONE - S(N1)
          IF (MODE.EQ.4) THEN
              L(N3) = TEN*L(N3)
              INCONS = INCONS + 1
              IF (INCONS.GT.5) GO TO 330
              GOTO 150
          ELSE IF (MODE.NE.1) THEN
              GOTO 330
          ENDIF
      ELSE IF (MODE.NE.1) THEN
          GOTO 330
      ENDIF

C   UPDATE MULTIPLIERS FOR L1-TEST

      DO 160 I=1,N
         V(I) = G(I) - DDOT(M,A(1,I),1,R,1)
  160 CONTINUE
      F0 = F
      CALL DCOPY(N, X, 1, X0, 1)
      GS = DDOT(N, G, 1, S, 1)
      H1 = ABS(GS)
      H2 = ZERO
      DO 170 J=1,M
         IF (J.LE.MEQ) THEN
             H3 = C(J)
         ELSE
             H3 = ZERO
         ENDIF
         H2 = H2 + MAX(-C(J),H3)
         H3 = ABS(R(J))
         MU(J) = MAX(H3,(MU(J)+H3)/TWO)
         H1 = H1 + H3*ABS(C(J))
  170 CONTINUE

C   CHECK CONVERGENCE

      MODE = 0
      IF (H1.LT.ACC .AND. H2.LT.ACC) GO TO 330
      H1 = ZERO
      DO 180 J=1,M
         IF (J.LE.MEQ) THEN
             H3 = C(J)
         ELSE
             H3 = ZERO
         ENDIF
         H1 = H1 + MU(J)*MAX(-C(J),H3)
  180 CONTINUE
      T0 = F + H1
      H3 = GS - H1*H4
      MODE = 8
      IF (H3.GE.ZERO) GO TO 110

C   LINE SEARCH WITH AN L1-TESTFUNCTION

      LINE = 0
      ALPHA = ONE
      IF (IEXACT.EQ.1) GOTO 210

C   INEXACT LINESEARCH

  190     LINE = LINE + 1
          H3 = ALPHA*H3
          CALL DSCAL(N, ALPHA, S, 1)
          CALL DCOPY(N, X0, 1, X, 1)
          CALL DAXPY(N, ONE, S, 1, X, 1)
          MODE = 1
          GO TO 330
  200         IF (H1.LE.H3/TEN .OR. LINE.GT.10) GO TO 240
              ALPHA = MAX(H3/(TWO*(H3-H1)),ALFMIN)
              GO TO 190

C   EXACT LINESEARCH

  210 IF (LINE.NE.3) THEN
          ALPHA = LINMIN(LINE,ALFMIN,ONE,T,TOL)
          CALL DCOPY(N, X0, 1, X, 1)
          CALL DAXPY(N, ALPHA, S, 1, X, 1)
          MODE = 1
          GOTO 330
      ENDIF
      CALL DSCAL(N, ALPHA, S, 1)
      GOTO 240

C   CALL FUNCTIONS AT CURRENT X

  220     T = F
          DO 230 J=1,M
             IF (J.LE.MEQ) THEN
                 H1 = C(J)
             ELSE
                 H1 = ZERO
             ENDIF
             T = T + MU(J)*MAX(-C(J),H1)
  230     CONTINUE
          H1 = T - T0
          GOTO (200, 210) IEXACT+1

C   CHECK CONVERGENCE

  240 H3 = ZERO
      DO 250 J=1,M
         IF (J.LE.MEQ) THEN
             H1 = C(J)
         ELSE
             H1 = ZERO
         ENDIF
         H3 = H3 + MAX(-C(J),H1)
  250 CONTINUE
      IF ((ABS(F-F0).LT.ACC .OR. DNRM2(N,S,1).LT.ACC) .AND. H3.LT.ACC)
     *   THEN
            MODE = 0
         ELSE
            MODE = -1
         ENDIF
      GO TO 330

C   CHECK relaxed CONVERGENCE in case of positive directional derivative

  255 CONTINUE
      IF ((ABS(F-F0).LT.TOL .OR. DNRM2(N,S,1).LT.TOL) .AND. H3.LT.TOL)
     *   THEN
            MODE = 0
         ELSE
            MODE = 8
         ENDIF
      GO TO 330

C   CALL JACOBIAN AT CURRENT X

C   UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA

  260 DO 270 I=1,N
         U(I) = G(I) - DDOT(M,A(1,I),1,R,1) - V(I)
  270 CONTINUE

C   L'*S

      K = 0
      DO 290 I=1,N
         H1 = ZERO
         K = K + 1
         DO 280 J=I+1,N
            K = K + 1
            H1 = H1 + L(K)*S(J)
  280    CONTINUE
         V(I) = S(I) + H1
  290 CONTINUE

C   D*L'*S

      K = 1
      DO 300 I=1,N
         V(I) = L(K)*V(I)
         K = K + N1 - I
  300 CONTINUE

C   L*D*L'*S

      DO 320 I=N,1,-1
         H1 = ZERO
         K = I
         DO 310 J=1,I - 1
            H1 = H1 + L(K)*V(J)
            K = K + N - J
  310    CONTINUE
         V(I) = V(I) + H1
  320 CONTINUE

      H1 = DDOT(N,S,1,U,1)
      H2 = DDOT(N,S,1,V,1)
      H3 = 0.2D0*H2
      IF (H1.LT.H3) THEN
          H4 = (H2-H3)/(H2-H1)
          H1 = H3
          CALL DSCAL(N, H4, U, 1)
          CALL DAXPY(N, ONE-H4, V, 1, U, 1)
      ENDIF
      CALL LDL(N, L, U, +ONE/H1, V)
      CALL LDL(N, L, V, -ONE/H2, U)

C   END OF MAIN ITERATION

      GO TO 130

C   END OF SLSQPB

  330 END
      