|
|
SUBROUTINE BAKSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE AX=B WHERE A IS UPPER TRIANGULAR MATRIX.
C NOTE THAT A IS INPUT AS A LOWER TRIANGULAR MATRIX AND
C THAT THIS ROUTINE TAKES ITS TRANSPOSE IMPLICITLY.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) --> LOWER TRIANGULAR MATRIX (PRESERVED)
C X(N) <-- SOLUTION VECTOR
C B(N) --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE,
C THEN VECTORS B AND X MAY SHARE THE SAME STORAGE.
C
DIMENSION A(NR,1),X(N),B(N)
C
C SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE)
C
I=N
X(I)=B(I)/A(I,I)
IF(N.EQ.1) RETURN
30 IP1=I
I=I-1
SUM=0.
DO 40 J=IP1,N
SUM=SUM+A(J,I)*X(J)
40 CONTINUE
X(I)=(B(I)-SUM)/A(I,I)
IF(I.GT.1) GO TO 30
RETURN
END
SUBROUTINE CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)
C
C PURPOSE
C -------
C FIND THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION
C OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO A IF
C NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) <--> ON ENTRY: MATRIX FOR WHICH TO FIND PERTURBED
C CHOLESKY DECOMPOSITION
C ON EXIT: CONTAINS L OF LL+ DECOMPOSITION
C IN LOWER TRIANGULAR PART AND DIAGONAL OF "A"
C DIAGMX --> MAXIMUM DIAGONAL ELEMENT OF "A"
C TOL --> TOLERANCE
C ADDMAX <-- MAXIMUM AMOUNT IMPLICITLY ADDED TO DIAGONAL OF "A"
C IN FORMING THE CHOLESKY DECOMPOSITION OF A+D
C INTERNAL VARIABLES
C ------------------
C AMINL SMALLEST ELEMENT ALLOWED ON DIAGONAL OF L
C AMNLSQ =AMINL**2
C OFFMAX MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN OF A
C
C
C DESCRIPTION
C -----------
C THE NORMAL CHOLESKY DECOMPOSITION IS PERFORMED. HOWEVER, IF AT ANY
C POINT THE ALGORITHM WOULD ATTEMPT TO SET L(I,I)=SQRT(TEMP)
C WITH TEMP < TOL*DIAGMX, THEN L(I,I) IS SET TO SQRT(TOL*DIAGMX)
C INSTEAD. THIS IS EQUIVALENT TO ADDING TOL*DIAGMX-TEMP TO A(I,I)
C
C
DIMENSION A(NR,1)
C
ADDMAX=0.
AMINL=SQRT(DIAGMX*TOL)
AMNLSQ=AMINL*AMINL
C
C FORM COLUMN J OF L
C
DO 100 J=1,N
C FIND DIAGONAL ELEMENTS OF L
SUM=0.
IF(J.EQ.1) GO TO 20
JM1=J-1
DO 10 K=1,JM1
SUM=SUM + A(J,K)*A(J,K)
10 CONTINUE
20 TEMP=A(J,J)-SUM
IF(TEMP.LT.AMNLSQ) GO TO 30
C IF(TEMP.GE.AMINL**2)
C THEN
A(J,J)=SQRT(TEMP)
GO TO 40
C ELSE
C
C FIND MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN
30 OFFMAX=0.
IF(J.EQ.N) GO TO 37
JP1=J+1
DO 35 I=JP1,N
IF(ABS(A(I,J)).GT.OFFMAX) OFFMAX=ABS(A(I,J))
35 CONTINUE
37 IF(OFFMAX.LE.AMNLSQ) OFFMAX=AMNLSQ
C
C ADD TO DIAGONAL ELEMENT TO ALLOW CHOLESKY DECOMPOSITION TO CONTINUE
A(J,J)=SQRT(OFFMAX)
ADDMAX=AMAX1(ADDMAX,OFFMAX-TEMP)
C ENDIF
C
C FIND I,J ELEMENT OF LOWER TRIANGULAR MATRIX
40 IF(J.EQ.N) GO TO 100
JP1=J+1
DO 70 I=JP1,N
SUM=0.0
IF(J.EQ.1) GO TO 60
JM1=J-1
DO 50 K=1,JM1
SUM=SUM+A(I,K)*A(J,K)
50 CONTINUE
60 A(I,J)=(A(I,J)-SUM)/A(J,J)
70 CONTINUE
100 CONTINUE
RETURN
END
SUBROUTINE CHLHSN(NR,N,A,EPSM,SX,UDIAG)
C
C PURPOSE
C -------
C FIND THE L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION OF THE PERTURBED
C MODEL HESSIAN MATRIX A+MU*I(WHERE MU\0 AND I IS THE IDENTITY MATRIX)
C WHICH IS SAFELY POSITIVE DEFINITE. IF A IS SAFELY POSITIVE DEFINITE
C UPON ENTRY, THEN MU=0.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) <--> ON ENTRY; "A" IS MODEL HESSIAN (ONLY LOWER
C TRIANGULAR PART AND DIAGONAL STORED)
C ON EXIT: A CONTAINS L OF LL+ DECOMPOSITION OF
C PERTURBED MODEL HESSIAN IN LOWER TRIANGULAR
C PART AND DIAGONAL AND CONTAINS HESSIAN IN UPPER
C TRIANGULAR PART AND UDIAG
C EPSM --> MACHINE EPSILON
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C UDIAG(N) <-- ON EXIT: CONTAINS DIAGONAL OF HESSIAN
C
C INTERNAL VARIABLES
C ------------------
C TOL TOLERANCE
C DIAGMN MINIMUM ELEMENT ON DIAGONAL OF A
C DIAGMX MAXIMUM ELEMENT ON DIAGONAL OF A
C OFFMAX MAXIMUM OFF-DIAGONAL ELEMENT OF A
C OFFROW SUM OF OFF-DIAGONAL ELEMENTS IN A ROW OF A
C EVMIN MINIMUM EIGENVALUE OF A
C EVMAX MAXIMUM EIGENVALUE OF A
C
C DESCRIPTION
C -----------
C 1. IF "A" HAS ANY NEGATIVE DIAGONAL ELEMENTS, THEN CHOOSE MU>0
C SUCH THAT THE DIAGONAL OF A:=A+MU*I IS ALL POSITIVE
C WITH THE RATIO OF ITS SMALLEST TO LARGEST ELEMENT ON THE
C ORDER OF SQRT(EPSM).
C
C 2. "A" UNDERGOES A PERTURBED CHOLESKY DECOMPOSITION WHICH
C RESULTS IN AN LL+ DECOMPOSITION OF A+D, WHERE D IS A
C NON-NEGATIVE DIAGONAL MATRIX WHICH IS IMPLICITLY ADDED TO
C "A" DURING THE DECOMPOSITION IF "A" IS NOT POSITIVE DEFINITE.
C "A" IS RETAINED AND NOT CHANGED DURING THIS PROCESS BY
C COPYING L INTO THE UPPER TRIANGULAR PART OF "A" AND THE
C DIAGONAL INTO UDIAG. THEN THE CHOLESKY DECOMPOSITION ROUTINE
C IS CALLED. ON RETURN, ADDMAX CONTAINS MAXIMUM ELEMENT OF D.
C
C 3. IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2
C AND RETURN IS MADE TO CALLING PROGRAM. OTHERWISE,
C THE MINIMUM NUMBER SDD WHICH MUST BE ADDED TO THE
C DIAGONAL OF A TO MAKE IT SAFELY STRICTLY DIAGONALLY DOMINANT
C IS CALCULATED. SINCE A+ADDMAX*I AND A+SDD*I ARE SAFELY
C POSITIVE DEFINITE, CHOOSE MU=MIN(ADDMAX,SDD) AND DECOMPOSE
C A+MU*I TO OBTAIN L.
C
DIMENSION A(NR,1),SX(N),UDIAG(N)
C
C SCALE HESSIAN
C PRE- AND POST- MULTIPLY "A" BY INV(SX)
C
DO 20 J=1,N
DO 10 I=J,N
A(I,J)=A(I,J)/(SX(I)*SX(J))
10 CONTINUE
20 CONTINUE
C
C STEP1
C -----
C NOTE: IF A DIFFERENT TOLERANCE IS DESIRED THROUGHOUT THIS
C ALGORITHM, CHANGE TOLERANCE HERE:
TOL=SQRT(EPSM)
C
DIAGMX=A(1,1)
DIAGMN=A(1,1)
IF(N.EQ.1) GO TO 35
DO 30 I=2,N
IF(A(I,I).LT.DIAGMN) DIAGMN=A(I,I)
IF(A(I,I).GT.DIAGMX) DIAGMX=A(I,I)
30 CONTINUE
35 POSMAX=AMAX1(DIAGMX,0.)
C
C DIAGMN .LE. 0
C
IF(DIAGMN.GT.POSMAX*TOL) GO TO 100
C IF(DIAGMN.LE.POSMAX*TOL)
C THEN
AMU=TOL*(POSMAX-DIAGMN)-DIAGMN
IF(AMU.NE.0.) GO TO 60
C IF(AMU.EQ.0.)
C THEN
C
C FIND LARGEST OFF-DIAGONAL ELEMENT OF A
OFFMAX=0.
IF(N.EQ.1) GO TO 50
DO 45 I=2,N
IM1=I-1
DO 40 J=1,IM1
IF(ABS(A(I,J)).GT.OFFMAX) OFFMAX=ABS(A(I,J))
40 CONTINUE
45 CONTINUE
50 AMU=OFFMAX
IF(AMU.NE.0.) GO TO 55
C IF(AMU.EQ.0.)
C THEN
AMU=1.0
GO TO 60
C ELSE
55 AMU=AMU*(1.0+TOL)
C ENDIF
C ENDIF
C
C A=A + MU*I
C
60 DO 65 I=1,N
A(I,I)=A(I,I)+AMU
65 CONTINUE
DIAGMX=DIAGMX+AMU
C ENDIF
C
C STEP2
C -----
C COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
C AND DIAGONAL OF "A" TO UDIAG
C
100 CONTINUE
DO 110 J=1,N
UDIAG(J)=A(J,J)
IF(J.EQ.N) GO TO 110
JP1=J+1
DO 105 I=JP1,N
A(J,I)=A(I,J)
105 CONTINUE
110 CONTINUE
C
CALL CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)
C
C
C STEP3
C -----
C IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2,
C THE LL+ DECOMPOSITION HAS BEEN DONE, AND WE RETURN.
C OTHERWISE, ADDMAX>0. PERTURB "A" SO THAT IT IS SAFELY
C DIAGONALLY DOMINANT AND FIND LL+ DECOMPOSITION
C
IF(ADDMAX.LE.0.) GO TO 170
C IF(ADDMAX.GT.0.)
C THEN
C
C RESTORE ORIGINAL "A" (LOWER TRIANGULAR PART AND DIAGONAL)
C
DO 120 J=1,N
A(J,J)=UDIAG(J)
IF(J.EQ.N) GO TO 120
JP1=J+1
DO 115 I=JP1,N
A(I,J)=A(J,I)
115 CONTINUE
120 CONTINUE
C
C FIND SDD SUCH THAT A+SDD*I IS SAFELY POSITIVE DEFINITE
C NOTE: EVMIN<0 SINCE A IS NOT POSITIVE DEFINITE;
C
EVMIN=0.
EVMAX=A(1,1)
DO 150 I=1,N
OFFROW=0.
IF(I.EQ.1) GO TO 135
IM1=I-1
DO 130 J=1,IM1
OFFROW=OFFROW+ABS(A(I,J))
130 CONTINUE
135 IF(I.EQ.N) GO TO 145
IP1=I+1
DO 140 J=IP1,N
OFFROW=OFFROW+ABS(A(J,I))
140 CONTINUE
145 EVMIN=AMIN1(EVMIN,A(I,I)-OFFROW)
EVMAX=AMAX1(EVMAX,A(I,I)+OFFROW)
150 CONTINUE
SDD=TOL*(EVMAX-EVMIN)-EVMIN
C
C PERTURB "A" AND DECOMPOSE AGAIN
C
AMU=AMIN1(SDD,ADDMAX)
DO 160 I=1,N
A(I,I)=A(I,I)+AMU
UDIAG(I)=A(I,I)
160 CONTINUE
C
C "A" NOW GUARANTEED SAFELY POSITIVE DEFINITE
C
CALL CHOLDC(NR,N,A,0.,TOL,ADDMAX)
C ENDIF
C
C UNSCALE HESSIAN AND CHOLESKY DECOMPOSITION MATRIX
C
170 DO 190 J=1,N
DO 175 I=J,N
A(I,J)=SX(I)*A(I,J)
175 CONTINUE
IF(J.EQ.1) GO TO 185
JM1=J-1
DO 180 I=1,JM1
A(I,J)=SX(I)*SX(J)*A(I,J)
180 CONTINUE
185 UDIAG(J)=UDIAG(J)*SX(J)*SX(J)
190 CONTINUE
RETURN
END
SUBROUTINE DFAULT(N,X,TYPSIZ,FSCALE,METHOD,IEXP,MSG,NDIGIT,
+ ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADTL,STEPMX,STEPTL)
C
C PURPOSE
C -------
C SET DEFAULT VALUES FOR EACH INPUT VARIABLE TO
C MINIMIZATION ALGORITHM.
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C X(N) --> INITIAL GUESS TO SOLUTION (TO COMPUTE MAX STEP SIZE)
C TYPSIZ(N) <-- TYPICAL SIZE FOR EACH COMPONENT OF X
C FSCALE <-- ESTIMATE OF SCALE OF MINIMIZATION FUNCTION
C METHOD <-- ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C IEXP <-- =0 IF MINIMIZATION FUNCTION NOT EXPENSIVE TO EVALUATE
C MSG <-- MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT
C NDIGIT <-- NUMBER OF GOOD DIGITS IN MINIMIZATION FUNCTION
C ITNLIM <-- MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IAGFLG <-- =0 IF ANALYTIC GRADIENT NOT SUPPLIED
C IAHFLG <-- =0 IF ANALYTIC HESSIAN NOT SUPPLIED
C IPR <-- DEVICE TO WHICH TO SEND OUTPUT
C DLT <-- TRUST REGION RADIUS
C GRADTL <-- TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE ENOUGH
C TO ZERO TO TERMINATE ALGORITHM
C STEPMX <-- VALUE OF ZERO TO TRIP DEFAULT MAXIMUM IN OPTCHK
C STEPTL <-- TOLERANCE AT WHICH SUCCESSIVE ITERATES CONSIDERED
C CLOSE ENOUGH TO TERMINATE ALGORITHM
C
DIMENSION TYPSIZ(N),X(N)
C
C SET TYPICAL SIZE OF X AND MINIMIZATION FUNCTION
DO 10 I=1,N
TYPSIZ(I)=1.0
10 CONTINUE
FSCALE=1.0
C
C SET TOLERANCES
DLT=-1.0
EPSM=EPSMCH(1.0)
GRADTL=EPSM**(1.0/3.0)
STEPMX=0.0
STEPTL=SQRT(EPSM)
C
C SET FLAGS
METHOD=1
IEXP=1
MSG=0
NDIGIT=-1
ITNLIM=150
IAGFLG=0
IAHFLG=0
IPR=6
C
RETURN
END
SUBROUTINE DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX,
+ STEPTL,DLT,IRETCD,MXTAKE,SC,WRK1,WRK2,WRK3,IPR)
C
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE (XPLS) BY THE DOUBLE DOGLEG METHOD
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE X[K-1]
C F --> FUNCTION VALUE AT OLD ITERATE, F(X)
C G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE
C A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN
C IN LOWER TRIANGULAR PART AND DIAGONAL
C P(N) --> NEWTON STEP
C XPLS(N) <-- NEW ITERATE X[K]
C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C DLT <--> TRUST REGION RADIUS
C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C IRETCD <-- RETURN CODE
C =0 SATISFACTORY XPLS FOUND
C =1 FAILED TO FIND SATISFACTORY XPLS SUFFICIENTLY
C DISTINCT FROM X
C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C SC(N) --> WORKSPACE [CURRENT STEP]
C WRK1(N) --> WORKSPACE (AND PLACE HOLDING ARGUMENT TO TREGUP)
C WRK2(N) --> WORKSPACE
C WRK3(N) --> WORKSPACE
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
DIMENSION X(N),XPLS(N),G(N),P(N)
DIMENSION SX(N)
DIMENSION SC(N),WRK1(N),WRK2(N),WRK3(N)
DIMENSION A(NR,1)
LOGICAL FSTDOG,NWTAKE,MXTAKE
EXTERNAL FCN
C
IRETCD=4
FSTDOG=.TRUE.
TMP=0.
DO 5 I=1,N
TMP=TMP+SX(I)*SX(I)*P(I)*P(I)
5 CONTINUE
RNWTLN=SQRT(TMP)
C$ WRITE(IPR,954) RNWTLN
C
100 CONTINUE
C
C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
CALL DOGSTP(NR,N,G,A,P,SX,RNWTLN,DLT,NWTAKE,FSTDOG,
+ WRK1,WRK2,CLN,ETA,SC,IPR,STEPMX)
C
C CHECK NEW POINT AND UPDATE TRUST REGION
CALL TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,DLT,
+ IRETCD,WRK3,FPLSP,XPLS,FPLS,MXTAKE,IPR,2,WRK1)
IF(IRETCD.LE.1) RETURN
GO TO 100
950 FORMAT(42H DOGDRV INITIAL TRUST REGION NOT GIVEN.,
+ 22H COMPUTE CAUCHY STEP.)
951 FORMAT(18H DOGDRV ALPHA =,E20.13/
+ 18H DOGDRV BETA =,E20.13/
+ 18H DOGDRV DLT =,E20.13/
+ 18H DOGDRV NWTAKE=,L1 )
952 FORMAT(28H DOGDRV CURRENT STEP (SC))
954 FORMAT(18H0DOGDRV RNWTLN=,E20.13)
955 FORMAT(14H DOGDRV ,5(E20.13,3X))
END
SUBROUTINE DOGSTP(NR,N,G,A,P,SX,RNWTLN,DLT,NWTAKE,FSTDOG,
+ SSD,V,CLN,ETA,SC,IPR,STEPMX)
C
C PURPOSE
C -------
C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C G(N) --> GRADIENT AT CURRENT ITERATE, G(X)
C A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN
C LOWER PART AND DIAGONAL
C P(N) --> NEWTON STEP
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C RNWTLN --> NEWTON STEP LENGTH
C DLT <--> TRUST REGION RADIUS
C NWTAKE <--> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C FSTDOG <--> BOOLEAN, =.TRUE. IF ON FIRST LEG OF DOGLEG
C SSD(N) <--> WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE
C QUADRATIC MODEL IN THE SCALED STEEPEST DESCENT
C DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C V(N) <--> WORKSPACE [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C CLN <--> CAUCHY LENGTH
C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C ETA [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C SC(N) <-- CURRENT STEP
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C
C INTERNAL VARIABLES
C ------------------
C CLN LENGTH OF CAUCHY STEP
C
DIMENSION G(N),P(N)
DIMENSION SX(N)
DIMENSION SC(N),SSD(N),V(N)
DIMENSION A(NR,1)
LOGICAL NWTAKE,FSTDOG
C
C CAN WE TAKE NEWTON STEP
C
IF(RNWTLN.GT.DLT) GO TO 100
C IF(RNWTLN.LE.DLT)
C THEN
NWTAKE=.TRUE.
DO 10 I=1,N
SC(I)=P(I)
10 CONTINUE
DLT=RNWTLN
C$ WRITE(IPR,951)
GO TO 700
C ELSE
C
C NEWTON STEP TOO LONG
C CAUCHY STEP IS ON DOUBLE DOGLEG CURVE
C
100 NWTAKE=.FALSE.
IF(.NOT.FSTDOG) GO TO 200
C IF(FSTDOG)
C THEN
C
C CALCULATE DOUBLE DOGLEG CURVE (SSD)
FSTDOG=.FALSE.
ALPHA=0.
DO 110 I=1,N
ALPHA=ALPHA + (G(I)*G(I))/(SX(I)*SX(I))
110 CONTINUE
BETA=0.
DO 130 I=1,N
TMP=0.
DO 120 J=I,N
TMP=TMP + (A(J,I)*G(J))/(SX(J)*SX(J))
120 CONTINUE
BETA=BETA+TMP*TMP
130 CONTINUE
DO 140 I=1,N
SSD(I)=-(ALPHA/BETA)*G(I)/SX(I)
140 CONTINUE
CLN=ALPHA*SQRT(ALPHA)/BETA
ETA=.2 + (.8*ALPHA*ALPHA)/(-BETA*SDOT(N,G,1,P,1))
DO 150 I=1,N
V(I)=ETA*SX(I)*P(I) - SSD(I)
150 CONTINUE
IF (DLT .EQ. (-1.0)) DLT = AMIN1(CLN, STEPMX)
C$ WRITE(IPR,954) ALPHA,BETA,CLN,ETA
C$ WRITE(IPR,955)
C$ WRITE(IPR,960) (SSD(I),I=1,N)
C$ WRITE(IPR,956)
C$ WRITE(IPR,960) (V(I),I=1,N)
C ENDIF
200 IF(ETA*RNWTLN.GT.DLT) GO TO 220
C IF(ETA*RNWTLN .LE. DLT)
C THEN
C
C TAKE PARTIAL STEP IN NEWTON DIRECTION
C
DO 210 I=1,N
SC(I)=(DLT/RNWTLN)*P(I)
210 CONTINUE
C$ WRITE(IPR,957)
GO TO 700
C ELSE
220 IF(CLN.LT.DLT) GO TO 240
C IF(CLN.GE.DLT)
C THEN
C TAKE STEP IN STEEPEST DESCENT DIRECTION
C
DO 230 I=1,N
SC(I)=(DLT/CLN)*SSD(I)/SX(I)
230 CONTINUE
C$ WRITE(IPR,958)
GO TO 700
C ELSE
C
C CALCULATE CONVEX COMBINATION OF SSD AND ETA*P
C WHICH HAS SCALED LENGTH DLT
C
240 DOT1=SDOT(N,V,1,SSD,1)
DOT2=SDOT(N,V,1,V,1)
ALAM=(-DOT1+SQRT((DOT1*DOT1)-DOT2*(CLN*CLN-DLT*DLT)))/DOT2
DO 250 I=1,N
SC(I)=(SSD(I) + ALAM*V(I))/SX(I)
250 CONTINUE
C$ WRITE(IPR,959)
C ENDIF
C ENDIF
C ENDIF
700 CONTINUE
C$ WRITE(IPR,952) FSTDOG,NWTAKE,RNWTLN,DLT
C$ WRITE(IPR,953)
C$ WRITE(IPR,960) (SC(I),I=1,N)
RETURN
C
951 FORMAT(27H0DOGSTP TAKE NEWTON STEP)
952 FORMAT(18H DOGSTP FSTDOG=,L1/
+ 18H DOGSTP NWTAKE=,L1/
+ 18H DOGSTP RNWTLN=,E20.13/
+ 18H DOGSTP DLT =,E20.13)
953 FORMAT(28H DOGSTP CURRENT STEP (SC))
954 FORMAT(18H DOGSTP ALPHA =,E20.13/
+ 18H DOGSTP BETA =,E20.13/
+ 18H DOGSTP CLN =,E20.13/
+ 18H DOGSTP ETA =,E20.13)
955 FORMAT(28H DOGSTP CAUCHY STEP (SSD))
956 FORMAT(12H DOGSTP V)
957 FORMAT(48H0DOGSTP TAKE PARTIAL STEP IN NEWTON DIRECTION)
958 FORMAT(50H0DOGSTP TAKE STEP IN STEEPEST DESCENT DIRECTION)
959 FORMAT(39H0DOGSTP TAKE CONVEX COMBINATION STEP)
960 FORMAT(14H DOGSTP ,5(E20.13,3X))
END
SUBROUTINE D1FCN(N,X,G)
C
C PURPOSE
C -------
C DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC
C WHEN SPECIFIC ANALYTIC GRADIENT FUNCTION NOT SUPPLIED.
C
DIMENSION X(N),G(N)
STOP
END
SUBROUTINE D2FCN(NR,N,X,H)
C
C PURPOSE
C -------
C DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC
C WHEN SPECIFIC ANALYTIC HESSIAN FUNCTION NOT SUPPLIED.
C
DIMENSION X(N),H(NR,1)
STOP
END
SUBROUTINE FORSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE AX=B WHERE A IS LOWER TRIANGULAR MATRIX
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) --> LOWER TRIANGULAR MATRIX (PRESERVED)
C X(N) <-- SOLUTION VECTOR
C B(N) --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE,
C THEN VECTORS B AND X MAY SHARE THE SAME STORAGE.
C
DIMENSION A(NR,1),X(N),B(N)
C
C SOLVE LX=B. (FOREWARD SOLVE)
C
X(1)=B(1)/A(1,1)
IF(N.EQ.1) RETURN
DO 20 I=2,N
SUM=0.0
IM1=I-1
DO 10 J=1,IM1
SUM=SUM+A(I,J)*X(J)
10 CONTINUE
X(I)=(B(I)-SUM)/A(I,I)
20 CONTINUE
RETURN
END
SUBROUTINE FSTOCD (N, X, FCN, SX, RNOISE, G)
C PURPOSE
C -------
C FIND CENTRAL DIFFERENCE APPROXIMATION G TO THE FIRST DERIVATIVE
C (GRADIENT) OF THE FUNCTION DEFINED BY FCN AT THE POINT X.
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C X --> POINT AT WHICH GRADIENT IS TO BE APPROXIMATED.
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION.
C SX --> DIAGONAL SCALING MATRIX FOR X.
C RNOISE --> RELATIVE NOISE IN FCN [F(X)].
C G <-- CENTRAL DIFFERENCE APPROXIMATION TO GRADIENT.
C
C
REAL X(N)
REAL SX(N)
REAL G(N)
EXTERNAL FCN
C
C FIND I TH STEPSIZE, EVALUATE TWO NEIGHBORS IN DIRECTION OF I TH
C UNIT VECTOR, AND EVALUATE I TH COMPONENT OF GRADIENT.
C
THIRD = 1.0/3.0
DO 10 I = 1, N
STEPI = RNOISE**THIRD * AMAX1(ABS(X(I)), 1.0/SX(I))
XTEMPI = X(I)
X(I) = XTEMPI + STEPI
CALL FCN (N, X, FPLUS)
X(I) = XTEMPI - STEPI
CALL FCN (N, X, FMINUS)
X(I) = XTEMPI
G(I) = (FPLUS - FMINUS)/(2.0*STEPI)
10 CONTINUE
RETURN
END
SUBROUTINE FSTOFD(NR,M,N,XPLS,FCN,FPLS,A,SX,RNOISE,FHAT,ICASE)
C PURPOSE
C -------
C FIND FIRST ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" TO THE
C FIRST DERIVATIVE OF THE FUNCTION DEFINED BY THE SUBPROGRAM "FNAME"
C EVALUATED AT THE NEW ITERATE "XPLS".
C
C
C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE:
C 1) THE FIRST DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN
C ANALYTIC USER ROUTINE HAS BEEN SUPPLIED;
C 2) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION
C IF NO ANALYTIC USER ROUTINE HAS BEEN SUPPLIED FOR THE HESSIAN BUT
C ONE HAS BEEN SUPPLIED FOR THE GRADIENT ("FCN") AND IF THE
C OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE
C
C NOTE
C ----
C _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE GRADIENT OF THE FUNCTION
C (FCN). FCN(X) # F: R(N)-->R(1)
C _M=N (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION
C FCN(X) # F: R(N)-->R(N).
C _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN OF THE OPTIMIZATIO
C FUNCTION, WHERE THE HESSIAN IS THE FIRST DERIVATIVE OF "FCN"
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C M --> NUMBER OF ROWS IN A
C N --> NUMBER OF COLUMNS IN A; DIMENSION OF PROBLEM
C XPLS(N) --> NEW ITERATE: X[K]
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C FPLS(M) --> _M=1 (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE:
C FCN(XPLS)
C _M=N (OPTIMIZATION) VALUE OF FIRST DERIVATIVE
C (GRADIENT) GIVEN BY USER FUNCTION FCN
C _M=N (SYSTEMS) FUNCTION VALUE OF ASSOCIATED
C MINIMIZATION FUNCTION
C A(NR,N) <-- FINITE DIFFERENCE APPROXIMATION (SEE NOTE). ONLY
C LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C RNOISE --> RELATIVE NOISE IN FCN [F(X)]
C FHAT(M) --> WORKSPACE
C ICASE --> =1 OPTIMIZATION (GRADIENT)
C =2 SYSTEMS
C =3 OPTIMIZATION (HESSIAN)
C
C INTERNAL VARIABLES
C ------------------
C STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION
C
DIMENSION XPLS(N),FPLS(M)
DIMENSION FHAT(M)
DIMENSION SX(N)
DIMENSION A(NR,1)
C
C FIND J-TH COLUMN OF A
C EACH COLUMN IS DERIVATIVE OF F(FCN) WITH RESPECT TO XPLS(J)
C
DO 30 J=1,N
STEPSZ=SQRT(RNOISE)*AMAX1(ABS(XPLS(J)),1./SX(J))
XTMPJ=XPLS(J)
XPLS(J)=XTMPJ+STEPSZ
CALL FCN(N,XPLS,FHAT)
XPLS(J)=XTMPJ
DO 20 I=1,M
A(I,J)=(FHAT(I)-FPLS(I))/STEPSZ
20 CONTINUE
30 CONTINUE
IF(ICASE.NE.3) RETURN
C
C IF COMPUTING HESSIAN, A MUST BE SYMMETRIC
C
IF(N.EQ.1) RETURN
NM1=N-1
DO 50 J=1,NM1
JP1=J+1
DO 40 I=JP1,M
A(I,J)=(A(I,J)+A(J,I))/2.0
40 CONTINUE
50 CONTINUE
RETURN
END
SUBROUTINE GRDCHK(N,X,FCN,F,G,TYPSIZ,SX,FSCALE,RNF,
+ ANALTL,WRK1,MSG,IPR)
C
C PURPOSE
C -------
C CHECK ANALYTIC GRADIENT AGAINST ESTIMATED GRADIENT
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C X(N) --> ESTIMATE TO A ROOT OF FCN
C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C FCN: R(N) --> R(1)
C F --> FUNCTION VALUE: FCN(X)
C G(N) --> GRADIENT: G(X)
C TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X
C SX(N) --> DIAGONAL SCALING MATRIX: SX(I)=1./TYPSIZ(I)
C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN
C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND
C ANALYTICAL GRADIENTS
C WRK1(N) --> WORKSPACE
C MSG <-- MESSAGE OR ERROR CODE
C ON OUTPUT: =-21, PROBABLE CODING ERROR OF GRADIENT
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
DIMENSION X(N),G(N)
DIMENSION SX(N),TYPSIZ(N)
DIMENSION WRK1(N)
EXTERNAL FCN
C
C COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO
C ANALYTIC GRADIENT.
C
CALL FSTOFD(1,1,N,X,FCN,F,WRK1,SX,RNF,WRK,1)
KER=0
DO 5 I=1,N
GS=AMAX1(ABS(F),FSCALE)/AMAX1(ABS(X(I)),TYPSIZ(I))
IF(ABS(G(I)-WRK1(I)).GT.AMAX1(ABS(G(I)),GS)*ANALTL) KER=1
5 CONTINUE
IF(KER.EQ.0) GO TO 20
WRITE(IPR,901)
WRITE(IPR,902) (I,G(I),WRK1(I),I=1,N)
MSG=-21
20 CONTINUE
RETURN
901 FORMAT(47H0GRDCHK PROBABLE ERROR IN CODING OF ANALYTIC,
+ 19H GRADIENT FUNCTION./
+ 16H GRDCHK COMP,12X,8HANALYTIC,12X,8HESTIMATE)
902 FORMAT(11H GRDCHK ,I5,3X,E20.13,3X,E20.13)
END
SUBROUTINE HESCHK(NR,N,X,FCN,D1FCN,D2FCN,F,G,A,TYPSIZ,SX,RNF,
+ ANALTL,IAGFLG,UDIAG,WRK1,WRK2,MSG,IPR)
C
C PURPOSE
C -------
C CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN
C (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN
C D2FCN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF A)
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> ESTIMATE TO A ROOT OF FCN
C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C FCN: R(N) --> R(1)
C D1FCN --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN.
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN.
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C F --> FUNCTION VALUE: FCN(X)
C G(N) <-- GRADIENT: G(X)
C A(N,N) <-- ON EXIT: HESSIAN IN LOWER TRIANGULAR PART AND DIAG
C TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X
C SX(N) --> DIAGONAL SCALING MATRIX: SX(I)=1./TYPSIZ(I)
C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND
C ANALYTICAL GRADIENTS
C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED
C UDIAG(N) --> WORKSPACE
C WRK1(N) --> WORKSPACE
C WRK2(N) --> WORKSPACE
C MSG <--> MESSAGE OR ERROR CODE
C ON INPUT : IF =1XX DO NOT COMPARE ANAL + EST HESS
C ON OUTPUT: =-22, PROBABLE CODING ERROR OF HESSIAN
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
DIMENSION X(N),G(N),A(NR,1)
DIMENSION TYPSIZ(N),SX(N)
DIMENSION UDIAG(N),WRK1(N),WRK2(N)
EXTERNAL FCN,D1FCN
C
C COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN.
C
IF(IAGFLG.EQ.1) CALL FSTOFD(NR,N,N,X,D1FCN,G,A,SX,RNF,WRK1,3)
IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,X,FCN,F,A,SX,RNF,WRK1,WRK2)
C
KER=0
C
C COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
C AND DIAGONAL OF "A" TO UDIAG
C
DO 30 J=1,N
UDIAG(J)=A(J,J)
IF(J.EQ.N) GO TO 30
JP1=J+1
DO 25 I=JP1,N
A(J,I)=A(I,J)
25 CONTINUE
30 CONTINUE
C
C COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE
C APPROXIMATION.
C
CALL D2FCN(NR,N,X,A)
DO 40 J=1,N
HS=AMAX1(ABS(G(J)),1.0)/AMAX1(ABS(X(J)),TYPSIZ(J))
IF(ABS(A(J,J)-UDIAG(J)).GT.AMAX1(ABS(UDIAG(J)),HS)*ANALTL)
+ KER=1
IF(J.EQ.N) GO TO 40
JP1=J+1
DO 35 I=JP1,N
IF(ABS(A(I,J)-A(J,I)).GT.AMAX1(ABS(A(I,J)),HS)*ANALTL) KER=1
35 CONTINUE
40 CONTINUE
C
IF(KER.EQ.0) GO TO 90
WRITE(IPR,901)
DO 50 I=1,N
IF(I.EQ.1) GO TO 45
IM1=I-1
DO 43 J=1,IM1
WRITE(IPR,902) I,J,A(I,J),A(J,I)
43 CONTINUE
45 WRITE(IPR,902) I,I,A(I,I),UDIAG(I)
50 CONTINUE
MSG=-22
C ENDIF
90 CONTINUE
RETURN
901 FORMAT(47H HESCHK PROBABLE ERROR IN CODING OF ANALYTIC,
+ 18H HESSIAN FUNCTION./
+ 21H HESCHK ROW COL,14X,8HANALYTIC,14X,10H(ESTIMATE))
902 FORMAT(11H HESCHK ,2I5,2X,E20.13,2X,1H(,E20.13,1H))
END
SUBROUTINE HOOKDR(NR,N,X,F,G,A,UDIAG,P,XPLS,FPLS,FCN,SX,STEPMX,
+ STEPTL,DLT,IRETCD,MXTAKE,AMU,DLTP,PHI,PHIP0,
+ SC,XPLSP,WRK0,EPSM,ITNCNT,IPR)
C
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE (XPLS) BY THE MORE-HEBDON METHOD
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE X[K-1]
C F --> FUNCTION VALUE AT OLD ITERATE, F(X)
C G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE
C A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN LOWER
C TRIANGULAR PART AND DIAGONAL.
C HESSIAN IN UPPER TRIANGULAR PART AND UDIAG.
C UDIAG(N) --> DIAGONAL OF HESSIAN IN A(.,.)
C P(N) --> NEWTON STEP
C XPLS(N) <-- NEW ITERATE X[K]
C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C DLT <--> TRUST REGION RADIUS
C IRETCD <-- RETURN CODE
C =0 SATISFACTORY XPLS FOUND
C =1 FAILED TO FIND SATISFACTORY XPLS SUFFICIENTLY
C DISTINCT FROM X
C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C AMU <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C DLTP <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHI <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHIP0 <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C SC(N) --> WORKSPACE
C XPLSP(N) --> WORKSPACE
C WRK0(N) --> WORKSPACE
C EPSM --> MACHINE EPSILON
C ITNCNT --> ITERATION COUNT
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
REAL X(N),G(N),P(N),XPLS(N),SX(N)
REAL A(NR,1),UDIAG(N)
REAL SC(N),XPLSP(N),WRK0(N)
LOGICAL MXTAKE,NWTAKE
LOGICAL FSTIME
EXTERNAL FCN
C
IRETCD=4
FSTIME=.TRUE.
TMP=0.
DO 5 I=1,N
TMP=TMP+SX(I)*SX(I)*P(I)*P(I)
5 CONTINUE
RNWTLN=SQRT(TMP)
C$ WRITE(IPR,954) RNWTLN
C
IF(ITNCNT.GT.1) GO TO 100
C IF(ITNCNT.EQ.1)
C THEN
AMU=0.
C
C IF FIRST ITERATION AND TRUST REGION NOT PROVIDED BY USER,
C COMPUTE INITIAL TRUST REGION.
C
IF(DLT.NE. (-1.)) GO TO 100
C IF(DLT.EQ. (-1.))
C THEN
ALPHA=0.
DO 10 I=1,N
ALPHA=ALPHA+(G(I)*G(I))/(SX(I)*SX(I))
10 CONTINUE
BETA=0.0
DO 30 I=1,N
TMP=0.
DO 20 J=I,N
TMP=TMP + (A(J,I)*G(J))/(SX(J)*SX(J))
20 CONTINUE
BETA=BETA+TMP*TMP
30 CONTINUE
DLT=ALPHA*SQRT(ALPHA)/BETA
DLT = AMIN1(DLT, STEPMX)
C$ WRITE(IPR,950)
C$ WRITE(IPR,951) ALPHA,BETA,DLT
C ENDIF
C ENDIF
C
100 CONTINUE
C
C FIND NEW STEP BY MORE-HEBDON ALGORITHM
CALL HOOKST(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU,
+ DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
DLTP=DLT
C
C CHECK NEW POINT AND UPDATE TRUST REGION
CALL TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,
+ DLT,IRETCD,XPLSP,FPLSP,XPLS,FPLS,MXTAKE,IPR,3,UDIAG)
IF(IRETCD.LE.1) RETURN
GO TO 100
C
950 FORMAT(43H HOOKDR INITIAL TRUST REGION NOT GIVEN. ,
+ 21H COMPUTE CAUCHY STEP.)
951 FORMAT(18H HOOKDR ALPHA =,E20.13/
+ 18H HOOKDR BETA =,E20.13/
+ 18H HOOKDR DLT =,E20.13)
952 FORMAT(28H HOOKDR CURRENT STEP (SC))
954 FORMAT(18H0HOOKDR RNWTLN=,E20.13)
955 FORMAT(14H HOOKDR ,5(E20.13,3X))
END
SUBROUTINE HOOKST(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU,
+ DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
C
C PURPOSE
C -------
C FIND NEW STEP BY MORE-HEBDON ALGORITHM
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C G(N) --> GRADIENT AT CURRENT ITERATE, G(X)
C A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN
C LOWER TRIANGULAR PART AND DIAGONAL.
C HESSIAN OR APPROX IN UPPER TRIANGULAR PART
C UDIAG(N) --> DIAGONAL OF HESSIAN IN A(.,.)
C P(N) --> NEWTON STEP
C SX(N) --> DIAGONAL SCALING MATRIX FOR N
C RNWTLN --> NEWTON STEP LENGTH
C DLT <--> TRUST REGION RADIUS
C AMU <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C DLTP --> TRUST REGION RADIUS AT LAST EXIT FROM THIS ROUTINE
C PHI <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHIP0 <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C FSTIME <--> BOOLEAN. =.TRUE. IF FIRST ENTRY TO THIS ROUTINE
C DURING K-TH ITERATION
C SC(N) <-- CURRENT STEP
C NWTAKE <-- BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C WRK0 --> WORKSPACE
C EPSM --> MACHINE EPSILON
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
REAL G(N),P(N),SX(N),SC(N),WRK0(N)
REAL A(NR,1),UDIAG(N)
LOGICAL NWTAKE,DONE
LOGICAL FSTIME
C
C HI AND ALO ARE CONSTANTS USED IN THIS ROUTINE.
C CHANGE HERE IF OTHER VALUES ARE TO BE SUBSTITUTED.
HI=1.5
ALO=.75
C -----
IF(RNWTLN.GT.HI*DLT) GO TO 15
C IF(RNWTLN.LE.HI*DLT)
C THEN
C
C TAKE NEWTON STEP
C
NWTAKE=.TRUE.
DO 10 I=1,N
SC(I)=P(I)
10 CONTINUE
DLT=AMIN1(DLT,RNWTLN)
AMU=0.
C$ WRITE(IPR,951)
RETURN
C ELSE
C
C NEWTON STEP NOT TAKEN
C
15 CONTINUE
C$ WRITE(IPR,952)
NWTAKE=.FALSE.
IF(AMU.LE.0.) GO TO 20
C IF(AMU.GT.0.)
C THEN
AMU=AMU- (PHI+DLTP) *((DLTP-DLT)+PHI)/(DLT*PHIP)
C$ WRITE(IPR,956) AMU
C ENDIF
20 CONTINUE
PHI=RNWTLN-DLT
IF(.NOT.FSTIME) GO TO 28
C IF(FSTIME)
C THEN
DO 25 I=1,N
WRK0(I)=SX(I)*SX(I)*P(I)
25 CONTINUE
C
C SOLVE L*Y = (SX**2)*P
C
CALL FORSLV(NR,N,A,WRK0,WRK0)
PHIP0=-SNRM2(N,WRK0,1)**2/RNWTLN
FSTIME=.FALSE.
C ENDIF
28 PHIP=PHIP0
AMULO=-PHI/PHIP
AMUUP=0.0
DO 30 I=1,N
AMUUP=AMUUP+(G(I)*G(I))/(SX(I)*SX(I))
30 CONTINUE
AMUUP=SQRT(AMUUP)/DLT
DONE=.FALSE.
C$ WRITE(IPR,956) AMU
C$ WRITE(IPR,959) PHI
C$ WRITE(IPR,960) PHIP
C$ WRITE(IPR,957) AMULO
C$ WRITE(IPR,958) AMUUP
C
C TEST VALUE OF AMU; GENERATE NEXT AMU IF NECESSARY
C
100 CONTINUE
IF(DONE) RETURN
C$ WRITE(IPR,962)
IF(AMU.GE.AMULO .AND. AMU.LE.AMUUP) GO TO 110
C IF(AMU.LT.AMULO .OR. AMU.GT.AMUUP)
C THEN
AMU=AMAX1(SQRT(AMULO*AMUUP),AMUUP*1.0E-3)
C$ WRITE(IPR,956) AMU
C ENDIF
110 CONTINUE
C
C COPY (H,UDIAG) TO L
C WHERE H <-- H+AMU*(SX**2) [DO NOT ACTUALLY CHANGE (H,UDIAG)]
DO 130 J=1,N
A(J,J)=UDIAG(J) + AMU*SX(J)*SX(J)
IF(J.EQ.N) GO TO 130
JP1=J+1
DO 120 I=JP1,N
A(I,J)=A(J,I)
120 CONTINUE
130 CONTINUE
C
C FACTOR H=L(L+)
C
CALL CHOLDC(NR,N,A,0.,SQRT(EPSM),ADDMAX)
C
C SOLVE H*P = L(L+)*SC = -G
C
DO 140 I=1,N
WRK0(I)=-G(I)
140 CONTINUE
CALL LLTSLV(NR,N,A,SC,WRK0)
C$ WRITE(IPR,955)
C$ WRITE(IPR,963) (SC(I),I=1,N)
C
C RESET H. NOTE SINCE UDIAG HAS NOT BEEN DESTROYED WE NEED DO
C NOTHING HERE. H IS IN THE UPPER PART AND IN UDIAG, STILL INTACT
C
STEPLN=0.
DO 150 I=1,N
STEPLN=STEPLN + SX(I)*SX(I)*SC(I)*SC(I)
150 CONTINUE
STEPLN=SQRT(STEPLN)
PHI=STEPLN-DLT
DO 160 I=1,N
WRK0(I)=SX(I)*SX(I)*SC(I)
160 CONTINUE
CALL FORSLV(NR,N,A,WRK0,WRK0)
PHIP=-SNRM2(N,WRK0,1)**2/STEPLN
C$ WRITE(IPR,961) DLT,STEPLN
C$ WRITE(IPR,959) PHI
C$ WRITE(IPR,960) PHIP
IF((ALO*DLT.GT.STEPLN .OR. STEPLN.GT.HI*DLT) .AND.
+ (AMUUP-AMULO.GT.0.)) GO TO 170
C IF((ALO*DLT.LE.STEPLN .AND. STEPLN.LE.HI*DLT) .OR.
C (AMUUP-AMULO.LE.0.))
C THEN
C
C SC IS ACCEPTABLE HOOKSTEP
C
C$ WRITE(IPR,954)
DONE=.TRUE.
GO TO 100
C ELSE
C
C SC NOT ACCEPTABLE HOOKSTEP. SELECT NEW AMU
C
170 CONTINUE
C$ WRITE(IPR,953)
AMULO=AMAX1(AMULO,AMU-(PHI/PHIP))
IF(PHI.LT.0.) AMUUP=AMIN1(AMUUP,AMU)
AMU=AMU-(STEPLN*PHI)/(DLT*PHIP)
C$ WRITE(IPR,956) AMU
C$ WRITE(IPR,957) AMULO
C$ WRITE(IPR,958) AMUUP
GO TO 100
C ENDIF
C ENDIF
C
951 FORMAT(27H0HOOKST TAKE NEWTON STEP)
952 FORMAT(32H0HOOKST NEWTON STEP NOT TAKEN)
953 FORMAT(31H HOOKST SC IS NOT ACCEPTABLE)
954 FORMAT(27H HOOKST SC IS ACCEPTABLE)
955 FORMAT(28H HOOKST CURRENT STEP (SC))
956 FORMAT(18H HOOKST AMU =,E20.13)
957 FORMAT(18H HOOKST AMULO =,E20.13)
958 FORMAT(18H HOOKST AMUUP =,E20.13)
959 FORMAT(18H HOOKST PHI =,E20.13)
960 FORMAT(18H HOOKST PHIP =,E20.13)
961 FORMAT(18H HOOKST DLT =,E20.13/
+ 18H HOOKST STEPLN=,E20.13)
962 FORMAT(23H0HOOKST FIND NEW AMU)
963 FORMAT(14H HOOKST ,5(E20.13,3X))
END
SUBROUTINE HSNINT(NR,N,A,SX,METHOD)
C
C PURPOSE
C -------
C PROVIDE INITIAL HESSIAN WHEN USING SECANT UPDATES
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) <-- INITIAL HESSIAN (LOWER TRIANGULAR MATRIX)
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C =1,2 FACTORED SECANT METHOD USED
C =3 UNFACTORED SECANT METHOD USED
C
DIMENSION A(NR,1),SX(N)
C
DO 100 J=1,N
IF(METHOD.EQ.3) A(J,J)=SX(J)*SX(J)
IF(METHOD.NE.3) A(J,J)=SX(J)
IF(J.EQ.N) GO TO 100
JP1=J+1
DO 90 I=JP1,N
A(I,J)=0.
90 CONTINUE
100 CONTINUE
RETURN
END
SUBROUTINE LLTSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE AX=B WHERE A HAS THE FORM L(L-TRANSPOSE)
C BUT ONLY THE LOWER TRIANGULAR PART, L, IS STORED.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) --> MATRIX OF FORM L(L-TRANSPOSE).
C ON RETURN A IS UNCHANGED.
C X(N) <-- SOLUTION VECTOR
C B(N) --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NOT REQUIRED BY CALLING PROGRAM, THEN
C B AND X MAY SHARE THE SAME STORAGE.
C
DIMENSION A(NR,1),X(N),B(N)
C
C FORWARD SOLVE, RESULT IN X
C
CALL FORSLV(NR,N,A,X,B)
C
C BACK SOLVE, RESULT IN X
C
CALL BAKSLV(NR,N,A,X,X)
RETURN
END
SUBROUTINE LNSRCH(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE,
+ IRETCD,STEPMX,STEPTL,SX,IPR)
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE BY LINE SEARCH.
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE: X[K-1]
C F --> FUNCTION VALUE AT OLD ITERATE, F(X)
C G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE
C P(N) --> NON-ZERO NEWTON STEP
C XPLS(N) <-- NEW ITERATE X[K]
C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C IRETCD <-- RETURN CODE
C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
C INTERNAL VARIABLES
C ------------------
C SLN NEWTON LENGTH
C RLN RELATIVE LENGTH OF NEWTON STEP
C
INTEGER N,IRETCD
REAL STEPMX,STEPTL
REAL SX(N)
REAL X(N),G(N),P(N)
REAL XPLS(N)
REAL F,FPLS
LOGICAL MXTAKE
C
MXTAKE=.FALSE.
IRETCD=2
C$ WRITE(IPR,954)
C$ WRITE(IPR,955) (P(I),I=1,N)
TMP=0.0
DO 5 I=1,N
TMP=TMP+SX(I)*SX(I)*P(I)*P(I)
5 CONTINUE
SLN=SQRT(TMP)
IF(SLN.LE.STEPMX) GO TO 10
C
C NEWTON STEP LONGER THAN MAXIMUM ALLOWED
SCL=STEPMX/SLN
CALL SCLMUL(N,SCL,P,P)
SLN=STEPMX
C$ WRITE(IPR,954)
C$ WRITE(IPR,955) (P(I),I=1,N)
10 CONTINUE
SLP=SDOT(N,G,1,P,1)
RLN=0.
DO 15 I=1,N
RLN=AMAX1(RLN,ABS(P(I))/AMAX1(ABS(X(I)),1./SX(I)))
15 CONTINUE
RMNLMB=STEPTL/RLN
ALMBDA=1.0
C$ WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL
C
C LOOP
C CHECK IF NEW ITERATE SATISFACTORY. GENERATE NEW LAMBDA IF NECESSARY.
C
100 CONTINUE
IF(IRETCD.LT.2) RETURN
DO 105 I=1,N
XPLS(I)=X(I) + ALMBDA*P(I)
105 CONTINUE
CALL FCN(N,XPLS,FPLS)
C$ WRITE(IPR,950) ALMBDA
C$ WRITE(IPR,951)
C$ WRITE(IPR,955) (XPLS(I),I=1,N)
C$ WRITE(IPR,953) FPLS
IF(FPLS.GT. F+SLP*1.E-4*ALMBDA) GO TO 130
C IF(FPLS.LE. F+SLP*1.E-4*ALMBDA)
C THEN
C
C SOLUTION FOUND
C
IRETCD=0
IF(ALMBDA.EQ.1.0 .AND. SLN.GT. .99*STEPMX) MXTAKE=.TRUE.
GO TO 100
C
C SOLUTION NOT (YET) FOUND
C
C ELSE
130 IF(ALMBDA .GE. RMNLMB) GO TO 140
C IF(ALMBDA .LT. RMNLMB)
C THEN
C
C NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X
C
IRETCD=1
GO TO 100
C ELSE
C
C CALCULATE NEW LAMBDA
C
140 IF(ALMBDA.NE.1.0) GO TO 150
C IF(ALMBDA.EQ.1.0)
C THEN
C
C FIRST BACKTRACK: QUADRATIC FIT
C
TLMBDA=-SLP/(2.*(FPLS-F-SLP))
GO TO 170
C ELSE
C
C ALL SUBSEQUENT BACKTRACKS: CUBIC FIT
C
150 T1=FPLS-F-ALMBDA*SLP
T2=PFPLS-F-PLMBDA*SLP
T3=1.0/(ALMBDA-PLMBDA)
A=T3*(T1/(ALMBDA*ALMBDA) - T2/(PLMBDA*PLMBDA))
B=T3*(T2*ALMBDA/(PLMBDA*PLMBDA)
+ - T1*PLMBDA/(ALMBDA*ALMBDA) )
DISC=B*B-3.0*A*SLP
IF(DISC.LE. B*B) GO TO 160
C IF(DISC.GT. B*B)
C THEN
C
C ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM
C
TLMBDA=(-B+SIGN(1.0,A)*SQRT(DISC))/(3.0*A)
GO TO 165
C ELSE
C
C BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM
C
160 TLMBDA=(-B-SIGN(1.0,A)*SQRT(DISC))/(3.0*A)
C ENDIF
165 IF(TLMBDA.GT. .5*ALMBDA) TLMBDA=.5*ALMBDA
C ENDIF
170 PLMBDA=ALMBDA
PFPLS=FPLS
IF(TLMBDA.GE. ALMBDA*.1) GO TO 180
C IF(TLMBDA.LT.ALMBDA/10.)
C THEN
ALMBDA=ALMBDA*.1
GO TO 190
C ELSE
180 ALMBDA=TLMBDA
C ENDIF
C ENDIF
C ENDIF
190 GO TO 100
950 FORMAT(18H LNSRCH ALMBDA=,E20.13)
951 FORMAT(29H LNSRCH NEW ITERATE (XPLS))
952 FORMAT(18H LNSRCH SLN =,E20.13/
+ 18H LNSRCH SLP =,E20.13/
+ 18H LNSRCH RMNLMB=,E20.13/
+ 18H LNSRCH STEPMX=,E20.13/
+ 18H LNSRCH STEPTL=,E20.13)
953 FORMAT(19H LNSRCH F(XPLS)=,E20.13)
954 FORMAT(26H0LNSRCH NEWTON STEP (P))
955 FORMAT(14H LNSRCH ,5(E20.13,3X))
END
SUBROUTINE OPTCHK(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM,
+ DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR)
C
C PURPOSE
C -------
C CHECK INPUT FOR REASONABLENESS
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C X(N) --> ON ENTRY, ESTIMATE TO ROOT OF FCN
C TYPSIZ(N) <--> TYPICAL SIZE OF EACH COMPONENT OF X
C SX(N) <-- DIAGONAL SCALING MATRIX FOR X
C FSCALE <--> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN
C GRADTL --> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE
C ENOUGH TO ZERO TO TERMINATE ALGORITHM
C ITNLIM <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C NDIGIT <--> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN
C EPSM --> MACHINE EPSILON
C DLT <--> TRUST REGION RADIUS
C METHOD <--> ALGORITHM INDICATOR
C IEXP <--> EXPENSE FLAG
C IAGFLG <--> =1 IF ANALYTIC GRADIENT SUPPLIED
C IAHFLG <--> =1 IF ANALYTIC HESSIAN SUPPLIED
C STEPMX <--> MAXIMUM STEP SIZE
C MSG <--> MESSAGE AND ERROR CODE
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
DIMENSION X(N),TYPSIZ(N),SX(N)
C
C CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES.
C IF NOT, SET THEM TO DEFAULT VALUES.
IF(METHOD.LT.1 .OR. METHOD.GT.3) METHOD=1
IF(IAGFLG.NE.1) IAGFLG=0
IF(IAHFLG.NE.1) IAHFLG=0
IF(IEXP.NE.0) IEXP=1
IF(MOD(MSG/2,2).EQ.1 .AND. IAGFLG.EQ.0) GO TO 830
IF(MOD(MSG/4,2).EQ.1 .AND. IAHFLG.EQ.0) GO TO 835
C
C CHECK DIMENSION OF PROBLEM
C
IF(N.LE.0) GO TO 805
IF(N.EQ.1 .AND. MOD(MSG,2).EQ.0) GO TO 810
C
C COMPUTE SCALE MATRIX
C
DO 10 I=1,N
IF(TYPSIZ(I).EQ.0.) TYPSIZ(I)=1.0
IF(TYPSIZ(I).LT.0.) TYPSIZ(I)=-TYPSIZ(I)
SX(I)=1.0/TYPSIZ(I)
10 CONTINUE
C
C CHECK MAXIMUM STEP SIZE
C
IF (STEPMX .GT. 0.0) GO TO 20
STPSIZ = 0.0
DO 15 I = 1, N
STPSIZ = STPSIZ + X(I)*X(I)*SX(I)*SX(I)
15 CONTINUE
STPSIZ = SQRT(STPSIZ)
STEPMX = AMAX1(1.0E3*STPSIZ, 1.0E3)
20 CONTINUE
C CHECK FUNCTION SCALE
IF(FSCALE.EQ.0.) FSCALE=1.0
IF(FSCALE.LT.0.) FSCALE=-FSCALE
C
C CHECK GRADIENT TOLERANCE
IF(GRADTL.LT.0.) GO TO 815
C
C CHECK ITERATION LIMIT
IF(ITNLIM.LE.0) GO TO 820
C
C CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN
IF(NDIGIT.EQ.0) GO TO 825
IF(NDIGIT.LT.0) NDIGIT=-ALOG10(EPSM)
C
C CHECK TRUST REGION RADIUS
IF(DLT.LE.0.) DLT=-1.0
IF (DLT .GT. STEPMX) DLT = STEPMX
RETURN
C
C ERROR EXITS
C
805 WRITE(IPR,901) N
MSG=-1
GO TO 895
810 WRITE(IPR,902)
MSG=-2
GO TO 895
815 WRITE(IPR,903) GRADTL
MSG=-3
GO TO 895
820 WRITE(IPR,904) ITNLIM
MSG=-4
GO TO 895
825 WRITE(IPR,905) NDIGIT
MSG=-5
GO TO 895
830 WRITE(IPR,906) MSG,IAGFLG
MSG=-6
GO TO 895
835 WRITE(IPR,907) MSG,IAHFLG
MSG=-7
895 RETURN
901 FORMAT(32H0OPTCHK ILLEGAL DIMENSION, N=,I5)
902 FORMAT(55H0OPTCHK +++ WARNING +++ THIS PACKAGE IS INEFFICIENT,
+ 26H FOR PROBLEMS OF SIZE N=1./
+ 48H OPTCHK CHECK INSTALLATION LIBRARIES FOR MORE,
+ 22H APPROPRIATE ROUTINES./
+ 41H OPTCHK IF NONE, SET MSG AND RESUBMIT.)
903 FORMAT(38H0OPTCHK ILLEGAL TOLERANCE. GRADTL=,E20.13)
904 FORMAT(44H0OPTCHK ILLEGAL ITERATION LIMIT. ITNLIM=,I5)
905 FORMAT(52H0OPTCHK MINIMIZATION FUNCTION HAS NO GOOD DIGITS.,
+ 9H NDIGIT=,I5)
906 FORMAT(50H0OPTCHK USER REQUESTS THAT ANALYTIC GRADIENT BE,
+ 33H ACCEPTED AS PROPERLY CODED (MSG=,I5, 2H),/
+ 45H OPTCHK BUT ANALYTIC GRADIENT NOT SUPPLIED,
+ 9H (IAGFLG=,I5, 2H).)
907 FORMAT(49H0OPTCHK USER REQUESTS THAT ANALYTIC HESSIAN BE,
+ 33H ACCEPTED AS PROPERLY CODED (MSG=,I5, 2H),/
+ 44H OPTCHK BUT ANALYTIC HESSIAN NOT SUPPLIED,
+ 9H (IAHFLG=,I5, 2H).)
END
SUBROUTINE OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE,
+ METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
+ DLT,GRADTL,STEPMX,STEPTL,
+ XPLS,FPLS,GPLS,ITRMCD,
+ A,UDIAG,G,P,SX,WRK0,WRK1,WRK2,WRK3)
C
C PURPOSE
C -------
C DRIVER FOR NON-LINEAR OPTIMIZATION PROBLEM
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> ON ENTRY: ESTIMATE TO A ROOT OF FCN
C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C FCN: R(N) --> R(1)
C D1FCN --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE GRADIENT
C OF FCN. MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE HESSIAN OF
C OF FCN. MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X
C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION
C METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C =1 LINE SEARCH
C =2 DOUBLE DOGLEG
C =3 MORE-HEBDON
C IEXP --> =1 IF OPTIMIZATION FUNCTION FCN IS EXPENSIVE TO
C EVALUATE, =0 OTHERWISE. IF SET THEN HESSIAN WILL
C BE EVALUATED BY SECANT UPDATE INSTEAD OF
C ANALYTICALLY OR BY FINITE DIFFERENCES
C MSG <--> ON INPUT: (.GT.0) MESSAGE TO INHIBIT CERTAIN
C AUTOMATIC CHECKS
C ON OUTPUT: (.LT.0) ERROR CODE; =0 NO ERROR
C NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN
C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED
C IAHFLG --> =1 IF ANALYTIC HESSIAN SUPPLIED
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C DLT --> TRUST REGION RADIUS
C GRADTL --> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE
C ENOUGH TO ZERO TO TERMINATE ALGORITHM
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C XPLS(N) <--> ON EXIT: XPLS IS LOCAL MINIMUM
C FPLS <--> ON EXIT: FUNCTION VALUE AT SOLUTION, XPLS
C GPLS(N) <--> ON EXIT: GRADIENT AT SOLUTION XPLS
C ITRMCD <-- TERMINATION CODE
C A(N,N) --> WORKSPACE FOR HESSIAN (OR ESTIMATE)
C AND ITS CHOLESKY DECOMPOSITION
C UDIAG(N) --> WORKSPACE [FOR DIAGONAL OF HESSIAN]
C G(N) --> WORKSPACE (FOR GRADIENT AT CURRENT ITERATE)
C P(N) --> WORKSPACE FOR STEP
C SX(N) --> WORKSPACE (FOR DIAGONAL SCALING MATRIX)
C WRK0(N) --> WORKSPACE
C WRK1(N) --> WORKSPACE
C WRK2(N) --> WORKSPACE
C WRK3(N) --> WORKSPACE
C
C
C INTERNAL VARIABLES
C ------------------
C ANALTL TOLERANCE FOR COMPARISON OF ESTIMATED AND
C ANALYTICAL GRADIENTS AND HESSIANS
C EPSM MACHINE EPSILON
C F FUNCTION VALUE: FCN(X)
C ITNCNT CURRENT ITERATION, K
C RNF RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN.
C NOISE=10.**(-NDIGIT)
C
DIMENSION X(N),XPLS(N),G(N),GPLS(N),P(N)
DIMENSION TYPSIZ(N),SX(N)
DIMENSION A(NR,1),UDIAG(N)
DIMENSION WRK0(N),WRK1(N),WRK2(N),WRK3(N)
LOGICAL MXTAKE,NOUPDT
EXTERNAL FCN,D1FCN,D2FCN
C
C INITIALIZATION
C --------------
DO 10 I=1,N
P(I)=0.
10 CONTINUE
ITNCNT=0
IRETCD=-1
EPSM=EPSMCH(1.0)
CALL OPTCHK(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM,
+ DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR)
IF(MSG.LT.0) RETURN
RNF=AMAX1(10.**(-NDIGIT),EPSM)
ANALTL=AMAX1(1.E-2,SQRT(RNF))
C
IF(MOD(MSG/8,2).EQ.1) GO TO 15
WRITE(IPR,901)
WRITE(IPR,900) (TYPSIZ(I),I=1,N)
WRITE(IPR,902)
WRITE(IPR,900) (SX(I),I=1,N)
WRITE(IPR,903) FSCALE
WRITE(IPR,904) NDIGIT,IAGFLG,IAHFLG,IEXP,METHOD,ITNLIM,EPSM
WRITE(IPR,905) STEPMX,STEPTL,GRADTL,DLT,RNF,ANALTL
15 CONTINUE
C
C EVALUATE FCN(X)
C
CALL FCN(N,X,F)
C
C EVALUATE ANALYTIC OR FINITE DIFFERENCE GRADIENT AND CHECK ANALYTIC
C GRADIENT, IF REQUESTED.
C
IF (IAGFLG .EQ. 1) GO TO 20
C IF (IAGFLG .EQ. 0)
C THEN
CALL FSTOFD (1, 1, N, X, FCN, F, G, SX, RNF, WRK, 1)
GO TO 25
C
20 CALL D1FCN (N, X, G)
IF (MOD(MSG/2,2) .EQ. 1) GO TO 25
C IF (MOD(MSG/2,2).EQ.0)
C THEN
CALL GRDCHK (N, X, FCN, F, G, TYPSIZ, SX, FSCALE,
1 RNF, ANALTL, WRK1, MSG, IPR)
IF (MSG .LT. 0) RETURN
25 CONTINUE
C
CALL OPTSTP(N,X,F,G,WRK1,ITNCNT,ICSCMX,
+ ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,
+ IPR,MSG)
IF(ITRMCD.NE.0) GO TO 700
C
IF(IEXP.NE.1) GO TO 80
C
C IF OPTIMIZATION FUNCTION EXPENSIVE TO EVALUATE (IEXP=1), THEN
C HESSIAN WILL BE OBTAINED BY SECANT UPDATES. GET INITIAL HESSIAN.
C
CALL HSNINT(NR,N,A,SX,METHOD)
GO TO 90
80 CONTINUE
C
C EVALUATE ANALYTIC OR FINITE DIFFERENCE HESSIAN AND CHECK ANALYTIC
C HESSIAN IF REQUESTED (ONLY IF USER-SUPPLIED ANALYTIC HESSIAN
C ROUTINE D2FCN FILLS ONLY LOWER TRIANGULAR PART AND DIAGONAL OF A).
C
IF (IAHFLG .EQ. 1) GO TO 82
C IF (IAHFLG .EQ. 0)
C THEN
IF (IAGFLG .EQ. 1) CALL FSTOFD (NR, N, N, X, D1FCN, G, A, SX,
1 RNF, WRK1, 3)
IF (IAGFLG .NE. 1) CALL SNDOFD (NR, N, X, FCN, F, A, SX, RNF,
1 WRK1, WRK2)
GO TO 88
C
C ELSE
82 IF (MOD(MSG/4,2).EQ.0) GO TO 85
C IF (MOD(MSG/4, 2) .EQ. 1)
C THEN
CALL D2FCN (NR, N, X, A)
GO TO 88
C
C ELSE
85 CALL HESCHK (NR, N, X, FCN, D1FCN, D2FCN, F, G, A, TYPSIZ,
1 SX, RNF, ANALTL, IAGFLG, UDIAG, WRK1, WRK2, MSG, IPR)
C
C HESCHK EVALUATES D2FCN AND CHECKS IT AGAINST THE FINITE
C DIFFERENCE HESSIAN WHICH IT CALCULATES BY CALLING FSTOFD
C (IF IAGFLG .EQ. 1) OR SNDOFD (OTHERWISE).
C
IF (MSG .LT. 0) RETURN
88 CONTINUE
C
90 IF(MOD(MSG/8,2).EQ.0)
+ CALL RESULT(NR,N,X,F,G,A,P,ITNCNT,1,IPR)
C
C
C ITERATION
C ---------
100 ITNCNT=ITNCNT+1
C
C FIND PERTURBED LOCAL MODEL HESSIAN AND ITS LL+ DECOMPOSITION
C (SKIP THIS STEP IF LINE SEARCH OR DOGSTEP TECHNIQUES BEING USED WITH
C SECANT UPDATES. CHOLESKY DECOMPOSITION L ALREADY OBTAINED FROM
C SECFAC.)
C
IF(IEXP.EQ.1 .AND. METHOD.NE.3) GO TO 105
103 CALL CHLHSN(NR,N,A,EPSM,SX,UDIAG)
105 CONTINUE
C
C SOLVE FOR NEWTON STEP: AP=-G
C
DO 110 I=1,N
WRK1(I)=-G(I)
110 CONTINUE
CALL LLTSLV(NR,N,A,P,WRK1)
C
C DECIDE WHETHER TO ACCEPT NEWTON STEP XPLS=X + P
C OR TO CHOOSE XPLS BY A GLOBAL STRATEGY.
C
IF (IAGFLG .NE. 0 .OR. METHOD .EQ. 1) GO TO 111
DLTSAV = DLT
IF (METHOD .EQ. 2) GO TO 111
AMUSAV = AMU
DLPSAV = DLTP
PHISAV = PHI
PHPSAV = PHIP0
111 IF(METHOD.EQ.1)
+ CALL LNSRCH(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE,IRETCD,
+ STEPMX,STEPTL,SX,IPR)
IF(METHOD.EQ.2)
+ CALL DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX,
+ STEPTL,DLT,IRETCD,MXTAKE,WRK0,WRK1,WRK2,WRK3,IPR)
IF(METHOD.EQ.3)
+ CALL HOOKDR(NR,N,X,F,G,A,UDIAG,P,XPLS,FPLS,FCN,SX,STEPMX,
+ STEPTL,DLT,IRETCD,MXTAKE,AMU,DLTP,PHI,PHIP0,WRK0,
+ WRK1,WRK2,EPSM,ITNCNT,IPR)
C
C IF COULD NOT FIND SATISFACTORY STEP AND FORWARD DIFFERENCE
C GRADIENT WAS USED, RETRY USING CENTRAL DIFFERENCE GRADIENT.
C
IF (IRETCD .NE. 1 .OR. IAGFLG .NE. 0) GO TO 112
C IF (IRETCD .EQ. 1 .AND. IAGFLG .EQ. 0)
C THEN
C
C SET IAGFLG FOR CENTRAL DIFFERENCES
C
IAGFLG = -1
WRITE(IPR,906) ITNCNT
C
CALL FSTOCD (N, X, FCN, SX, RNF, G)
IF (METHOD .EQ. 1) GO TO 105
DLT = DLTSAV
IF (METHOD .EQ. 2) GO TO 105
AMU = AMUSAV
DLTP = DLPSAV
PHI = PHISAV
PHIP0 = PHPSAV
GO TO 103
C ENDIF
C
C CALCULATE STEP FOR OUTPUT
C
112 CONTINUE
DO 114 I = 1, N
P(I) = XPLS(I) - X(I)
114 CONTINUE
C
C CALCULATE GRADIENT AT XPLS
C
IF (IAGFLG .EQ. (-1)) GO TO 116
IF (IAGFLG .EQ. 0) GO TO 118
C
C ANALYTIC GRADIENT
CALL D1FCN (N, XPLS, GPLS)
GO TO 120
C
C CENTRAL DIFFERENCE GRADIENT
116 CALL FSTOCD (N, XPLS, FCN, SX, RNF, GPLS)
GO TO 120
C
C FORWARD DIFFERENCE GRADIENT
118 CALL FSTOFD (1, 1, N, XPLS, FCN, FPLS, GPLS, SX, RNF, WRK, 1)
120 CONTINUE
C
C CHECK WHETHER STOPPING CRITERIA SATISFIED
C
CALL OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX,
+ ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,
+ IPR,MSG)
IF(ITRMCD.NE.0) GO TO 690
C
C EVALUATE HESSIAN AT XPLS
C
IF(IEXP.EQ.0) GO TO 130
IF(METHOD.EQ.3)
+ CALL SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT,RNF,
+ IAGFLG,NOUPDT,WRK1,WRK2,WRK3)
IF(METHOD.NE.3)
+ CALL SECFAC(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF,IAGFLG,
+ NOUPDT,WRK0,WRK1,WRK2,WRK3)
GO TO 150
130 IF(IAHFLG.EQ.1) GO TO 140
IF(IAGFLG.EQ.1)
+ CALL FSTOFD(NR,N,N,XPLS,D1FCN,GPLS,A,SX,RNF,WRK1,3)
IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNF,WRK1,WRK2)
GO TO 150
140 CALL D2FCN(NR,N,XPLS,A)
150 CONTINUE
IF(MOD(MSG/16,2).EQ.1)
+ CALL RESULT(NR,N,XPLS,FPLS,GPLS,A,P,ITNCNT,1,IPR)
C
C X <-- XPLS AND G <-- GPLS AND F <-- FPLS
C
F=FPLS
DO 160 I=1,N
X(I)=XPLS(I)
G(I)=GPLS(I)
160 CONTINUE
GO TO 100
C
C TERMINATION
C -----------
C RESET XPLS,FPLS,GPLS, IF PREVIOUS ITERATE SOLUTION
C
690 IF(ITRMCD.NE.3) GO TO 710
700 CONTINUE
FPLS=F
DO 705 I=1,N
XPLS(I)=X(I)
GPLS(I)=G(I)
705 CONTINUE
C
C PRINT RESULTS
C
710 CONTINUE
IF(MOD(MSG/8,2).EQ.0)
+ CALL RESULT(NR,N,XPLS,FPLS,GPLS,A,P,ITNCNT,0,IPR)
MSG=0
RETURN
C
900 FORMAT(14H OPTDRV ,5(E20.13,3X))
901 FORMAT(20H0OPTDRV TYPICAL X)
902 FORMAT(40H OPTDRV DIAGONAL SCALING MATRIX FOR X)
903 FORMAT(22H OPTDRV TYPICAL F =,E20.13)
904 FORMAT(40H0OPTDRV NUMBER OF GOOD DIGITS IN FCN=,I5/
+ 27H OPTDRV GRADIENT FLAG =,I5,18H (=1 IF ANALYTIC,
+ 19H GRADIENT SUPPLIED)/
+ 27H OPTDRV HESSIAN FLAG =,I5,18H (=1 IF ANALYTIC,
+ 18H HESSIAN SUPPLIED)/
+ 27H OPTDRV EXPENSE FLAG =,I5, 9H (=1 IF,
+ 45H MINIMIZATION FUNCTION EXPENSIVE TO EVALUATE)/
+ 27H OPTDRV METHOD TO USE =,I5,19H (=1,2,3 FOR LINE,
+ 49H SEARCH, DOUBLE DOGLEG, MORE-HEBDON RESPECTIVELY)/
+ 27H OPTDRV ITERATION LIMIT=,I5/
+ 27H OPTDRV MACHINE EPSILON=,E20.13)
905 FORMAT(30H0OPTDRV MAXIMUM STEP SIZE =,E20.13/
+ 30H OPTDRV STEP TOLERANCE =,E20.13/
+ 30H OPTDRV GRADIENT TOLERANCE=,E20.13/
+ 30H OPTDRV TRUST REG RADIUS =,E20.13/
+ 30H OPTDRV REL NOISE IN FCN =,E20.13/
+ 30H OPTDRV ANAL-FD TOLERANCE =,E20.13)
906 FORMAT(52H OPTDRV SHIFT FROM FORWARD TO CENTRAL DIFFERENCES,
1 14H IN ITERATION , I5)
END
SUBROUTINE OPTIF0(NR,N,X,FCN,XPLS,FPLS,GPLS,ITRMCD,A,WRK)
C
C PURPOSE
C -------
C PROVIDE SIMPLEST INTERFACE TO MINIMIZATION PACKAGE.
C USER HAS NO CONTROL OVER OPTIONS.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> INITIAL ESTIMATE OF MINIMUM
C FCN --> NAME OF ROUTINE TO EVALUATE MINIMIZATION FUNCTION.
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE.
C XPLS(N) <-- LOCAL MINIMUM
C FPLS <-- FUNCTION VALUE AT LOCAL MINIMUM XPLS
C GPLS(N) <-- GRADIENT AT LOCAL MINIMUM XPLS
C ITRMCD <-- TERMINATION CODE
C A(N,N) --> WORKSPACE
C WRK(N,9) --> WORKSPACE
C
DIMENSION X(N),XPLS(N),GPLS(N)
DIMENSION A(NR,1),WRK(NR,1)
EXTERNAL FCN,D1FCN,D2FCN
C
C EQUIVALENCE WRK(N,1) = UDIAG(N)
C WRK(N,2) = G(N)
C WRK(N,3) = P(N)
C WRK(N,4) = TYPSIZ(N)
C WRK(N,5) = SX(N)
C WRK(N,6) = WRK0(N)
C WRK(N,7) = WRK1(N)
C WRK(N,8) = WRK2(N)
C WRK(N,9) = WRK3(N)
C
CALL DFAULT(N,X,WRK(1,4),FSCALE,METHOD,IEXP,MSG,NDIGIT,
+ ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADTL,STEPMX,STEPTL)
CALL OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,WRK(1,4),FSCALE,
+ METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
+ DLT,GRADTL,STEPMX,STEPTL,
+ XPLS,FPLS,GPLS,ITRMCD,
+ A,WRK(1,1),WRK(1,2),WRK(1,3),WRK(1,5),WRK(1,6),
+ WRK(1,7),WRK(1,8),WRK(1,9))
RETURN
END
SUBROUTINE OPTIF9(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE,
+ METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
+ DLT,GRADTL,STEPMX,STEPTL,
+ XPLS,FPLS,GPLS,ITRMCD,A,WRK)
C
C PURPOSE
C -------
C PROVIDE COMPLETE INTERFACE TO MINIMIZATION PACKAGE.
C USER HAS FULL CONTROL OVER OPTIONS.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> ON ENTRY: ESTIMATE TO A ROOT OF FCN
C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C FCN: R(N) --> R(1)
C D1FCN --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE GRADIENT
C OF FCN. MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE HESSIAN OF
C OF FCN. MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X
C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION
C METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C =1 LINE SEARCH
C =2 DOUBLE DOGLEG
C =3 MORE-HEBDON
C IEXP --> =1 IF OPTIMIZATION FUNCTION FCN IS EXPENSIVE TO
C EVALUATE, =0 OTHERWISE. IF SET THEN HESSIAN WILL
C BE EVALUATED BY SECANT UPDATE INSTEAD OF
C ANALYTICALLY OR BY FINITE DIFFERENCES
C MSG <--> ON INPUT: (.GT.0) MESSAGE TO INHIBIT CERTAIN
C AUTOMATIC CHECKS
C ON OUTPUT: (.LT.0) ERROR CODE; =0 NO ERROR
C NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN
C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED
C IAHFLG --> =1 IF ANALYTIC HESSIAN SUPPLIED
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C DLT --> TRUST REGION RADIUS
C GRADTL --> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE
C ENOUGH TO ZERO TO TERMINATE ALGORITHM
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C XPLS(N) <--> ON EXIT: XPLS IS LOCAL MINIMUM
C FPLS <--> ON EXIT: FUNCTION VALUE AT SOLUTION, XPLS
C GPLS(N) <--> ON EXIT: GRADIENT AT SOLUTION XPLS
C ITRMCD <-- TERMINATION CODE
C A(N,N) --> WORKSPACE FOR HESSIAN (OR ESTIMATE)
C AND ITS CHOLESKY DECOMPOSITION
C WRK(N,8) --> WORKSPACE
C
DIMENSION X(N),XPLS(N),GPLS(N),TYPSIZ(N)
DIMENSION A(NR,1),WRK(NR,1)
EXTERNAL FCN,D1FCN,D2FCN
C
C EQUIVALENCE WRK(N,1) = UDIAG(N)
C WRK(N,2) = G(N)
C WRK(N,3) = P(N)
C WRK(N,4) = SX(N)
C WRK(N,5) = WRK0(N)
C WRK(N,6) = WRK1(N)
C WRK(N,7) = WRK2(N)
C WRK(N,8) = WRK3(N)
C
CALL OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE,
+ METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
+ DLT,GRADTL,STEPMX,STEPTL,
+ XPLS,FPLS,GPLS,ITRMCD,
+ A,WRK(1,1),WRK(1,2),WRK(1,3),WRK(1,4),WRK(1,5),
+ WRK(1,6),WRK(1,7),WRK(1,8))
RETURN
END
SUBROUTINE OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX,
+ ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,IPR,MSG)
C
C UNCONSTRAINED MINIMIZATION STOPPING CRITERIA
C --------------------------------------------
C FIND WHETHER THE ALGORITHM SHOULD TERMINATE, DUE TO ANY
C OF THE FOLLOWING:
C 1) PROBLEM SOLVED WITHIN USER TOLERANCE
C 2) CONVERGENCE WITHIN USER TOLERANCE
C 3) ITERATION LIMIT REACHED
C 4) DIVERGENCE OR TOO RESTRICTIVE MAXIMUM STEP (STEPMX) SUSPECTED
C
C
C PARAMETERS
C ----------
C N --> DIMENSION OF PROBLEM
C XPLS(N) --> NEW ITERATE X[K]
C FPLS --> FUNCTION VALUE AT NEW ITERATE F(XPLS)
C GPLS(N) --> GRADIENT AT NEW ITERATE, G(XPLS), OR APPROXIMATE
C X(N) --> OLD ITERATE X[K-1]
C ITNCNT --> CURRENT ITERATION K
C ICSCMX <--> NUMBER CONSECUTIVE STEPS .GE. STEPMX
C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C ITRMCD <-- TERMINATION CODE
C GRADTL --> TOLERANCE AT WHICH RELATIVE GRADIENT CONSIDERED CLOSE
C ENOUGH TO ZERO TO TERMINATE ALGORITHM
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION
C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IRETCD --> RETURN CODE
C MXTAKE --> BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C MSG --> IF MSG INCLUDES A TERM 8, SUPPRESS OUTPUT
C
C
INTEGER N,ITNCNT,ICSCMX,ITRMCD,ITNLIM
REAL GRADTL,STEPTL,FSCALE
REAL SX(N)
REAL XPLS(N),GPLS(N),X(N)
REAL FPLS
LOGICAL MXTAKE
C
ITRMCD=0
C
C LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X
IF(IRETCD.NE.1) GO TO 50
C IF(IRETCD.EQ.1)
C THEN
JTRMCD=3
GO TO 600
C ENDIF
50 CONTINUE
C
C FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM.
C CHECK WHETHER WITHIN TOLERANCE
C
D=AMAX1(ABS(FPLS),FSCALE)
RGX=0.0
DO 100 I=1,N
RELGRD=ABS(GPLS(I))*AMAX1(ABS(XPLS(I)),1./SX(I))/D
RGX=AMAX1(RGX,RELGRD)
100 CONTINUE
JTRMCD=1
IF(RGX.LE.GRADTL) GO TO 600
C
IF(ITNCNT.EQ.0) RETURN
C
C FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM
C CHECK WHETHER WITHIN TOLERANCE.
C
RSX=0.0
DO 120 I=1,N
RELSTP=ABS(XPLS(I)-X(I))/AMAX1(ABS(XPLS(I)),1./SX(I))
RSX=AMAX1(RSX,RELSTP)
120 CONTINUE
JTRMCD=2
IF(RSX.LE.STEPTL) GO TO 600
C
C CHECK ITERATION LIMIT
C
JTRMCD=4
IF(ITNCNT.GE.ITNLIM) GO TO 600
C
C CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX
C
IF(MXTAKE) GO TO 140
C IF(.NOT.MXTAKE)
C THEN
ICSCMX=0
RETURN
C ELSE
140 CONTINUE
IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900)
ICSCMX=ICSCMX+1
IF(ICSCMX.LT.5) RETURN
JTRMCD=5
C ENDIF
C
C
C PRINT TERMINATION CODE
C
600 ITRMCD=JTRMCD
IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD
GO TO 700
601 WRITE(IPR,901)
GO TO 700
602 WRITE(IPR,902)
GO TO 700
603 WRITE(IPR,903)
GO TO 700
604 WRITE(IPR,904)
GO TO 700
605 WRITE(IPR,905)
C
700 RETURN
C
900 FORMAT(48H0OPTSTP STEP OF MAXIMUM LENGTH (STEPMX) TAKEN)
901 FORMAT(43H0OPTSTP RELATIVE GRADIENT CLOSE TO ZERO./
+ 48H OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.)
902 FORMAT(48H0OPTSTP SUCCESSIVE ITERATES WITHIN TOLERANCE./
+ 48H OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.)
903 FORMAT(52H0OPTSTP LAST GLOBAL STEP FAILED TO LOCATE A POINT,
+ 14H LOWER THAN X./
+ 51H OPTSTP EITHER X IS AN APPROXIMATE LOCAL MINIMUM,
+ 17H OF THE FUNCTION,/
+ 50H OPTSTP THE FUNCTION IS TOO NON-LINEAR FOR THIS,
+ 11H ALGORITHM,/
+ 34H OPTSTP OR STEPTL IS TOO LARGE.)
904 FORMAT(36H0OPTSTP ITERATION LIMIT EXCEEDED./
+ 28H OPTSTP ALGORITHM FAILED.)
905 FORMAT(39H0OPTSTP MAXIMUM STEP SIZE EXCEEDED 5,
+ 19H CONSECUTIVE TIMES./
+ 50H OPTSTP EITHER THE FUNCTION IS UNBOUNDED BELOW,/
+ 47H OPTSTP BECOMES ASYMPTOTIC TO A FINITE VALUE,
+ 30H FROM ABOVE IN SOME DIRECTION,/
+ 33H OPTSTP OR STEPMX IS TOO SMALL)
END
SUBROUTINE QRAUX1(NR,N,R,I)
C
C PURPOSE
C -------
C INTERCHANGE ROWS I,I+1 OF THE UPPER HESSENBERG MATRIX R,
C COLUMNS I TO N
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF MATRIX
C R(N,N) <--> UPPER HESSENBERG MATRIX
C I --> INDEX OF ROW TO INTERCHANGE (I.LT.N)
C
DIMENSION R(NR,1)
DO 10 J=I,N
TMP=R(I,J)
R(I,J)=R(I+1,J)
R(I+1,J)=TMP
10 CONTINUE
RETURN
END
SUBROUTINE QRAUX2(NR,N,R,I,A,B)
C
C PURPOSE
C -------
C PRE-MULTIPLY R BY THE JACOBI ROTATION J(I,I+1,A,B)
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF MATRIX
C R(N,N) <--> UPPER HESSENBERG MATRIX
C I --> INDEX OF ROW
C A --> SCALAR
C B --> SCALAR
C
DIMENSION R(NR,1)
DEN=SQRT(A*A + B*B)
C=A/DEN
S=B/DEN
DO 10 J=I,N
Y=R(I,J)
Z=R(I+1,J)
R(I,J)=C*Y - S*Z
R(I+1,J)=S*Y + C*Z
10 CONTINUE
RETURN
END
SUBROUTINE QRUPDT(NR,N,A,U,V)
C
C PURPOSE
C -------
C FIND AN ORTHOGONAL (N*N) MATRIX (Q*) AND AN UPPER TRIANGULAR (N*N)
C MATRIX (R*) SUCH THAT (Q*)(R*)=R+U(V+)
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) <--> ON INPUT: CONTAINS R
C ON OUTPUT: CONTAINS (R*)
C U(N) --> VECTOR
C V(N) --> VECTOR
C
DIMENSION A(NR,1)
DIMENSION U(N),V(N)
C
C DETERMINE LAST NON-ZERO IN U(.)
C
K=N
10 IF(U(K).NE.0. .OR. K.EQ.1) GO TO 20
C IF(U(K).EQ.0. .AND. K.GT.1)
C THEN
K=K-1
GO TO 10
C ENDIF
C
C (K-1) JACOBI ROTATIONS TRANSFORM
C R + U(V+) --> (R*) + (U(1)*E1)(V+)
C WHICH IS UPPER HESSENBERG
C
20 IF(K.LE.1) GO TO 40
KM1=K-1
DO 30 II=1,KM1
I=KM1-II+1
IF(U(I).NE.0.) GO TO 25
C IF(U(I).EQ.0.)
C THEN
CALL QRAUX1(NR,N,A,I)
U(I)=U(I+1)
GO TO 30
C ELSE
25 CALL QRAUX2(NR,N,A,I,U(I),-U(I+1))
U(I)=SQRT(U(I)*U(I) + U(I+1)*U(I+1))
C ENDIF
30 CONTINUE
C ENDIF
C
C R <-- R + (U(1)*E1)(V+)
C
40 DO 50 J=1,N
A(1,J)=A(1,J) +U(1)*V(J)
50 CONTINUE
C
C (K-1) JACOBI ROTATIONS TRANSFORM UPPER HESSENBERG R
C TO UPPER TRIANGULAR (R*)
C
IF(K.LE.1) GO TO 100
KM1=K-1
DO 80 I=1,KM1
IF(A(I,I).NE.0.) GO TO 70
C IF(A(I,I).EQ.0.)
C THEN
CALL QRAUX1(NR,N,A,I)
GO TO 80
C ELSE
70 T1=A(I,I)
T2=-A(I+1,I)
CALL QRAUX2(NR,N,A,I,T1,T2)
C ENDIF
80 CONTINUE
C ENDIF
100 RETURN
END
SUBROUTINE RESULT(NR,N,X,F,G,A,P,ITNCNT,IFLG,IPR)
C
C PURPOSE
C -------
C PRINT INFORMATION
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> ITERATE X[K]
C F --> FUNCTION VALUE AT X[K]
C G(N) --> GRADIENT AT X[K]
C A(N,N) --> HESSIAN AT X[K]
C P(N) --> STEP TAKEN
C ITNCNT --> ITERATION NUMBER K
C IFLG --> FLAG CONTROLLING INFO TO PRINT
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C
DIMENSION X(N),G(N),P(N),A(NR,1)
C PRINT ITERATION NUMBER
WRITE(IPR,903) ITNCNT
IF(IFLG.EQ.0) GO TO 120
C
C PRINT STEP
WRITE(IPR,907)
WRITE(IPR,905) (P(I),I=1,N)
C
C PRINT CURRENT ITERATE
120 CONTINUE
WRITE(IPR,904)
WRITE(IPR,905) (X(I),I=1,N)
C
C PRINT FUNCTION VALUE
WRITE(IPR,906)
WRITE(IPR,905) F
C
C PRINT GRADIENT
WRITE(IPR,908)
WRITE(IPR,905) (G(I),I=1,N)
C
C PRINT HESSIAN FROM ITERATION K
IF(IFLG.EQ.0) GO TO 140
WRITE(IPR,901)
DO 130 I=1,N
WRITE(IPR,900) I
WRITE(IPR,902) (A(I,J),J=1,I)
130 CONTINUE
C
140 RETURN
900 FORMAT(15H RESULT ROW,I5)
901 FORMAT(29H RESULT HESSIAN AT X(K))
902 FORMAT(14H RESULT ,5(2X,E20.13))
903 FORMAT(/21H0RESULT ITERATE K=,I5)
904 FORMAT(18H RESULT X(K))
905 FORMAT(22H RESULT ,5(2X,E20.13) )
906 FORMAT(30H RESULT FUNCTION AT X(K))
907 FORMAT(18H RESULT STEP)
908 FORMAT(30H RESULT GRADIENT AT X(K))
END
SUBROUTINE SECFAC(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF,
+ IAGFLG,NOUPDT,S,Y,U,W)
C
C PURPOSE
C -------
C UPDATE HESSIAN BY THE BFGS FACTORED METHOD
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE, X[K-1]
C G(N) --> GRADIENT OR APPROXIMATE AT OLD ITERATE
C A(N,N) <--> ON ENTRY: CHOLESKY DECOMPOSITION OF HESSIAN IN
C LOWER PART AND DIAGONAL.
C ON EXIT: UPDATED CHOLESKY DECOMPOSITION OF HESSIAN
C IN LOWER TRIANGULAR PART AND DIAGONAL
C XPLS(N) --> NEW ITERATE, X[K]
C GPLS(N) --> GRADIENT OR APPROXIMATE AT NEW ITERATE
C EPSM --> MACHINE EPSILON
C ITNCNT --> ITERATION COUNT
C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 ITHERWISE
C NOUPDT <--> BOOLEAN: NO UPDATE YET
C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C S(N) --> WORKSPACE
C Y(N) --> WORKSPACE
C U(N) --> WORKSPACE
C W(N) --> WORKSPACE
C
REAL X(N),XPLS(N),G(N),GPLS(N)
REAL A(NR,1)
REAL S(N),Y(N),U(N),W(N)
LOGICAL NOUPDT,SKPUPD
C
IF(ITNCNT.EQ.1) NOUPDT=.TRUE.
DO 10 I=1,N
S(I)=XPLS(I)-X(I)
Y(I)=GPLS(I)-G(I)
10 CONTINUE
DEN1=SDOT(N,S,1,Y,1)
SNORM2=SNRM2(N,S,1)
YNRM2=SNRM2(N,Y,1)
IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 110
C IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2)
C THEN
CALL MVMLTU(NR,N,A,S,U)
DEN2=SDOT(N,U,1,U,1)
C
C L <-- SQRT(DEN1/DEN2)*L
C
ALP=SQRT(DEN1/DEN2)
IF(.NOT.NOUPDT) GO TO 50
C IF(NOUPDT)
C THEN
DO 30 J=1,N
U(J)=ALP*U(J)
DO 20 I=J,N
A(I,J)=ALP*A(I,J)
20 CONTINUE
30 CONTINUE
NOUPDT=.FALSE.
DEN2=DEN1
ALP=1.0
C ENDIF
50 SKPUPD=.TRUE.
C
C W = L(L+)S = HS
C
CALL MVMLTL(NR,N,A,U,W)
I=1
IF(IAGFLG.NE.0) GO TO 55
C IF(IAGFLG.EQ.0)
C THEN
RELTOL=SQRT(RNF)
GO TO 60
C ELSE
55 RELTOL=RNF
C ENDIF
60 IF(I.GT.N .OR. .NOT.SKPUPD) GO TO 70
C IF(I.LE.N .AND. SKPUPD)
C THEN
IF(ABS(Y(I)-W(I)) .LT. RELTOL*AMAX1(ABS(G(I)),ABS(GPLS(I))))
+ GO TO 65
C IF(ABS(Y(I)-W(I)) .GE. RELTOL*AMAX1(ABS(G(I)),ABS(GPLS(I))))
C THEN
SKPUPD=.FALSE.
GO TO 60
C ELSE
65 I=I+1
GO TO 60
C ENDIF
C ENDIF
70 IF(SKPUPD) GO TO 110
C IF(.NOT.SKPUPD)
C THEN
C
C W=Y-ALP*L(L+)S
C
DO 75 I=1,N
W(I)=Y(I)-ALP*W(I)
75 CONTINUE
C
C ALP=1/SQRT(DEN1*DEN2)
C
ALP=ALP/DEN1
C
C U=(L+)/SQRT(DEN1*DEN2) = (L+)S/SQRT((Y+)S * (S+)L(L+)S)
C
DO 80 I=1,N
U(I)=ALP*U(I)
80 CONTINUE
C
C COPY L INTO UPPER TRIANGULAR PART. ZERO L.
C
IF(N.EQ.1) GO TO 93
DO 90 I=2,N
IM1=I-1
DO 85 J=1,IM1
A(J,I)=A(I,J)
A(I,J)=0.
85 CONTINUE
90 CONTINUE
C
C FIND Q, (L+) SUCH THAT Q(L+) = (L+) + U(W+)
C
93 CALL QRUPDT(NR,N,A,U,W)
C
C UPPER TRIANGULAR PART AND DIAGONAL OF A NOW CONTAIN UPDATED
C CHOLESKY DECOMPOSITION OF HESSIAN. COPY BACK TO LOWER
C TRIANGULAR PART.
C
IF(N.EQ.1) GO TO 110
DO 100 I=2,N
IM1=I-1
DO 95 J=1,IM1
A(I,J)=A(J,I)
95 CONTINUE
100 CONTINUE
C ENDIF
C ENDIF
110 RETURN
END
SUBROUTINE SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT,
+ RNF,IAGFLG,NOUPDT,S,Y,T)
C
C PURPOSE
C -------
C UPDATE HESSIAN BY THE BFGS UNFACTORED METHOD
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE, X[K-1]
C G(N) --> GRADIENT OR APPROXIMATE AT OLD ITERATE
C A(N,N) <--> ON ENTRY: APPROXIMATE HESSIAN AT OLD ITERATE
C IN UPPER TRIANGULAR PART (AND UDIAG)
C ON EXIT: UPDATED APPROX HESSIAN AT NEW ITERATE
C IN LOWER TRIANGULAR PART AND DIAGONAL
C [LOWER TRIANGULAR PART OF SYMMETRIC MATRIX]
C UDIAG --> ON ENTRY: DIAGONAL OF HESSIAN
C XPLS(N) --> NEW ITERATE, X[K]
C GPLS(N) --> GRADIENT OR APPROXIMATE AT NEW ITERATE
C EPSM --> MACHINE EPSILON
C ITNCNT --> ITERATION COUNT
C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 OTHERWISE
C NOUPDT <--> BOOLEAN: NO UPDATE YET
C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C S(N) --> WORKSPACE
C Y(N) --> WORKSPACE
C T(N) --> WORKSPACE
C
REAL X(N),G(N),XPLS(N),GPLS(N)
REAL A(NR,1)
REAL UDIAG(N)
REAL S(N),Y(N),T(N)
LOGICAL NOUPDT,SKPUPD
C
C COPY HESSIAN IN UPPER TRIANGULAR PART AND UDIAG TO
C LOWER TRIANGULAR PART AND DIAGONAL
C
DO 5 J=1,N
A(J,J)=UDIAG(J)
IF(J.EQ.N) GO TO 5
JP1=J+1
DO 4 I=JP1,N
A(I,J)=A(J,I)
4 CONTINUE
5 CONTINUE
C
IF(ITNCNT.EQ.1) NOUPDT=.TRUE.
DO 10 I=1,N
S(I)=XPLS(I)-X(I)
Y(I)=GPLS(I)-G(I)
10 CONTINUE
DEN1=SDOT(N,S,1,Y,1)
SNORM2=SNRM2(N,S,1)
YNRM2=SNRM2(N,Y,1)
IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 100
C IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2)
C THEN
CALL MVMLTS(NR,N,A,S,T)
DEN2=SDOT(N,S,1,T,1)
IF(.NOT. NOUPDT) GO TO 50
C IF(NOUPDT)
C THEN
C
C H <-- [(S+)Y/(S+)HS]H
C
GAM=DEN1/DEN2
DEN2=GAM*DEN2
DO 30 J=1,N
T(J)=GAM*T(J)
DO 20 I=J,N
A(I,J)=GAM*A(I,J)
20 CONTINUE
30 CONTINUE
NOUPDT=.FALSE.
C ENDIF
50 SKPUPD=.TRUE.
C
C CHECK UPDATE CONDITION ON ROW I
C
DO 60 I=1,N
TOL=RNF*AMAX1(ABS(G(I)),ABS(GPLS(I)))
IF(IAGFLG.EQ.0) TOL=TOL/SQRT(RNF)
IF(ABS(Y(I)-T(I)).LT.TOL) GO TO 60
C IF(ABS(Y(I)-T(I)).GE.TOL)
C THEN
SKPUPD=.FALSE.
GO TO 70
C ENDIF
60 CONTINUE
70 IF(SKPUPD) GO TO 100
C IF(.NOT.SKPUPD)
C THEN
C
C BFGS UPDATE
C
DO 90 J=1,N
DO 80 I=J,N
A(I,J)=A(I,J)+Y(I)*Y(J)/DEN1-T(I)*T(J)/DEN2
80 CONTINUE
90 CONTINUE
C ENDIF
C ENDIF
100 RETURN
END
SUBROUTINE SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNOISE,STEPSZ,ANBR)
C PURPOSE
C -------
C FIND SECOND ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A"
C TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE SUBP
C "FCN" EVALUATED AT THE NEW ITERATE "XPLS"
C
C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE
C 1) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION
C IF NO ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER
C THE GRADIENT OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION
C "FCN" IS INEXPENSIVE TO EVALUATE.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C XPLS(N) --> NEW ITERATE: X[K]
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C FPLS --> FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C A(N,N) <-- FINITE DIFFERENCE APPROXIMATION TO HESSIAN
C ONLY LOWER TRIANGULAR MATRIX AND DIAGONAL
C ARE RETURNED
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C RNOISE --> RELATIVE NOISE IN FNAME [F(X)]
C STEPSZ(N) --> WORKSPACE (STEPSIZE IN I-TH COMPONENT DIRECTION)
C ANBR(N) --> WORKSPACE (NEIGHBOR IN I-TH DIRECTION)
C
C
DIMENSION XPLS(N)
DIMENSION SX(N)
DIMENSION STEPSZ(N),ANBR(N)
DIMENSION A(NR,1)
C
C FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION
C OF I-TH UNIT VECTOR.
C
OV3 = 1.0/3.0
DO 10 I=1,N
STEPSZ(I)=RNOISE**OV3 * AMAX1(ABS(XPLS(I)),1./SX(I))
XTMPI=XPLS(I)
XPLS(I)=XTMPI+STEPSZ(I)
CALL FCN(N,XPLS,ANBR(I))
XPLS(I)=XTMPI
10 CONTINUE
C
C CALCULATE COLUMN I OF A
C
DO 30 I=1,N
XTMPI=XPLS(I)
XPLS(I)=XTMPI+2.0*STEPSZ(I)
CALL FCN(N,XPLS,FHAT)
A(I,I)=((FPLS-ANBR(I))+(FHAT-ANBR(I)))/(STEPSZ(I)*STEPSZ(I))
C
C CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN
IF(I.EQ.N) GO TO 25
XPLS(I)=XTMPI+STEPSZ(I)
IP1=I+1
DO 20 J=IP1,N
XTMPJ=XPLS(J)
XPLS(J)=XTMPJ+STEPSZ(J)
CALL FCN(N,XPLS,FHAT)
A(J,I)=((FPLS-ANBR(I))+(FHAT-ANBR(J)))/(STEPSZ(I)*STEPSZ(J))
XPLS(J)=XTMPJ
20 CONTINUE
25 XPLS(I)=XTMPI
30 CONTINUE
RETURN
END
SUBROUTINE TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,
+ DLT,IRETCD,XPLSP,FPLSP,XPLS,FPLS,MXTAKE,IPR,METHOD,UDIAG)
C
C PURPOSE
C -------
C DECIDE WHETHER TO ACCEPT XPLS=X+SC AS THE NEXT ITERATE AND UPDATE THE
C TRUST REGION DLT.
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C X(N) --> OLD ITERATE X[K-1]
C F --> FUNCTION VALUE AT OLD ITERATE, F(X)
C G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE
C A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN
C LOWER TRIANGULAR PART AND DIAGONAL.
C HESSIAN OR APPROX IN UPPER TRIANGULAR PART
C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C SC(N) --> CURRENT STEP
C SX(N) --> DIAGONAL SCALING MATRIX FOR X
C NWTAKE --> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C DLT <--> TRUST REGION RADIUS
C IRETCD <--> RETURN CODE
C =0 XPLS ACCEPTED AS NEXT ITERATE;
C DLT TRUST REGION FOR NEXT ITERATION.
C =1 XPLS UNSATISFACTORY BUT ACCEPTED AS NEXT ITERATE
C BECAUSE XPLS-X .LT. SMALLEST ALLOWABLE
C STEP LENGTH.
C =2 F(XPLS) TOO LARGE. CONTINUE CURRENT ITERATION
C WITH NEW REDUCED DLT.
C =3 F(XPLS) SUFFICIENTLY SMALL, BUT QUADRATIC MODEL
C PREDICTS F(XPLS) SUFFICIENTLY WELL TO CONTINUE
C CURRENT ITERATION WITH NEW DOUBLED DLT.
C XPLSP(N) <--> WORKSPACE [VALUE NEEDS TO BE RETAINED BETWEEN
C SUCCESIVE CALLS OF K-TH GLOBAL STEP]
C FPLSP <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C XPLS(N) <-- NEW ITERATE X[K]
C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C IPR --> DEVICE TO WHICH TO SEND OUTPUT
C METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C =1 LINE SEARCH
C =2 DOUBLE DOGLEG
C =3 MORE-HEBDON
C UDIAG(N) --> DIAGONAL OF HESSIAN IN A(.,.)
C
DIMENSION X(N),XPLS(N),G(N)
DIMENSION SX(N),SC(N),XPLSP(N)
DIMENSION A(NR,1)
LOGICAL NWTAKE,MXTAKE
REAL UDIAG(N)
C
MXTAKE=.FALSE.
DO 100 I=1,N
XPLS(I)=X(I)+SC(I)
100 CONTINUE
CALL FCN(N,XPLS,FPLS)
DLTF=FPLS-F
SLP=SDOT(N,G,1,SC,1)
C
C NEXT STATEMENT ADDED FOR CASE OF COMPILERS WHICH DO NOT OPTIMIZE
C EVALUATION OF NEXT "IF" STATEMENT (IN WHICH CASE FPLSP COULD BE
C UNDEFINED).
IF(IRETCD.EQ.4) FPLSP=0.0
C$ WRITE(IPR,961) IRETCD,FPLS,FPLSP,DLTF,SLP
IF(IRETCD.NE.3 .OR. (FPLS.LT.FPLSP .AND. DLTF.LE. 1.E-4*SLP))
+ GO TO 130
C IF(IRETCD.EQ.3 .AND. (FPLS.GE.FPLSP .OR. DLTF.GT. 1.E-4*SLP))
C THEN
C
C RESET XPLS TO XPLSP AND TERMINATE GLOBAL STEP
C
IRETCD=0
DO 110 I=1,N
XPLS(I)=XPLSP(I)
110 CONTINUE
FPLS=FPLSP
DLT=.5*DLT
C$ WRITE(IPR,951)
GO TO 230
C ELSE
C
C FPLS TOO LARGE
C
130 IF(DLTF.LE. 1.E-4*SLP) GO TO 170
C IF(DLTF.GT. 1.E-4*SLP)
C THEN
C$ WRITE(IPR,952)
RLN=0.
DO 140 I=1,N
RLN=AMAX1(RLN,ABS(SC(I))/AMAX1(ABS(XPLS(I)),1./SX(I)))
140 CONTINUE
C$ WRITE(IPR,962) RLN
IF(RLN.GE.STEPTL) GO TO 150
C IF(RLN.LT.STEPTL)
C THEN
C
C CANNOT FIND SATISFACTORY XPLS SUFFICIENTLY DISTINCT FROM X
C
IRETCD=1
C$ WRITE(IPR,954)
GO TO 230
C ELSE
C
C REDUCE TRUST REGION AND CONTINUE GLOBAL STEP
C
150 IRETCD=2
DLTMP=-SLP*DLT/(2.*(DLTF-SLP))
C$ WRITE(IPR,963) DLTMP
IF(DLTMP.GE. .1*DLT) GO TO 155
C IF(DLTMP.LT. .1*DLT)
C THEN
DLT=.1*DLT
GO TO 160
C ELSE
155 DLT=DLTMP
C ENDIF
160 CONTINUE
C$ WRITE(IPR,955)
GO TO 230
C ENDIF
C ELSE
C
C FPLS SUFFICIENTLY SMALL
C
170 CONTINUE
C$ WRITE(IPR,958)
DLTFP=0.
IF (METHOD .EQ. 3) GO TO 180
C
C IF (METHOD .EQ. 2)
C THEN
C
DO 177 I = 1, N
TEMP = 0.0
DO 173 J = I, N
TEMP = TEMP + (A(J, I)*SC(J))
173 CONTINUE
DLTFP = DLTFP + TEMP*TEMP
177 CONTINUE
GO TO 190
C
C ELSE
C
180 DO 187 I = 1, N
DLTFP = DLTFP + UDIAG(I)*SC(I)*SC(I)
IF (I .EQ. N) GO TO 187
TEMP = 0
IP1 = I + 1
DO 183 J = IP1, N
TEMP = TEMP + A(I, J)*SC(I)*SC(J)
183 CONTINUE
DLTFP = DLTFP + 2.0*TEMP
187 CONTINUE
C
C END IF
C
190 DLTFP = SLP + DLTFP/2.0
C$ WRITE(IPR,964) DLTFP,NWTAKE
IF(IRETCD.EQ.2 .OR. (ABS(DLTFP-DLTF).GT. .1*ABS(DLTF))
+ .OR. NWTAKE .OR. (DLT.GT. .99*STEPMX)) GO TO 210
C IF(IRETCD.NE.2 .AND. (ABS(DLTFP-DLTF) .LE. .1*ABS(DLTF))
C + .AND. (.NOT.NWTAKE) .AND. (DLT.LE. .99*STEPMX))
C THEN
C
C DOUBLE TRUST REGION AND CONTINUE GLOBAL STEP
C
IRETCD=3
DO 200 I=1,N
XPLSP(I)=XPLS(I)
200 CONTINUE
FPLSP=FPLS
DLT=AMIN1(2.*DLT,STEPMX)
C$ WRITE(IPR,959)
GO TO 230
C ELSE
C
C ACCEPT XPLS AS NEXT ITERATE. CHOOSE NEW TRUST REGION.
C
210 CONTINUE
C$ WRITE(IPR,960)
IRETCD=0
IF(DLT.GT. .99*STEPMX) MXTAKE=.TRUE.
IF(DLTF.LT. .1*DLTFP) GO TO 220
C IF(DLTF.GE. .1*DLTFP)
C THEN
C
C DECREASE TRUST REGION FOR NEXT ITERATION
C
DLT=.5*DLT
GO TO 230
C ELSE
C
C CHECK WHETHER TO INCREASE TRUST REGION FOR NEXT ITERATION
C
220 IF(DLTF.LE. .75*DLTFP) DLT=AMIN1(2.*DLT,STEPMX)
C ENDIF
C ENDIF
C ENDIF
C ENDIF
230 CONTINUE
C$ WRITE(IPR,953)
C$ WRITE(IPR,956) IRETCD,MXTAKE,DLT,FPLS
C$ WRITE(IPR,957)
C$ WRITE(IPR,965) (XPLS(I),I=1,N)
RETURN
C
951 FORMAT(55H TREGUP RESET XPLS TO XPLSP. TERMINATION GLOBAL STEP)
952 FORMAT(26H TREGUP FPLS TOO LARGE.)
953 FORMAT(38H0TREGUP VALUES AFTER CALL TO TREGUP)
954 FORMAT(54H TREGUP CANNOT FIND SATISFACTORY XPLS DISTINCT FROM,
+ 27H X. TERMINATE GLOBAL STEP.)
955 FORMAT(53H TREGUP REDUCE TRUST REGION. CONTINUE GLOBAL STEP.)
956 FORMAT(21H TREGUP IRETCD=,I3/
+ 21H TREGUP MXTAKE=,L1/
+ 21H TREGUP DLT =,E20.13/
+ 21H TREGUP FPLS =,E20.13)
957 FORMAT(32H TREGUP NEW ITERATE (XPLS))
958 FORMAT(35H TREGUP FPLS SUFFICIENTLY SMALL.)
959 FORMAT(54H TREGUP DOUBLE TRUST REGION. CONTINUE GLOBAL STEP.)
960 FORMAT(50H TREGUP ACCEPT XPLS AS NEW ITERATE. CHOOSE NEW,
+ 38H TRUST REGION. TERMINATE GLOBAL STEP.)
961 FORMAT(18H TREGUP IRETCD=,I5/
+ 18H TREGUP FPLS =,E20.13/
+ 18H TREGUP FPLSP =,E20.13/
+ 18H TREGUP DLTF =,E20.13/
+ 18H TREGUP SLP =,E20.13)
962 FORMAT(18H TREGUP RLN =,E20.13)
963 FORMAT(18H TREGUP DLTMP =,E20.13)
964 FORMAT(18H TREGUP DLTFP =,E20.13/
+ 18H TREGUP NWTAKE=,L1)
965 FORMAT(14H TREGUP ,5(E20.13,3X))
END
FUNCTION EPSMCH(DUM)
C
C PURPOSE
C -------
C CALCULATE MACHINE EPSILON
C
C PARAMETERS
C ----------
C EPSMCH <-- MACHINE EPSILON
C
DATA ONE/1.0/
DATA HALF/.5/
C
DEL=ONE
10 CONTINUE
EPSMCH=DEL
DEL=DEL*HALF
IF(ONE+DEL.GT.ONE) GO TO 10
RETURN
END
SUBROUTINE MVMLTL(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=LX
C WHERE L IS A LOWER TRIANGULAR MATRIX STORED IN A
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) --> LOWER TRIANGULAR (N*N) MATRIX
C X(N) --> OPERAND VECTOR
C Y(N) <-- RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE
C
DIMENSION A(NR,1),X(N),Y(N)
DO 30 I=1,N
SUM=0.
DO 10 J=1,I
SUM=SUM+A(I,J)*X(J)
10 CONTINUE
Y(I)=SUM
30 CONTINUE
RETURN
END
SUBROUTINE MVMLTS(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=AX
C WHERE "A" IS A SYMMETRIC (N*N) MATRIX STORED IN ITS LOWER
C TRIANGULAR PART AND X,Y ARE N-VECTORS
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(N,N) --> SYMMETRIC (N*N) MATRIX STORED IN
C LOWER TRIANGULAR PART AND DIAGONAL
C X(N) --> OPERAND VECTOR
C Y(N) <-- RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE.
C
DIMENSION A(NR,1),X(N),Y(N)
DO 30 I=1,N
SUM=0.
DO 10 J=1,I
SUM=SUM+A(I,J)*X(J)
10 CONTINUE
IF(I.EQ.N) GO TO 25
IP1=I+1
DO 20 J=IP1,N
SUM=SUM+A(J,I)*X(J)
20 CONTINUE
25 Y(I)=SUM
30 CONTINUE
RETURN
END
SUBROUTINE MVMLTU(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=(L+)X
C WHERE L IS A LOWER TRIANGULAR MATRIX STORED IN A
C (L-TRANSPOSE (L+) IS TAKEN IMPLICITLY)
C
C PARAMETERS
C ----------
C NR --> ROW DIMENSION OF MATRIX
C N --> DIMENSION OF PROBLEM
C A(NR,1) --> LOWER TRIANGULAR (N*N) MATRIX
C X(N) --> OPERAND VECTOR
C Y(N) <-- RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE
C
DIMENSION A(NR,1),X(N),Y(N)
DO 30 I=1,N
SUM=0.
DO 10 J=I,N
SUM=SUM+A(J,I)*X(J)
10 CONTINUE
Y(I)=SUM
30 CONTINUE
RETURN
END
SUBROUTINE SCLMUL(N,S,V,Z)
C
C PURPOSE
C -------
C MULTIPLY VECTOR BY SCALAR
C RESULT VECTOR MAY BE OPERAND VECTOR
C
C PARAMETERS
C ----------
C N --> DIMENSION OF VECTORS
C S --> SCALAR
C V(N) --> OPERAND VECTOR
C Z(N) <-- RESULT VECTOR
DIMENSION V(N),Z(N)
DO 100 I=1,N
Z(I)=S*V(I)
100 CONTINUE
RETURN
END
REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
C RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY.
C SDOT = SUM FOR I = 0 TO N-1 OF SX(LX+I*INCX) * SY(LY+I*INCY),
C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C DEFINED IN A SIMILAR WAY USING INCY.
C
REAL SX(1),SY(1)
SDOT = 0.0E0
IF(N.LE.0)RETURN
IF(INCX.EQ.INCY) IF(INCX-1)5,20,60
5 CONTINUE
C
C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS.
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
SDOT = SDOT + SX(IX)*SY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
20 M = MOD(N,5)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
SDOT = SDOT + SX(I)*SY(I)
30 CONTINUE
IF( N .LT. 5 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
* SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
50 CONTINUE
RETURN
C
C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
60 CONTINUE
NS=N*INCX
DO 70 I=1,NS,INCX
SDOT = SDOT + SX(I)*SY(I)
70 CONTINUE
RETURN
END
REAL FUNCTION SNRM2 ( N, SX, INCX)
INTEGER NEXT
REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
DATA ZERO, ONE /0.0E0, 1.0E0/
C
C EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE
C INCREMENT INCX .
C IF N .LE. 0 RETURN WITH RESULT = 0.
C IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C C.L.LAWSON, 1978 JAN 08
C
C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
C HOPEFULLY APPLICABLE TO ALL MACHINES.
C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES.
C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES.
C WHERE
C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
C V = LARGEST NO. (OVERFLOW LIMIT)
C
C BRIEF OUTLINE OF ALGORITHM..
C
C PHASE 1 SCANS ZERO COMPONENTS.
C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C VALUES FOR CUTLO AND CUTHI..
C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
C UNIVAC AND DEC AT 2**(-103)
C THUS CUTLO = 2**(-51) = 4.44089E-16
C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C THUS CUTHI = 2**(63.5) = 1.30438E19
C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C THUS CUTLO = 2**(-33.5) = 8.23181D-11
C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19
C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
C
IF(N .GT. 0) GO TO 10
SNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C BEGIN MAIN LOOP
I = 1
20 GO TO NEXT,(30, 50, 70, 110)
30 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
XMAX = ZERO
C
C PHASE 1. SUM IS ZERO
C
50 IF( SX(I) .EQ. ZERO) GO TO 200
IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
C
C PREPARE FOR PHASE 2.
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
ASSIGN 110 TO NEXT
SUM = (SUM / SX(I)) / SX(I)
105 XMAX = ABS(SX(I))
GO TO 115
C
C PHASE 2. SUM IS SMALL.
C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70 IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75
C
C COMMON CODE FOR PHASES 2 AND 4.
C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
C
110 IF( ABS(SX(I)) .LE. XMAX ) GO TO 115
SUM = ONE + SUM * (XMAX / SX(I))**2
XMAX = ABS(SX(I))
GO TO 200
C
115 SUM = SUM + (SX(I)/XMAX)**2
GO TO 200
C
C
C PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI/FLOAT( N )
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J =I,NN,INCX
IF(ABS(SX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + SX(J)**2
SNRM2 = SQRT( SUM )
GO TO 300
C
200 CONTINUE
I = I + INCX
IF ( I .LE. NN ) GO TO 20
C
C END OF MAIN LOOP.
C
C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
SNRM2 = XMAX * SQRT(SUM)
300 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
INTEGER NEXT
DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
DATA ZERO, ONE /0.0D0, 1.0D0/
C
C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C INCREMENT INCX .
C IF N .LE. 0 RETURN WITH RESULT = 0.
C IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C C.L.LAWSON, 1978 JAN 08
C
C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
C HOPEFULLY APPLICABLE TO ALL MACHINES.
C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES.
C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES.
C WHERE
C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
C V = LARGEST NO. (OVERFLOW LIMIT)
C
C BRIEF OUTLINE OF ALGORITHM..
C
C PHASE 1 SCANS ZERO COMPONENTS.
C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C VALUES FOR CUTLO AND CUTHI..
C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
C UNIVAC AND DEC AT 2**(-103)
C THUS CUTLO = 2**(-51) = 4.44089E-16
C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C THUS CUTHI = 2**(63.5) = 1.30438E19
C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C THUS CUTLO = 2**(-33.5) = 8.23181D-11
C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19
C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C
IF(N .GT. 0) GO TO 10
DNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C BEGIN MAIN LOOP
I = 1
20 GO TO NEXT,(30, 50, 70, 110)
30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
XMAX = ZERO
C
C PHASE 1. SUM IS ZERO
C
50 IF( DX(I) .EQ. ZERO) GO TO 200
IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C PREPARE FOR PHASE 2.
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
ASSIGN 110 TO NEXT
SUM = (SUM / DX(I)) / DX(I)
105 XMAX = DABS(DX(I))
GO TO 115
C
C PHASE 2. SUM IS SMALL.
C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C COMMON CODE FOR PHASES 2 AND 4.
C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
C
110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
SUM = ONE + SUM * (XMAX / DX(I))**2
XMAX = DABS(DX(I))
GO TO 200
C
115 SUM = SUM + (DX(I)/XMAX)**2
GO TO 200
C
C
C PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI/FLOAT( N )
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J =I,NN,INCX
IF(DABS(DX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + DX(J)**2
DNRM2 = DSQRT( SUM )
GO TO 300
C
200 CONTINUE
I = I + INCX
IF ( I .LE. NN ) GO TO 20
C
C END OF MAIN LOOP.
C
C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
DNRM2 = XMAX * DSQRT(SUM)
300 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C RETURNS THE DOT PRODUCT OF DOUBLE PRECISION DX AND DY.
C DDOT = SUM FOR I = 0 TO N-1 OF DX(LX+I*INCX) * DY(LY+I*INCY)
C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C DEFINED IN A SIMILAR WAY USING INCY.
C
DOUBLE PRECISION DX(1),DY(1)
DDOT = 0.D0
IF(N.LE.0)RETURN
IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
5 CONTINUE
C
C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DDOT = DDOT + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1.
C
C
C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
20 M = MOD(N,5)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DDOT = DDOT + DX(I)*DY(I)
30 CONTINUE
IF( N .LT. 5 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
$ DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
50 CONTINUE
RETURN
C
C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
60 CONTINUE
NS = N*INCX
DO 70 I=1,NS,INCX
DDOT = DDOT + DX(I)*DY(I)
70 CONTINUE
RETURN
END
|
