SUBROUTINE LDL (N,A,Z,SIGMA,W)
C   LDL     LDL' - RANK-ONE - UPDATE

C   PURPOSE:
C           UPDATES THE LDL' FACTORS OF MATRIX A BY RANK-ONE MATRIX
C           SIGMA*Z*Z'

C   INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION)
C     N     : ORDER OF THE COEFFICIENT MATRIX A
C   * A     : POSITIVE DEFINITE MATRIX OF DIMENSION N;
C             ONLY THE LOWER TRIANGLE IS USED AND IS STORED COLUMN BY
C             COLUMN AS ONE DIMENSIONAL ARRAY OF DIMENSION N*(N+1)/2.
C   * Z     : VECTOR OF DIMENSION N OF UPDATING ELEMENTS
C     SIGMA : SCALAR FACTOR BY WHICH THE MODIFYING DYADE Z*Z' IS
C             MULTIPLIED

C   OUTPUT ARGUMENTS:
C     A     : UPDATED LDL' FACTORS

C   WORKING ARRAY:
C     W     : VECTOR OP DIMENSION N (USED ONLY IF SIGMA .LT. ZERO)

C   METHOD:
C     THAT OF FLETCHER AND POWELL AS DESCRIBED IN :
C     FLETCHER,R.,(1974) ON THE MODIFICATION OF LDL' FACTORIZATION.
C     POWELL,M.J.D.      MATH.COMPUTATION 28, 1067-1078.

C   IMPLEMENTED BY:
C     KRAFT,D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME
C               D-8031  OBERPFAFFENHOFEN

C   STATUS: 15. JANUARY 1980

C   SUBROUTINES REQUIRED: NONE

      INTEGER          I, IJ, J, N
      DOUBLE PRECISION A(*), T, V, W(*), Z(*), U, TP, ONE, BETA, FOUR,
     *                 ZERO, ALPHA, DELTA, GAMMA, SIGMA, EPMACH
      DATA ZERO, ONE, FOUR, EPMACH /0.0D0, 1.0D0, 4.0D0, 2.22D-16/

      IF(SIGMA.EQ.ZERO)               GOTO 280
      IJ=1
      T=ONE/SIGMA
      IF(SIGMA.GT.ZERO)               GOTO 220
C PREPARE NEGATIVE UPDATE
      DO 150 I=1,N
  150     W(I)=Z(I)
      DO 170 I=1,N
          V=W(I)
          T=T+V*V/A(IJ)
          DO 160 J=I+1,N
              IJ=IJ+1
  160         W(J)=W(J)-V*A(IJ)
  170     IJ=IJ+1
      IF(T.GE.ZERO) T=EPMACH/SIGMA
      DO 210 I=1,N
          J=N+1-I
          IJ=IJ-I
          U=W(J)
          W(J)=T
  210     T=T-U*U/A(IJ)
  220 CONTINUE
C HERE UPDATING BEGINS
      DO 270 I=1,N
          V=Z(I)
          DELTA=V/A(IJ)
          IF(SIGMA.LT.ZERO) TP=W(I)
          IF(SIGMA.GT.ZERO) TP=T+DELTA*V
          ALPHA=TP/T
          A(IJ)=ALPHA*A(IJ)
          IF(I.EQ.N)                  GOTO 280
          BETA=DELTA/TP
          IF(ALPHA.GT.FOUR)           GOTO 240
          DO 230 J=I+1,N
              IJ=IJ+1
              Z(J)=Z(J)-V*A(IJ)
  230         A(IJ)=A(IJ)+BETA*Z(J)
                                      GOTO 260
  240     GAMMA=T/TP
          DO 250 J=I+1,N
              IJ=IJ+1
              U=A(IJ)
              A(IJ)=GAMMA*U+BETA*Z(J)
  250         Z(J)=Z(J)-V*U
  260     IJ=IJ+1
  270     T=TP
  280 RETURN
C END OF LDL
      END

      DOUBLE PRECISION FUNCTION LINMIN (MODE, AX, BX, F, TOL)
C   LINMIN  LINESEARCH WITHOUT DERIVATIVES

C   PURPOSE:

C  TO FIND THE ARGUMENT LINMIN WHERE THE FUNCTION F TAKES IT'S MINIMUM
C  ON THE INTERVAL AX, BX.
C  COMBINATION OF GOLDEN SECTION AND SUCCESSIVE QUADRATIC INTERPOLATION.

C   INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION)

C *MODE   SEE OUTPUT ARGUMENTS
C  AX     LEFT ENDPOINT OF INITIAL INTERVAL
C  BX     RIGHT ENDPOINT OF INITIAL INTERVAL
C  F      FUNCTION VALUE AT LINMIN WHICH IS TO BE BROUGHT IN BY
C         REVERSE COMMUNICATION CONTROLLED BY MODE
C  TOL    DESIRED LENGTH OF INTERVAL OF UNCERTAINTY OF FINAL RESULT

