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