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