C   OUTPUT ARGUMENTS:

C  LINMIN ABSCISSA APPROXIMATING THE POINT WHERE F ATTAINS A MINIMUM
C  MODE   CONTROLS REVERSE COMMUNICATION
C         MUST BE SET TO 0 INITIALLY, RETURNS WITH INTERMEDIATE
C         VALUES 1 AND 2 WHICH MUST NOT BE CHANGED BY THE USER,
C         ENDS WITH CONVERGENCE WITH VALUE 3.

C   WORKING ARRAY:

C  NONE

C   METHOD:

C  THIS FUNCTION SUBPROGRAM IS A SLIGHTLY MODIFIED VERSION OF THE
C  ALGOL 60 PROCEDURE LOCALMIN GIVEN IN
C  R.P. BRENT: ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES,
C              PRENTICE-HALL (1973).

C   IMPLEMENTED BY:

C     KRAFT, D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME
C                D-8031  OBERPFAFFENHOFEN

C   STATUS: 31. AUGUST  1984

C   SUBROUTINES REQUIRED: NONE

      INTEGER          MODE
      DOUBLE PRECISION F, TOL, A, B, C, D, E, P, Q, R, U, V, W, X, M,
     &                 FU, FV, FW, FX, EPS, TOL1, TOL2, ZERO, AX, BX
      DATA             C /0.381966011D0/, EPS /1.5D-8/, ZERO /0.0D0/

C  EPS = SQUARE - ROOT OF MACHINE PRECISION
C  C = GOLDEN SECTION RATIO = (3-SQRT(5))/2

      GOTO (10, 55), MODE

C  INITIALIZATION

      A = AX
      B = BX
      E = ZERO
      V = A + C*(B - A)
      W = V
      X = W
      LINMIN = X
      MODE = 1
      GOTO 100

C  MAIN LOOP STARTS HERE

   10 FX = F
      FV = FX
      FW = FV
   20 M = 0.5D0*(A + B)
      TOL1 = EPS*ABS(X) + TOL
      TOL2 = TOL1 + TOL1

C  TEST CONVERGENCE

      IF (ABS(X - M) .LE. TOL2 - 0.5D0*(B - A)) GOTO 90
      R = ZERO
      Q = R
      P = Q
      IF (ABS(E) .LE. TOL1) GOTO 30

C  FIT PARABOLA

      R = (X - W)*(FX - FV)
      Q = (X - V)*(FX - FW)
      P = (X - V)*Q - (X - W)*R
      Q = Q - R
      Q = Q + Q
      IF (Q .GT. ZERO) P = -P
      IF (Q .LT. ZERO) Q = -Q
      R = E
      E = D

C  IS PARABOLA ACCEPTABLE

   30 IF (ABS(P) .GE. 0.5D0*ABS(Q*R) .OR.
     &    P .LE. Q*(A - X) .OR. P .GE. Q*(B-X)) GOTO 40

C  PARABOLIC INTERPOLATION STEP

      D = P/Q

C  F MUST NOT BE EVALUATED TOO CLOSE TO A OR B

      IF (U - A .LT. TOL2) D = SIGN(TOL1, M - X)
      IF (B - U .LT. TOL2) D = SIGN(TOL1, M - X)
      GOTO 50

C  GOLDEN SECTION STEP

   40 IF (X .GE. M) E = A - X
      IF (X .LT. M) E = B - X
      D = C*E

C  F MUST NOT BE EVALUATED TOO CLOSE TO X

   50 IF (ABS(D) .LT. TOL1) D = SIGN(TOL1, D)
      U = X + D
      LINMIN = U
      MODE = 2
      GOTO 100
   55 FU = F

C  UPDATE A, B, V, W, AND X

      IF (FU .GT. FX) GOTO 60
      IF (U .GE. X) A = X
      IF (U .LT. X) B = X
      V = W
      FV = FW
      W = X
      FW = FX
      X = U
      FX = FU
      GOTO 85
   60 IF (U .LT. X) A = U
      IF (U .GE. X) B = U
      IF (FU .LE. FW .OR. W .EQ. X) GOTO 70
      IF (FU .LE. FV .OR. V .EQ. X .OR. V .EQ. W) GOTO 80
      GOTO 85
   70 V = W
      FV = FW
      W = U
      FW = FU
      GOTO 85
   80 V = U
      FV = FU
   85 GOTO 20

C  END OF MAIN LOOP

   90 LINMIN = X
      MODE = 3
  100 RETURN

C  END OF LINMIN

      END