1 SUBROUTINE SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT, 2 + RNF,IAGFLG,NOUPDT,S,Y,T) 3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 4C 5C PURPOSE 6C ------- 7C UPDATE HESSIAN BY THE BFGS UNFACTORED METHOD 8C 9C PARAMETERS 10C ---------- 11C NR --> ROW DIMENSION OF MATRIX 12C N --> DIMENSION OF PROBLEM 13C X(N) --> OLD ITERATE, X[K-1] 14C G(N) --> GRADIENT OR APPROXIMATE AT OLD ITERATE 15C A(N,N) <--> ON ENTRY: APPROXIMATE HESSIAN AT OLD ITERATE 16C IN UPPER TRIANGULAR PART (AND UDIAG) 17C ON EXIT: UPDATED APPROX HESSIAN AT NEW ITERATE 18C IN LOWER TRIANGULAR PART AND DIAGONAL 19C [LOWER TRIANGULAR PART OF SYMMETRIC MATRIX] 20C UDIAG --> ON ENTRY: DIAGONAL OF HESSIAN 21C XPLS(N) --> NEW ITERATE, X[K] 22C GPLS(N) --> GRADIENT OR APPROXIMATE AT NEW ITERATE 23C EPSM --> MACHINE EPSILON 24C ITNCNT --> ITERATION COUNT 25C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN 26C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 OTHERWISE 27C NOUPDT <--> BOOLEAN: NO UPDATE YET 28C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] 29C S(N) --> WORKSPACE 30C Y(N) --> WORKSPACE 31C T(N) --> WORKSPACE 32C 33C***REVISION HISTORY (YYMMDD) 34C 000330 Modified array declarations. (JEC) 35C 36 DIMENSION X(N),G(N),XPLS(N),GPLS(N) 37 DIMENSION A(NR,*) 38 DIMENSION UDIAG(N) 39 DIMENSION S(N),Y(N),T(N) 40 LOGICAL NOUPDT,SKPUPD 41C 42C COPY HESSIAN IN UPPER TRIANGULAR PART AND UDIAG TO 43C LOWER TRIANGULAR PART AND DIAGONAL 44C 45 DO 5 J=1,N 46 A(J,J)=UDIAG(J) 47 IF(J.EQ.N) GO TO 5 48 JP1=J+1 49 DO 4 I=JP1,N 50 A(I,J)=A(J,I) 51 4 CONTINUE 52 5 CONTINUE 53C 54 IF(ITNCNT.EQ.1) NOUPDT=.TRUE. 55 DO 10 I=1,N 56 S(I)=XPLS(I)-X(I) 57 Y(I)=GPLS(I)-G(I) 58 10 CONTINUE 59 DEN1=DDOT(N,S,1,Y,1) 60 SNORM2=DNRM2(N,S,1) 61 YNRM2=DNRM2(N,Y,1) 62 IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 100 63C IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2) 64C THEN 65 CALL MVMLTS(NR,N,A,S,T) 66 DEN2=DDOT(N,S,1,T,1) 67 IF(.NOT. NOUPDT) GO TO 50 68C IF(NOUPDT) 69C THEN 70C 71C H <-- [(S+)Y/(S+)HS]H 72C 73 GAM=DEN1/DEN2 74 DEN2=GAM*DEN2 75 DO 30 J=1,N 76 T(J)=GAM*T(J) 77 DO 20 I=J,N 78 A(I,J)=GAM*A(I,J) 79 20 CONTINUE 80 30 CONTINUE 81 NOUPDT=.FALSE. 82C ENDIF 83 50 SKPUPD=.TRUE. 84C 85C CHECK UPDATE CONDITION ON ROW I 86C 87 DO 60 I=1,N 88 TOL=RNF*MAX(ABS(G(I)),ABS(GPLS(I))) 89 IF(IAGFLG.EQ.0) TOL=TOL/SQRT(RNF) 90 IF(ABS(Y(I)-T(I)).LT.TOL) GO TO 60 91C IF(ABS(Y(I)-T(I)).GE.TOL) 92C THEN 93 SKPUPD=.FALSE. 94 GO TO 70 95C ENDIF 96 60 CONTINUE 97 70 IF(SKPUPD) GO TO 100 98C IF(.NOT.SKPUPD) 99C THEN 100C 101C BFGS UPDATE 102C 103 DO 90 J=1,N 104 DO 80 I=J,N 105 A(I,J)=A(I,J)+Y(I)*Y(J)/DEN1-T(I)*T(J)/DEN2 106 80 CONTINUE 107 90 CONTINUE 108C ENDIF 109C ENDIF 110 100 RETURN 111 END 112