SUBROUTINE UNCMND (N,X0,FCN,X,F,INFO,W,LW) C***BEGIN PROLOGUE UNCMND C***DATE WRITTEN 870923 (YYMMDD) C***REVISION DATE 871222 (YYMMDD) C***CATEGORY NO. G1B1A1 C***KEYWORDS UNCONSTRAINED MINIMIZATION C***AUTHOR NASH, S.G., (GEORGE MASON UNIVERSITY) C***PURPOSE UNCMND minimizes a smooth nonlinear function of n variables. C A subroutine that computes the function value at any point C must be supplied, but derivative values are not required. C UNCMND provides a simple interface to more flexible lower C level routines. User has no control over options. C C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C This routine uses a quasi-Newton algorithm with line search C to minimize the function represented by the subroutine FCN. C At each iteration, the nonlinear function is approximated C by a quadratic function derived from a Taylor series. C The quadratic function is minimized to obtain a search direction, C and an approximate minimum of the nonlinear function along C the search direction is found using a line search. The C algorithm computes an approximation to the second derivative C matrix of the nonlinear function using quasi-Newton techniques. C C The UNCMND package is quite general, and provides many options C for the user. However, this subroutine is designed to be C easy to use, with few choices allowed. For example: C C 1. Only function values need be computed. First derivative C values are obtained by finite-differencing. This can be C very costly when the number of variables is large. C C 2. It is assumed that the function values can be obtained C accurately (to an accuracy comparable to the precision of C the computer arithmetic). C C 3. At most 150 iterations are allowed. C C 4. It is assumed that the function values are well-scaled, C that is, that the optimal function value is not pathologically C large or small. C C For more information, see the reference listed below. C C PARAMETERS C ---------- C N --> INTEGER C Dimension of problem C X0(N) --> DOUBLE PRECISION C Initial estimate of minimum C FCN --> Name of routine to evaluate minimization function. C Must be declared EXTERNAL in calling routine, and C have calling sequence C SUBROUTINE FCN(N, X, F) C with N and X as here, F the computed function value. C X(N) <-- DOUBLE PRECISION C Local minimum C F <-- DOUBLE PRECISION C Function value at local minimum X C INFO <-- INTEGER C Termination code C INFO = 0: Optimal solution found C INFO = 1: Terminated with gradient small, C X is probably optimal C INFO = 2: Terminated with stepsize small, C X is probably optimal C INFO = 3: Lower point cannot be found, C X is probably optimal C INFO = 4: Iteration limit (150) exceeded C INFO = 5: Too many large steps, C function may be unbounded C INFO = -1: Insufficient workspace C W(LW) --> DOUBLE PRECISION C Workspace C LW --> INTEGER C Size of workspace, at least N*(N+10) C C***REFERENCES R.B. SCHNABEL, J.E. KOONTZ, AND BE.E. WEISS, A MODULAR C SYSTEM OF ALGORITHMS FOR UNCONSTRAINED MINIMIZATION, C REPORT CU-CS-240-82, COMP. SCI. DEPT., UNIV. OF C COLORADO AT BOULDER, 1982. C***ROUTINES CALLED OPTDRD, XERROR C***END PROLOGUE UNCMND IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X0(N),X(N),W(LW) CHARACTER ERRMSG*68 EXTERNAL FCN, D1FCND, D2FCND C---------------------------------------------------------------- C SUBDIVIDE WORKSPACE C---------------------------------------------------------------- C***FIRST EXECUTABLE STATEMENT UNCMND IG = 1 IT = IG + N IW1 = IT + N IW2 = IW1 + N IW3 = IW2 + N IW4 = IW3 + N IW5 = IW4 + N IW6 = IW5 + N IW7 = IW6 + N IW8 = IW7 + N IA = IW8 + N LWMIN = IA + N*N-1 IF (LWMIN .GT. LW) THEN INFO = -1 WRITE(ERRMSG, '( * ''UNCMND ERROR (INFO=-1) -- INSUFFICIENT WORKSPACE'', * '', LW = '', I5 )' ) LW CALL XERROR(ERRMSG(1:60), 60, -1, 0) RETURN END IF C---------------------------------------------------------------- C SET UP PARAMETERS FOR OPTIF9 C---------------------------------------------------------------- C PARAMETERS THAT SHOULD NOT BE CHANGED C C NR = PARAMETER USED TO DIVIDE WORKSPACE C METHOD = 1 (LINE SEARCH) -- DO NOT CHANGE C MSG = 9 => NO PRINTING, N=1 ALLOWED C IAGFLG = 1 => ANALYTIC GRADIENT SUPPLIED (0 OTHERWISE) C IAHFLG = 1 => ANALYTIC HESSIAN SUPPLIED (0 OTHERWISE) C IPR = DEVICE FOR OUTPUT (IRRELEVANT IN CURRENT VERSION) C DLT = (IRRELEVANT PARAMETER FOR METHOD = 1) C EPSM = MACHINE EPSILON C IEXP = 1 => FUNCTION EXPENSIVE TO EVALUATE (IEXP = 0 => CHEAP) C NR = N METHOD = 1 MSG = 9 IAGFLG = 0 IAHFLG = 0 IPR = 0 DLT = -1.0D0 EPSM = D1MACH(4) IEXP = 1 C C PARAMETERS THAT MAY BE CHANGED: C C NDIGIT = -1 => OPTIF9 ASSUMES F IS FULLY ACCURATE C ITNLIM = 150 = MAXIMUM NUMBER OF ITERATIONS ALLOWED C GRADTL = ZERO TOLERANCE FOR GRADIENT, FOR CONVERGENCE TESTS C STEPMX = MAXIMUM ALLOWABLE STEP SIZE C STEPTL = ZERO TOLERANCE FOR STEP, FOR CONVERGENCE TESTS C FSCALE = TYPICAL ORDER-OF-MAGNITUDE SIZE OF FUNCTION C TYPSIZ = TYPICAL ORDER-OF-MAGNITUDE SIZE OF X (STORED IN W(LT)) C NDIGIT = -1 ITNLIM = 150 GRADTL = EPSM**(1.0D0/3.0D0) STEPMX = 1000.0D0 STEPTL = SQRT(EPSM) FSCALE = 1.0D0 DO 10 LT = IT,IT+N-1 W(LT) = 1.0D0 10 CONTINUE C C MINIMIZE FUNCTION C CALL OPTDRD (NR, N, X0, FCN, D1FCND, D2FCND, W(IT), FSCALE, + METHOD, IEXP, MSG, NDIGIT, ITNLIM, IAGFLG, IAHFLG, + IPR, DLT, GRADTL, STEPMX, STEPTL, + X, F, W(IG), INFO, W(IA), + W(IW1), W(IW2), W(IW3), W(IW4), + W(IW5), W(IW6), W(IW7), W(IW8)) C IF (INFO .EQ. 1) THEN WRITE(ERRMSG, '( * ''UNCMND WARNING -- INFO = 1'', * '': PROBABLY CONVERGED, GRADIENT SMALL'')' ) CALL XERROR(ERRMSG(1:62), 62, INFO, 0) END IF IF (INFO .EQ. 2) THEN WRITE(ERRMSG, '( * ''UNCMND WARNING -- INFO = 2'', * '': PROBABLY CONVERGED, STEPSIZE SMALL'')' ) CALL XERROR(ERRMSG(1:62), 62, INFO, 0) END IF IF (INFO .EQ. 3) THEN WRITE(ERRMSG, '( * ''UNCMND WARNING -- INFO = 3'', * '': CANNOT FIND LOWER POINT'')' ) CALL XERROR(ERRMSG(1:51), 51, INFO, 0) END IF IF (INFO .EQ. 4) THEN WRITE(ERRMSG, '( * ''UNCMND WARNING -- INFO = 4'', * '': TOO MANY ITERATIONS'')' ) CALL XERROR(ERRMSG(1:47), 47, INFO, 0) END IF IF (INFO .EQ. 5) THEN WRITE(ERRMSG, '( * ''UNCMND WARNING -- INFO = 5'', * '': TOO MANY LARGE STEPS, POSSIBLY UNBOUNDED'')' ) CALL XERROR(ERRMSG(1:68), 68, INFO, 0) END IF C RETURN END SUBROUTINE DUMMYD C*DUMMY SUBROUTINE (REPRESENTS UNUSED ROUTINES FROM FULL VERSION)* RETURN END SUBROUTINE BAKSLD(NR,N,A,X,B) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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.D0 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 D1FCND C*DUMMY* RETURN END SUBROUTINE D2FCND C*DUMMY* RETURN END SUBROUTINE FORSLD(NR,N,A,X,B) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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.0D0 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 FSTCDD (N, X, FCN, SX, RNOISE, G) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N) DIMENSION SX(N) DIMENSION 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.0D0/3.0D0 DO 10 I = 1, N STEPI = RNOISE**THIRD * MAX(ABS(X(I)), 1.0D0/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.0D0*STEPI) 10 CONTINUE RETURN END SUBROUTINE FSTFDD(NR,M,N,XPLS,FCN,FPLS,A,SX,RNOISE,FHAT,ICASE) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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)*MAX(ABS(XPLS(J)),1.D0/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.0D0 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE HSNNTD(NR,N,A,SX,METHOD) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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.D0 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE LLTSLD(NR,N,A,X,B) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NR,1),X(N),B(N) C C FORWARD SOLVE, RESULT IN X C CALL FORSLD(NR,N,A,X,B) C C BACK SOLVE, RESULT IN X C CALL BAKSLD(NR,N,A,X,X) RETURN END SUBROUTINE LNSRCD(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE, + IRETCD,STEPMX,STEPTL,SX,IPR) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N,IRETCD DIMENSION SX(N) DIMENSION X(N),G(N),P(N) DIMENSION XPLS(N) LOGICAL MXTAKE C IPR=IPR MXTAKE=.FALSE. IRETCD=2 TMP=0.0D0 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 SCLMLD(N,SCL,P,P) SLN=STEPMX 10 CONTINUE SLP=DDOT(N,G,1,P,1) RLN=0.D0 DO 15 I=1,N RLN=MAX(RLN,ABS(P(I))/MAX(ABS(X(I)),1.D0/SX(I))) 15 CONTINUE RMNLMB=STEPTL/RLN ALMBDA=1.0D0 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) IF(FPLS.GT. F+SLP*1.D-4*ALMBDA) GO TO 130 C IF(FPLS.LE. F+SLP*1.D-4*ALMBDA) C THEN C C SOLUTION FOUND C IRETCD=0 IF(ALMBDA.EQ.1.0D0 .AND. SLN.GT. .99D0*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.0D0) GO TO 150 C IF(ALMBDA.EQ.1.0D0) C THEN C C FIRST BACKTRACK: QUADRATIC FIT C TLMBDA=-SLP/(2.D0*(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.0D0/(ALMBDA-PLMBDA) A=T3*(T1/(ALMBDA*ALMBDA) - T2/(PLMBDA*PLMBDA)) B=T3*(T2*ALMBDA/(PLMBDA*PLMBDA) + - T1*PLMBDA/(ALMBDA*ALMBDA) ) DISC=B*B-3.0D0*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.0D0,A)*SQRT(DISC))/(3.0D0*A) GO TO 165 C ELSE C C BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM C 160 TLMBDA=(-B-SIGN(1.0D0,A)*SQRT(DISC))/(3.0D0*A) C ENDIF 165 IF(TLMBDA.GT. .5D0*ALMBDA) TLMBDA=.5D0*ALMBDA C ENDIF 170 PLMBDA=ALMBDA PFPLS=FPLS IF(TLMBDA.GE. ALMBDA*.1D0) GO TO 180 C IF(TLMBDA.LT.ALMBDA/10.) C THEN ALMBDA=ALMBDA*.1D0 GO TO 190 C ELSE 180 ALMBDA=TLMBDA C ENDIF C ENDIF C ENDIF 190 GO TO 100 END SUBROUTINE MVMLLD(NR,N,A,X,Y) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NR,1),X(N),Y(N) DO 30 I=1,N SUM=0.D0 DO 10 J=1,I SUM=SUM+A(I,J)*X(J) 10 CONTINUE Y(I)=SUM 30 CONTINUE RETURN END SUBROUTINE MVMLUD(NR,N,A,X,Y) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NR,1),X(N),Y(N) DO 30 I=1,N SUM=0.D0 DO 10 J=I,N SUM=SUM+A(J,I)*X(J) 10 CONTINUE Y(I)=SUM 30 CONTINUE RETURN END SUBROUTINE OPTCHD(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM, + DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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. IPR=IPR 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.D0) TYPSIZ(I)=1.0D0 IF(TYPSIZ(I).LT.0.D0) TYPSIZ(I)=-TYPSIZ(I) SX(I)=1.0D0/TYPSIZ(I) 10 CONTINUE C C CHECK MAXIMUM STEP SIZE C IF (STEPMX .GT. 0.0D0) GO TO 20 STPSIZ = 0.0D0 DO 15 I = 1, N STPSIZ = STPSIZ + X(I)*X(I)*SX(I)*SX(I) 15 CONTINUE STPSIZ = SQRT(STPSIZ) STEPMX = MAX(1.0D3*STPSIZ, 1.0D3) 20 CONTINUE C CHECK FUNCTION SCALE IF(FSCALE.EQ.0.D0) FSCALE=1.0D0 IF(FSCALE.LT.0.D0) FSCALE=-FSCALE C C CHECK GRADIENT TOLERANCE IF(GRADTL.LT.0.D0) 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=-LOG10(EPSM) C C CHECK TRUST REGION RADIUS IF(DLT.LE.0.D0) DLT=-1.0D0 IF (DLT .GT. STEPMX) DLT = STEPMX RETURN C C ERROR EXITS C C INSERT CALLS TO XERROR C OPTCHD: ILLEGAL DIMENSION, N 805 MSG=-1 GO TO 895 C OPTCHD: WARNING +++ THIS PACKAGE IS INEFFICIENT FOR N=1 C (CANNOT BE GENERATED WITH CURRENT VALUE OF MSG) 810 MSG=-2 GO TO 895 C OPTCHD: ILLEGAL TOLERANCE. GRADTL 815 MSG=-3 GO TO 895 C OPTCHD: ILLEGAL ITERATION LIMIT. ITNLIM 820 MSG=-4 GO TO 895 C OPTCHD: MINIMIZATION FUNCTION HAS NO GOOD DIGITS 825 MSG=-5 GO TO 895 C OPTCHD: USER REQUESTS THAT ANALYTIC GRADIENT BE USED C BUT ANALYTIC GRADIENT NOT SUPPLIED 830 MSG=-6 GO TO 895 C OPTCHD: USER REQUESTS THAT ANALYTIC HESSIAN BE USED C BUT ANALYTIC HESSIAN NOT SUPPLIED 835 MSG=-7 895 RETURN END SUBROUTINE OPTDRD(NR,N,X,FCN,D1FCND,D2FCND,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 THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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,D1FCND,D2FCND C C INITIALIZATION C -------------- DO 10 I=1,N P(I)=0.D0 10 CONTINUE ITNCNT=0 IRETCD=-1 EPSM=D1MACH(4) CALL OPTCHD(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM, + DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR) IF(MSG.LT.0) RETURN RNF=MAX(10.0D0**(-NDIGIT),EPSM) WRITE(*,*) 'EPSM',EPSM WRITE(*,*) 'RNF',RNF WRITE(*,*) 'NDIGIT',NDIGIT WRITE(*,*) 'GRADTL',GRADTL WRITE(*,*) 'FSCALE',FSCALE WRITE(*,*) 'STEPMX',STEPMX WRITE(*,*) 'TYPSIZ(1)',TYPSIZ(1) WRITE(*,*) 'TYPSIZ(2)',TYPSIZ(2) WRITE(*,*) 'SX(1)',SX(1) WRITE(*,*) 'SX(2)',SX(2) 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 FSTFDD (1, 1, N, X, FCN, F, G, SX, RNF, WRK, 1) GO TO 25 C 20 CALL D1FCND IF (MOD(MSG/2,2) .EQ. 1) GO TO 25 C IF (MOD(MSG/2,2).EQ.0) C THEN CALL DUMMYD IF (MSG .LT. 0) RETURN 25 CONTINUE C CALL OPTSTD(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 HSNNTD(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 D2FCND 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 FSTFDD (NR, N, N, X, D1FCND, G, A, SX, 1 RNF, WRK1, 3) IF (IAGFLG .NE. 1) CALL DUMMYD 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 D2FCND GO TO 88 C C ELSE 85 CALL DUMMYD C C HESCHK EVALUATES D2FCND AND CHECKS IT AGAINST THE FINITE C DIFFERENCE HESSIAN WHICH IT CALCULATES BY CALLING FSTFDD 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 DUMMYD 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 SECFCD.) C IF(IEXP.EQ.1 .AND. METHOD.NE.3) GO TO 105 103 CALL DUMMYD 105 CONTINUE C C SOLVE FOR NEWTON STEP: AP=-G C DO 110 I=1,N WRK1(I)=-G(I) 110 CONTINUE CALL LLTSLD(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 LNSRCD(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE,IRETCD, + STEPMX,STEPTL,SX,IPR) IF(METHOD.EQ.2) + CALL DUMMYD IF(METHOD.EQ.3) + CALL DUMMYD 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 C CALL FSTCDD (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 D1FCND GO TO 120 C C CENTRAL DIFFERENCE GRADIENT 116 CALL FSTCDD (N, XPLS, FCN, SX, RNF, GPLS) GO TO 120 C C FORWARD DIFFERENCE GRADIENT 118 CALL FSTFDD (1, 1, N, XPLS, FCN, FPLS, GPLS, SX, RNF, WRK, 1) 120 CONTINUE C C CHECK WHETHER STOPPING CRITERIA SATISFIED C CALL OPTSTD(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 DUMMYD IF(METHOD.NE.3) + CALL SECFCD(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 FSTFDD(NR,N,N,XPLS,D1FCND,GPLS,A,SX,RNF,WRK1,3) IF(IAGFLG.NE.1) CALL DUMMYD GO TO 150 140 CALL D2FCND 150 CONTINUE IF(MOD(MSG/16,2).EQ.1) + CALL DUMMYD 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 DUMMYD MSG=0 RETURN C END SUBROUTINE OPTSTD(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX, + ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,IPR,MSG) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N,ITNCNT,ICSCMX,ITRMCD,ITNLIM DIMENSION SX(N) DIMENSION XPLS(N),GPLS(N),X(N) LOGICAL MXTAKE C IPR=IPR 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=MAX(ABS(FPLS),FSCALE) RGX=0.0D0 DO 100 I=1,N RELGRD=ABS(GPLS(I))*MAX(ABS(XPLS(I)),1.D0/SX(I))/D RGX=MAX(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.0D0 DO 120 I=1,N RELSTP=ABS(XPLS(I)-X(I))/MAX(ABS(XPLS(I)),1.D0/SX(I)) RSX=MAX(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 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 C INSERT CALLS TO XERROR C OPTSTD: RELATIVE GRADIENT CLOSE TO ZERO. C CURRENT ITERATE IS PROBABLY SOLUTION. 601 CONTINUE GO TO 700 C OPTSTD: SUCCESSIVE ITERATES WITHIN TOLERANCE. C CURRENT ITERATE IS PROBABLY SOLUTION. 602 CONTINUE GO TO 700 C OPTSTD: LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X. C EITHER X IS AN APPROXIMATE LOCAL MINIMUM OF THE FUNCTION, C THE FUNCTION IS TOO NON-LINEAR FOR THIS ALGORITHM, C OR STEPTL IS TOO LARGE. 603 CONTINUE GO TO 700 C OPTSTD: ITERATION LIMIT EXCEEDED. ALGORITHM FAILED. 604 CONTINUE GO TO 700 C OPTSTD: MAXIMUM STEP SIZE EXCEEDED 5 CONSECUTIVE TIMES. C EITHER THE FUNCTION IS UNBOUNDED BELOW, C BECOMES ASYMPTOTIC TO A FINITE VALUE FROM ABOVE, C OR STEPMX IS TOO SMALL 605 CONTINUE C 700 RETURN C END SUBROUTINE QRAX1D(NR,N,R,I) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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 QRAX2D(NR,N,R,I,A,B) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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 QRUPDD(NR,N,A,U,V) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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.D0 .OR. K.EQ.1) GO TO 20 C IF(U(K).EQ.0.D0 .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.D0) GO TO 25 C IF(U(I).EQ.0.D0) C THEN CALL QRAX1D(NR,N,A,I) U(I)=U(I+1) GO TO 30 C ELSE 25 CALL QRAX2D(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.D0) GO TO 70 C IF(A(I,I).EQ.0.D0) C THEN CALL QRAX1D(NR,N,A,I) GO TO 80 C ELSE 70 T1=A(I,I) T2=-A(I+1,I) CALL QRAX2D(NR,N,A,I,T1,T2) C ENDIF 80 CONTINUE C ENDIF 100 RETURN END SUBROUTINE SCLMLD(N,S,V,Z) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION V(N),Z(N) DO 100 I=1,N Z(I)=S*V(I) 100 CONTINUE RETURN END SUBROUTINE SECFCD(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF, + IAGFLG,NOUPDT,S,Y,U,W) C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),XPLS(N),G(N),GPLS(N) DIMENSION A(NR,1) DIMENSION 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=DDOT(N,S,1,Y,1) SNORM2=DNRM2(N,S,1) YNRM2=DNRM2(N,Y,1) IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 110 C IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2) C THEN CALL MVMLUD(NR,N,A,S,U) DEN2=DDOT(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.0D0 C ENDIF 50 SKPUPD=.TRUE. C C W = L(L+)S = HS C CALL MVMLLD(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*MAX(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.D0 85 CONTINUE 90 CONTINUE C C FIND Q, (L+) SUCH THAT Q(L+) = (L+) + U(W+) C 93 CALL QRUPDD(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