C ALGORITHM 739, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 20, NO. 4, DECEMBER, 1994, P.518-530. c*** driver.f PROGRAM DRIVER IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION GPLS(50),H(50,4), Z TYPX(50),X(50),XPLS(50), Z WRK(50,8) INTEGER IWRK(4),GRDFLG,HESFLG EXTERNAL FCN,GRD,HSN NR=50 OPEN(5,FILE='wood.inp') CALL PRTFCN(N) READ (5,*)(X(I),I=1,N) PRINT *, "INITIAL APPROXIMATION TO THE SOLUTION X*:" PRINT *,(X(I),I=1,N) C SHORT FORM PRINT *," " PRINT *,"CALLING TENSRD - ALL INPUT PARAMETERS ARE SET" PRINT *," TO THEIR DEFAULT VALUES." PRINT *," OUTPUT WILL BE WRITTEN TO THE STANDARD OUTPUT." PRINT*," " CALL TENSRD(NR,N,X,FCN,MSG,XPLS,FPLS,GPLS,H,ITNNO,WRK,IWRK) PRINT *,"RESULTS OF TENSRD:" WRITE (*,211)(GPLS(I),I=1,N) PRINT *,'XPLS=',(XPLS(I),I=1,N) PRINT *,'FPLS=',FPLS,' ITNNO=',ITNNO,' MSG=',MSG C LONG FORM CALL PRTFCN(N) REWIND(5) READ (5,*)(X(I),I=1,N) PRINT *, "INITIAL APPROXIMATION TO THE SOLUTION X*:" PRINT *,(X(I),I=1,N) PRINT *," " PRINT *,"CALLING TENSOR AFTER RESETTING INPUT PARAMETERS" PRINT *,"IPR, MSG AND METHOD." PRINT *,"THE STANDARD METHOD IS USED IN THIS EXAMPLE." PRINT *,"ADDITIONAL OUTPUT WILL BE WRITTEN TO FILE 'FORT8'." PRINT *," " 997 CALL DFAULT(N,TYPX,FSCALE,GRADTL,STEPTL,ITNLIM,STEPMX, Z IPR,METHOD,GRDFLG,HESFLG,NDIGIT,MSG) OPEN(8,FILE='FORT8') IPR=8 MSG=2 METHOD=0 CALL TENSOR(NR,N,X,FCN,GRD,HSN,TYPX,FSCALE,GRADTL,STEPTL, Z ITNLIM,STEPMX,IPR,METHOD,GRDFLG,HESFLG, Z NDIGIT,MSG,XPLS,FPLS,GPLS,H,ITNNO,WRK,IWRK) PRINT *,"RESULTS OF TENSOR:" WRITE (*,211)(GPLS(I),I=1,N) 211 FORMAT(' G=',999E20.13) PRINT *,'XPLS=',(XPLS(I),I=1,N) PRINT *,'FPLS=',FPLS,' ITNNO=',ITNNO,' MSG=',MSG ENDFILE(8) CLOSE(8) CLOSE(5) STOP END subroutine fcn(n,x,f) implicit double precision (a-h,o-z) dimension x(n) if ( n .le. 0 )then n=4 else f=100d0*(x(1)*x(1)-x(2))**2d0+(1d0-x(1))**2d0 Z +90d0*(x(3)*x(3)-x(4))**2d0 Z +(1d0-x(3))**2d0+10.1d0*((1d0-x(2))**2d0+(1d0-x(4))**2d0) Z +19.8d0*(1d0-x(2))*(1d0-x(4)) endif return end subroutine prtfcn(n) n=4 print *,'__________________________________________' print *,'f(x)= Wood function' print *,'__________________________________________' return end C ---------------------- C | G R D | C ---------------------- SUBROUTINE GRD(N,X,G) C DUMMY ROUTINE END C ---------------------- C | H S N | C ---------------------- SUBROUTINE HSN(NR,N,X,H) C DUMMY ROUTINE END c*** tensrd.f C ---------------------- C | T E N S O R | C ---------------------- SUBROUTINE TENSOR(NR,N,X,FCN,GRD,HSN,TYPX,FSCALE,GRADTL,STEPTL, Z ITNLIM,STEPMX,IPR,METHOD,GRDFLG,HESFLG,NDIGIT, Z MSG,XPLS,FPLS,GPLS,H,ITNNO,WRK,IWRK) C C AUTHORS: TA-TUNG CHOW, ELIZABETH ESKOW AND ROBERT B. SCHNABEL C C DATE: MARCH 29, 1991 C C PURPOSE: C DERIVATIVE TENSOR METHOD FOR UNCONSTRAINED OPTIMIZATION C THE METHOD BASES EACH ITERATION ON A SPECIALLY CONSTRUCTED C FOURTH ORDER MODEL OF THE OBJECTIVE FUNCTION. THE MODEL C INTERPOLATES THE FUNCTION VALUE AND GRADIENT FROM THE PREVIOUS C ITERATE AND THE CURRENT FUNCTION VALUE, GRADIENT AND HESSIAN C MATRIX. C C BLAS SUBROUTINES: DCOPY,DDOT,DSCAL C UNCMIN SUBROUTINES: DFAULT, OPTCHK, GRDCHK, HESCHK, LNSRCH, FSTOFD, C SNDOFD, BAKSLV, FORSLV, OPTSTP C MODCHL SUBROUTINES: MODCHL, INIT, GERCH, FIN2X2 C C----------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE C GRD --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT C HSN --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN C HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART C AND DIAGONAL OF H C TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN C GRADTL --> GRADIENT TOLERANCE C STEPTL --> STEP TOLERANCE C ITNLIM --> ITERATION LIMIT C STEPMX --> MAXIMUM STEP LENGTH ALLOWED C IPR --> OUTPUT UNIT NUMBER C METHOD --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT C EACH ITERATION C IF VALUE IS 1 THEN TRY BOTH TENSOR AND NEWTON C STEPS AT EACH ITERATION C GRDFLG --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED C HESFLG --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED C NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN C MSG --> OUTPUT MESSAGE CONTROL C XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) C FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) C GPLS(N) <-- CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) C H(N,N) --> HESSIAN C ITNNO <-- NUMBER OF ITERATIONS C WRK(N,8)--> WORK SPACE C IWRK(N) --> WORK SPACE C INTEGER NR,N,ITNLIM,IPR,METHOD,GRDFLG,HESFLG,NDIGIT,MSG INTEGER ITNNO,IWRK(N) DOUBLE PRECISION X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX DOUBLE PRECISION XPLS(N),FPLS,GPLS(N),H(NR,N),WRK(NR,8) EXTERNAL FCN,GRD,HSN C C EQUIVALENCE WRK(N,1) = G(N) C WRK(N,2) = S(N) C WRK(N,3) = D(N) C WRK(N,4) = DN(N) C WRK(N,6) = E(N) C WRK(N,7) = WK1(N) C WRK(N,8) = WK2(N) C CALL OPT(NR,N,X,FCN,GRD,HSN,TYPX,FSCALE,GRADTL,STEPTL, Z ITNLIM,STEPMX,IPR,METHOD,GRDFLG,HESFLG,NDIGIT, Z MSG,XPLS,FPLS,GPLS,H,ITNNO, Z WRK(1,1),WRK(1,2),WRK(1,3),WRK(1,4), Z WRK(1,6),WRK(1,7),WRK(1,8),IWRK) RETURN END C ---------------------- C | T E N S R D | C ---------------------- SUBROUTINE TENSRD(NR,N,X,FCN,MSG,XPLS,FPLS,GPLS,H,ITNNO, Z WRK,IWRK) C C PURPOSE: C SHORT CALLING SEQUENCE SUBROUTINE C C----------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE C MSG --> OUTPUT MESSAGE CONTROL C XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) C FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) C GPLS(N) <-- GRADIENT AT FINAL POINT (OUTPUT) C H(N,N) --> HESSIAN C ITNNO <-- NUMBER OF ITERATIONS C WRK(N,8) --> WORK SPACE C IWRK(N) --> WORK SPACE C INTEGER NR,N,MSG,ITNNO,IWRK(N) DOUBLE PRECISION X(N),XPLS(N),FPLS,GPLS(N),H(NR,N),WRK(NR,8) DOUBLE PRECISION FSCALE,GRADTL,STEPTL,STEPMX INTEGER GRDFLG,HESFLG,ITNLIM,IPR,METHOD,NDIGIT EXTERNAL FCN,GRD,HSN C C EQUIVALENCE WRK(N,1) = G(N) C WRK(N,2) = S(N) C WRK(N,3) = D(N) C WRK(N,4) = DN(N) C WRK(N,5) = TYPX(N) C WRK(N,6) = E(N) C WRK(N,7) = WK1(N) C WRK(N,8) = WK2(N) C CALL DFAULT(N,WRK(1,5),FSCALE,GRADTL,STEPTL,ITNLIM,STEPMX,IPR, Z METHOD,GRDFLG,HESFLG,NDIGIT,MSG) CALL TENSOR(NR,N,X,FCN,GRD,HSN,WRK(1,5),FSCALE,GRADTL,STEPTL, Z ITNLIM,STEPMX,IPR,METHOD,GRDFLG,HESFLG,NDIGIT, Z MSG,XPLS,FPLS,GPLS,H,ITNNO,WRK,IWRK) END C ---------------------- C | O P T | C ---------------------- SUBROUTINE OPT(NR,N,X,FCN,GRD,HSN,TYPX,FSCALE,GRADTL,STEPTL, Z ITNLIM,STEPMX,IPR,METHOD,GRDFLG,HESFLG,NDIGIT, Z MSG,XPLS,FPLS,GPLS,H,ITNNO,G,S,D,DN,E,WK1,WK2, Z PIVOT) C C PURPOSE: C DERIVATIVE TENSOR METHODS FOR UNCONSTRAINED OPTIMIZATION C C----------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE C FCN: R(N) --> R(1) C GRD --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT C FCN: R(N) --> R(N) C HSN --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN C FCN: R(N) --> R(N X N) C HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART C AND DIAGONAL OF H C TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN C GRADTL --> GRADIENT TOLERANCE C STEPTL --> STEP TOLERANCE C ITNLIM --> ITERATION LIMIT C STEPMX --> MAXIMUM STEP LENGTH ALLOWED C IPR --> OUTPUT UNIT NUMBER C METHOD --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT C EACH ITERATION C GRDFLG --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED C HESFLG --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED C NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN C MSG --> OUTPUT MESSAGE CONTRO C XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) C FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) C GPLS(N) <-- CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) C H(N,N) --> HESSIAN C ITNNO <-- NUMBER OF ITERATIONS C G(N) --> PREVIOUS GRADIENT C S --> STEP TO PREVIOUS ITERATE (FOR TENSOR MODEL) C D --> CURRENT STEP (TENSOR OR NEWTON) C DN --> NEWTON STEP C E(N) --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION C WK1(N) --> WORKSPACE C WK2(N) --> WORKSPACE C PIVOT(N)--> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION C INTEGER NR,N,ITNLIM,IPR,METHOD,GRDFLG,HESFLG,NDIGIT,MSG INTEGER ITNNO,PIVOT(N),ICSCMX,I,ITRMCD,IRETCD,IMSG,NFCNT DOUBLE PRECISION F,FP,GPLS(N),G(N),H(NR,N),S(N) DOUBLE PRECISION D(N),WK1(N),TYPX(N),FSCALE,GRADTL,STEPTL DOUBLE PRECISION X(N),STEPMX,FINIT DOUBLE PRECISION WK2(N),FPLS,XPLS(N),E(N) DOUBLE PRECISION DN(N),EPS,RNF,ANALTL,GNORM,TWONRM DOUBLE PRECISION ADDMAX,TEMP,ALPHA,BETA,GD,FN,DDOT DOUBLE PRECISION RGX,RSX LOGICAL NOMIN,MXTAKE EXTERNAL FCN,GRD,HSN,DCOPY,DDOT C-------------------------------- C INITIALIZATION | C-------------------------------- C ITNNO=0 ICSCMX=0 NFCNT=0 CALL MCHEPS(EPS) CALL OPTCHK(NR,N,MSG,X,TYPX,FSCALE,GRADTL,STEPTL,ITNLIM, + NDIGIT,EPS,METHOD,GRDFLG,HESFLG,STEPMX,IPR) IF (MSG.LT.0) RETURN C C SCALE X DO 10 I=1,N X(I)=X(I)/TYPX(I) 10 CONTINUE C C-------------------------------- C INITIAL ITERATION | C-------------------------------- C C COMPUTE TYPICAL SIZE OF X DO 15 I=1,N DN(I)=1D0/TYPX(I) 15 CONTINUE RNF=MAX(10.0D0**(-NDIGIT),EPS) ANALTL=MAX(1D-2,SQRT(RNF)) C UNSCALE X AND COMPUTE F AND G DO 20 I=1,N WK1(I)=X(I)*TYPX(I) 20 CONTINUE CALL FCN(N,WK1,F) NFCNT=NFCNT+1 IF(GRDFLG .GE. 1)THEN CALL GRD(N,WK1,G) IF(GRDFLG .EQ. 1)THEN FSCALE=1D0 CALL GRDCHK(N,WK1,FCN,F,G,DN,TYPX,FSCALE,RNF,ANALTL,WK2, Z MSG,IPR,NFCNT) END IF ELSE CALL FSTOFD(1,1,N,WK1,FCN,F,G,TYPX,RNF,WK2,1,NFCNT) END IF IF(MSG .LT. -20) THEN RETURN END IF C GNORM=TWONRM(N,G) C C PRINT OUT 1ST ITERATION? IF(MSG .GE. 1)THEN WRITE(IPR,25)F 25 FORMAT(' INITIAL FUNCTION VALUE F=',E20.13) WRITE(IPR,30)(G(I),I=1,N) 30 FORMAT(' INITIAL GRADIENT G=',999E20.13) END IF C C TEST WHETHER THE INITIAL GUESS SATISFIES THE STOPPING CRITERIA IF(GNORM .LE. GRADTL)THEN FPLS=F DO 40,I=1,N XPLS(I)=X(I) GPLS(I)=G(I) 40 CONTINUE CALL OPTSTP(N,XPLS,FPLS,GPLS,X,ITNNO,ICSCMX, Z ITRMCD,GRADTL,STEPTL,FSCALE,ITNLIM,IRETCD, Z MXTAKE,IPR,MSG,RGX,RSX) GO TO 350 END IF FINIT=F C C------------------------ C ITERATION 1 | C------------------------ C ITNNO=ITNNO+1 C C COMPUTE HESSIAN IF(HESFLG .EQ. 0)THEN IF (GRDFLG .EQ. 1) THEN CALL FSTOFD(NR,N,N,WK1,GRD,G,H,TYPX,RNF,WK2,3,NFCNT) ELSE CALL SNDOFD(NR,N,WK1,FCN,F,H,TYPX,RNF,WK2,D,NFCNT) END IF ELSE IF(HESFLG .EQ. 2)THEN CALL HSN(NR,N,WK1,H) ELSE IF(HESFLG .EQ. 1)THEN C IN HESCHK GPLS,XPLS AND E ARE USED AS WORK SPACE CALL HESCHK(NR,N,WK1,FCN,GRD,HSN,F,G,H,DN,TYPX,RNF, Z ANALTL,GRDFLG,GPLS,XPLS,E,MSG,IPR,NFCNT) END IF END IF END IF IF(MSG .LT. -20)RETURN DO 50 I=2,N CALL DCOPY(I-1,H(I,1),NR,H(1,I),1) 50 CONTINUE C C CHOLESKY DECOMPOSITION FOR H (H=LLT) CALL CHOLDR(NR,N,H,D,EPS,PIVOT,E,WK1,ADDMAX) C C SOLVE FOR NEWTON STEP D DO 60 I=1,N 60 WK2(I)=-G(I) CALL FORSLV(NR,N,H,WK1,WK2) CALL BAKSLV(NR,N,H,D,WK1) C C APPLY LINESEARCH TO THE NEWTON STEP CALL LNSRCH(NR,N,X,F,G,D,XPLS,FPLS,MXTAKE,IRETCD,STEPMX, Z STEPTL,TYPX,FCN,WK1,NFCNT) C C UPDATE G C CALL DCOPY(N,GPLS(1),1,GP(1),1) C C UNSCALE XPLS AND COMPUTE GPLS DO 80 I=1,N WK1(I)=XPLS(I)*TYPX(I) 80 CONTINUE IF(GRDFLG .GE. 1)THEN CALL GRD(N,WK1,GPLS) ELSE CALL FSTOFD(1,1,N,WK1,FCN,FPLS,GPLS,TYPX,RNF,WK2,1,NFCNT) END IF C C CHECK STOPPING CONDITIONS CALL OPTSTP(N,XPLS,FPLS,GPLS,X,ITNNO,ICSCMX, Z ITRMCD,GRADTL,STEPTL,FSCALE,ITNLIM,IRETCD, Z MXTAKE,IPR,MSG,RGX,RSX) C C IF ITRMCD > 0 THEN STOPPING CONDITIONS SATISFIED IF(ITRMCD .GT. 0)GO TO 350 C C UPDATE X,F AND S FOR TENSOR MODEL FP=F F=FPLS DO 90 I=1,N TEMP=XPLS(I) S(I)=X(I)-TEMP X(I)=TEMP 90 CONTINUE C C IF MSG >= 2 THEN PRINT OUT EACH ITERATION IF(MSG .GE. 2)THEN WRITE (IPR,103)ITNNO 103 FORMAT(' ----------- ITERATION ',I4,' ----------------') WRITE(IPR,104)(X(I),I=1,N) 104 FORMAT(' X=',999E20.13) WRITE(IPR,105)FPLS 105 FORMAT(' F(X)=',E20.13) WRITE(IPR,106)(GPLS(I),I=1,N) 106 FORMAT(' G(X)=',999E20.13) END IF IF (MSG .GE. 3 )THEN WRITE (IPR,108) NFCNT,RGX 108 FORMAT('FUNCTION EVAL COUNT:',I6,' REL. GRAD. MAX:',E20.13) WRITE (IPR,110) RSX 110 FORMAT('REL. STEP MAX :',E20.13) ENDIF C C------------------------ C ITERATION > 1 | C------------------------ C C UNSCALE X AND COMPUTE H 200 DO 210 I=1,N WK1(I)=X(I)*TYPX(I) 210 CONTINUE IF(HESFLG .EQ. 0)THEN IF (GRDFLG .EQ. 1) THEN CALL FSTOFD(NR,N,N,WK1,GRD,G,H,TYPX,RNF,WK2,3,NFCNT) ELSE CALL SNDOFD(NR,N,WK1,FCN,F,H,TYPX,RNF,WK2,D,NFCNT) END IF ELSE CALL HSN(NR,N,WK1,H) END IF DO 230 I=2,N CALL DCOPY(I-1,H(I,1),NR,H(1,I),1) 230 CONTINUE C C IF METHOD = 0 THEN USE NEWTON STEP ONLY IF (METHOD .EQ. 0) THEN C C CHOLESKY DECOMPOSITION FOR H CALL CHOLDR(NR,N,H,WK2,EPS,PIVOT,E,WK1,ADDMAX) C C COMPUTE NEWTON STEP DO 240 I=1,N WK1(I)=-GPLS(I) 240 CONTINUE CALL FORSLV(NR,N,H,WK2,WK1) CALL BAKSLV(NR,N,H,D,WK2) C C NO TENSOR STEP NOMIN=.TRUE. GO TO 300 C END IF C C FORM TENSOR MODEL CALL MKMDL(NR,N,F,FP,GPLS,G,S,H,ALPHA,BETA,WK1,D) C C SOLVE TENSOR MODEL AND COMPUTE NEWTON STEP C ON INPUT : SH IS STORED IN WK1 C A=(G-GPLS-SH-S*BETA/(6*STS)) IS STORED IN D C ON OUTPUT: NEWTON STEP IS STORED IN DN C TENSOR STEP IS STORED IN D CALL SLVMDL(NR,N,H,XPLS,WK2,E,G,S,GPLS,PIVOT,D,WK1,DN, Z ALPHA,BETA,NOMIN,EPS) C C IF TENSOR MODEL HAS NO MINIMIZER THEN USE NEWTON STEP IF(NOMIN)THEN CALL DCOPY(N,DN(1),1,D(1),1) GO TO 300 END IF C C IF TENSOR STEP IS NOT IN DESCENT DIRECTION THEN USE NEWTON STEP GD = DDOT(N,GPLS(1),1,D(1),1) IF(GD .GT. 0D0)THEN CALL DCOPY(N,DN(1),1,D(1),1) NOMIN=.TRUE. END IF C 300 ITNNO=ITNNO+1 CALL DCOPY(N,GPLS(1),1,G(1),1) C C APPLY LINESEARCH TO TENSOR (OR NEWTON) STEP CALL LNSRCH(NR,N,X,F,G,D,XPLS,FPLS,MXTAKE,IRETCD,STEPMX, Z STEPTL,TYPX,FCN,WK1,NFCNT) C IF(.NOT. NOMIN)THEN C TENSOR STEP IS FOUND AND IN DESCENT DIRECTION, C APPLY LINESEARCH TO NEWTON STEP C NEW NEWTON POINT IN WK2 CALL LNSRCH(NR,N,X,F,G,DN,WK2,FN,MXTAKE,IRETCD,STEPMX, Z STEPTL,TYPX,FCN,WK1,NFCNT) C C COMPARE TENSOR STEP TO NEWTON STEP C IF NEWTON STEP IS BETTER, SET NEXT ITERATE TO NEW NEWTON POINT IF(FN .LT. FPLS)THEN FPLS=FN CALL DCOPY(N,DN(1),1,D(1),1) CALL DCOPY(N,WK2(1),1,XPLS(1),1) END IF END IF DO 320 I=1,N D(I)=XPLS(I)-X(I) 320 CONTINUE C C UNSCALE XPLS, AND COMPUTE FPLS AND GPLS DO 330 I=1,N WK1(I)=XPLS(I)*TYPX(I) 330 CONTINUE CALL FCN(N,WK1,FPLS) NFCNT=NFCNT+1 IF(GRDFLG .GE. 1)THEN CALL GRD(N,WK1,GPLS) ELSE CALL FSTOFD(1,1,N,WK1,FCN,FPLS,GPLS,TYPX,RNF,WK2,1,NFCNT) END IF C C CHECK STOPPING CONDITIONS IMSG=MSG CALL OPTSTP(N,XPLS,FPLS,GPLS,X,ITNNO,ICSCMX, Z ITRMCD,GRADTL,STEPTL,FSCALE,ITNLIM,IRETCD, Z MXTAKE,IPR,MSG,RGX,RSX) C C IF ITRMCD = 0 THEN NOT OVER YET IF(ITRMCD .EQ. 0)GO TO 500 C C IF MSG >= 1 THEN PRINT OUT FINAL ITERATION 350 IF(IMSG .GE. 1)THEN C C TRANSFORM X BACK TO ORIGINAL SPACE DO 360 I=1,N XPLS(I)=XPLS(I)*TYPX(I) 360 CONTINUE WRITE(IPR,370)(XPLS(I),I=1,N) 370 FORMAT(' FINAL X=',999E20.13) WRITE(IPR,380)(GPLS(I),I=1,N) 380 FORMAT(' GRADIENT G=',999E20.13) WRITE(IPR,390)FPLS,ITNNO 390 FORMAT(' FUNCTION VALUE F(X)=',E20.13, Z ' AT ITERATION ',I4) IF (IMSG .GE. 3) THEN WRITE (IPR,400) NFCNT,RGX 400 FORMAT('FUNCTION EVAL COUNT:',I6,' REL. GRAD. MAX:',E20.13) WRITE (IPR,410) RSX 410 FORMAT('REL. STEP MAX :',E20.13) ENDIF C UPDATE THE HESSIAN IF(HESFLG .EQ. 0)THEN IF (GRDFLG .EQ. 1) THEN CALL FSTOFD(NR,N,N,XPLS,GRD,GPLS,H,TYPX,RNF,WK2,3,NFCNT) ELSE CALL SNDOFD(NR,N,XPLS,FCN,FPLS,H,TYPX,RNF,WK2,D,NFCNT) END IF ELSE CALL HSN(NR,N,XPLS,H) END IF END IF RETURN C C UPDATE INFORMATION AT THE CURRENT POINT 500 CALL DCOPY(N,XPLS(1),1,X(1),1) DO 550 I=1,N S(I)=-D(I) 550 CONTINUE C C IF TOO MANY ITERATIONS THEN RETURN IF(ITNNO .GT. ITNLIM)GO TO 350 C C IF MSG >= 2 THEN PRINT OUT EACH ITERATION IF(MSG .GE. 2)THEN WRITE (IPR,560)ITNNO 560 FORMAT(' ----------- ITERATION ',I4,' ----------------') WRITE(IPR,570)(X(I),I=1,N) 570 FORMAT(' X=',999E20.13) WRITE(IPR,580)FPLS 580 FORMAT(' F(X)=',E20.13) WRITE(IPR,590)(GPLS(I),I=1,N) 590 FORMAT(' G(X)=',999E20.13) END IF IF (MSG .GE. 3) THEN WRITE (IPR,600) NFCNT,RGX 600 FORMAT('FUNCTION EVAL COUNT:',I6,' REL. GRAD. MAX:',E20.13) WRITE (IPR,610) RSX 610 FORMAT('REL. STEP MAX :',E20.13) ENDIF C UPDATE F FP=F F=FPLS C C PERFORM NEXT ITERATION GO TO 200 C C END OF ITERATION > 1 C END C ---------------------- C | D F A U L T | C ---------------------- SUBROUTINE DFAULT(N,TYPX,FSCALE,GRADTL,STEPTL,ITNLIM, Z STEPMX,IPR,METHOD,GRDFLG,HESFLG,NDIGIT,MSG) INTEGER N,ITNLIM,IPR,METHOD,NDIGIT,MSG,GRDFLG,HESFLG,I DOUBLE PRECISION TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX,EPS,TEMP CALL MCHEPS(EPS) METHOD=1 FSCALE=1D0 GRDFLG=0 HESFLG=0 DO 1 I=1,N TYPX(I)=1D0 1 CONTINUE TEMP=EPS**(1D0/3D0) GRADTL=TEMP STEPTL=TEMP*TEMP NDIGIT=-LOG10(EPS) C SET ACTUAL DFAULT VALUE OF STEPMX IN OPTCHK STEPMX=0D0 ITNLIM=100 IPR=6 MSG=1 END C ---------------------- C | O P T C H K | C ---------------------- SUBROUTINE OPTCHK(NR,N,MSG,X,TYPX,FSCALE,GRADTL,STEPTL, + ITNLIM,NDIGIT,EPS,METHOD,GRDFLG,HESFLG,STEPMX,IPR) C C PURPOSE C ------- C CHECK INPUT FOR REASONABLENESS C C PARAMETERS C ---------- C NR <--> ROW DIMENSION OF H AND WRK C N --> DIMENSION OF PROBLEM C X(N) --> ON ENTRY, ESTIMATE TO ROOT OF FCN C TYPX(N) <--> TYPICAL SIZE OF EACH COMPONENT OF 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 STEPTL <--> TOLERANCE AT WHICH STEP LENGTH 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 EPS --> MACHINE PRECISION C METHOD <--> ALGORITHM INDICATOR C GRDFLG <--> =1 IF ANALYTIC GRADIENT SUPPLIED C HESFLG <--> =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 INTEGER NR,N,ITNLIM,IPR,METHOD,NDIGIT,MSG,GRDFLG,HESFLG,I DOUBLE PRECISION X(N),TYPX(N),FSCALE,GRADTL,STEPTL,STEPMX,EPS DOUBLE PRECISION STPSIZ,TEMP INTRINSIC LOG10, MAX, SQRT C C CHECK DIMENSION OF PROBLEM C IF(N.LE.0) GO TO 805 IF(NR.LT.N) GO TO 806 C C CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES. C IF NOT, SET THEM TO DEFAULT VALUES. IF(METHOD.NE.0) METHOD=1 IF(GRDFLG.NE.2 .AND. GRDFLG.NE.1) GRDFLG=0 IF(HESFLG.NE.2 .AND. HESFLG.NE.1) HESFLG=0 IF(MSG.GT.3 .OR. MSG.LT.0 .) MSG=1 C C COMPUTE SCALE MATRIX C DO 10 I=1,N IF(TYPX(I).EQ.0.) TYPX(I)=1D0 IF(TYPX(I).LT.0.) TYPX(I)=-TYPX(I) 10 CONTINUE C C CHECK MAXIMUM STEP SIZE C IF (STEPMX .GT. 0D0) GO TO 20 STPSIZ = 0D0 DO 15 I = 1, N STPSIZ = STPSIZ + X(I)*X(I)/TYPX(I)*TYPX(I) 15 CONTINUE STPSIZ = SQRT(STPSIZ) STEPMX = MAX(1D3*STPSIZ, 1D3) 20 CONTINUE C CHECK FUNCTION SCALE IF(FSCALE.EQ.0.) FSCALE=1D0 IF(FSCALE.LT.0.) FSCALE=-FSCALE C CHECK GRADIENT TOLERANCE IF(GRADTL.LT.0.) GRADTL=EPS**(1D0/3D0) C CHECK STEP TOLERANCE IF(STEPTL.LT.0.) THEN TEMP=EPS**(1D0/3D0) STEPTL=TEMP*TEMP END IF C C CHECK ITERATION LIMIT IF(ITNLIM.LE.0) ITNLIM=100 C C CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN IF(NDIGIT.LE.0) NDIGIT=-LOG10(EPS) IF(10D0**(-NDIGIT).LE.EPS) NDIGIT=-LOG10(EPS) RETURN C C ERROR EXITS C 805 WRITE(IPR,901) N MSG=-20 RETURN 806 WRITE(IPR,902) NR,N MSG=-20 RETURN 901 FORMAT(32H0OPTCHK ILLEGAL DIMENSION, N=,I5) 902 FORMAT('OPTCHK ILLEGAL INPUT VALUES OF: NR=',I5,', N=',I5, + ', NR MUST BE <= N.') END C ---------------------- C | G R D C H K | C ---------------------- SUBROUTINE GRDCHK(N,X,FCN,F,G,TYPSIZ,TYPX,FSCALE,RNF, + ANALTL,WRK1,MSG,IPR,IFN) 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 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 IFN <--> NUMBER OF FUNCTION EVALUATIONS INTEGER N,MSG,IPR,IFN,KER,I DOUBLE PRECISION F,X(N),G(N),TYPX(N),FSCALE,ANALTL,RNF DOUBLE PRECISION TYPSIZ(N),GS,WRK DOUBLE PRECISION WRK1(N) EXTERNAL FCN INTRINSIC ABS,MAX C C COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO C ANALYTIC GRADIENT. C CALL FSTOFD(1,1,N,X,FCN,F,WRK1,TYPX,RNF,WRK,1,IFN) KER=0 DO 5 I=1,N GS=MAX(ABS(F),FSCALE)/MAX(ABS(X(I)),TYPSIZ(I)) IF(ABS(G(I)-WRK1(I)).GT.MAX(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 C ---------------------- C | H E S C H K | C ---------------------- SUBROUTINE HESCHK(NR,N,X,FCN,GRD,HSN,F,G,A,TYPSIZ,TYPX,RNF, + ANALTL,IAGFLG,UDIAG,WRK1,WRK2,MSG,IPR,IFN) C C PURPOSE C ------- C CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN C (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN C HSN 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 GRD --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN. C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE C HSN --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN. C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE C HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART C AND DIAGONAL OF A 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 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 IFN <--> NUMBER OF FUNCTION EVALUTATIONS INTEGER NR,N,MSG,IPR,KER,I,J,IAGFLG,JP1,IM1,IFN DOUBLE PRECISION F,X(N),G(N),TYPX(N),ANALTL,RNF DOUBLE PRECISION TYPSIZ(N) DOUBLE PRECISION WRK1(N) DOUBLE PRECISION A(NR,1),UDIAG(N),WRK2(N),HS EXTERNAL FCN,GRD,HSN INTRINSIC ABS,MAX C C COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN. C IF(IAGFLG.EQ.1) CALL FSTOFD(NR,N,N,X,GRD,G,A,TYPX,RNF, Z WRK1,3,IFN) IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,X,FCN,F,A,TYPX,RNF,WRK1, Z WRK2,IFN) 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 DO 32, I=1,N DO 32, J=I,N A(J,I)=0D0 32 CONTINUE C CALL HSN(NR,N,X,A) DO 40 J=1,N HS=MAX(ABS(G(J)),1.0D0)/MAX(ABS(X(J)),TYPSIZ(J)) IF(ABS(A(J,J)-UDIAG(J)).GT.MAX(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.MAX(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 C ----------------------- C | M C H E P S | C ----------------------- SUBROUTINE MCHEPS(EPS) C C PURPOSE: C COMPUTE MACHINE PRECISION C C------------------------------------------------------------------------- C C PARAMETERS: C C EPS <-- MACHINE PRECISION C DOUBLE PRECISION EPS,TEMP,TEMP1 TEMP = 1.D0 20 CONTINUE TEMP = TEMP / 2.D0 TEMP1 = TEMP + 1.D0 IF (TEMP1 .NE. 1.D0) GOTO 20 EPS = TEMP * 2.D0 RETURN END C FUNCTION TWONRM(N,V) C C PURPOSE: C COMPUTER L-2 NORM C C-------------------------------------------------------------------------- C C PARAMETER: C C N --> DIMENSION OF PROBLEM C V(N) --> VECTOR WHICH L-2 NORM IS EVALUATED C INTEGER N DOUBLE PRECISION V(N),TEMP,TWONRM,DDOT EXTERNAL DDOT INTRINSIC SQRT TEMP = DDOT(N,V(1),1,V(1),1) TWONRM=SQRT(TEMP) RETURN END C ---------------------- C | L N S R C H | C ---------------------- SUBROUTINE LNSRCH(NR,N,X,F,G,P,XPLS,FPLS,MXTAKE, + IRETCD,STEPMX,STEPTL,TYPX,FCN,W2,NFCNT) C C THE ALPHA CONDITION ONLY LINE SEARCH C 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 TYPX(N) --> DIAGONAL SCALING MATRIX FOR X (NOT IN UNCMIN) C IPR --> DEVICE TO WHICH TO SEND OUTPUT C W2 --> WORKING SPACE C NFCNT <--> NUMBER OF FUNCTION EVALUTIONS C C INTERNAL VARIABLES C ------------------ C SLN NEWTON LENGTH C RLN RELATIVE LENGTH OF NEWTON STEP C INTEGER NR,N,IRETCD,I,K,NFCNT DOUBLE PRECISION STEPMX,STEPTL,TYPX(N),W2(N),X(N),G(N) DOUBLE PRECISION XPLS(N),F,FPLS ,P(N) DOUBLE PRECISION ALPHA,TMP,SLN,SCL,SLP,DDOT DOUBLE PRECISION RLN,TEMP,TEMP1,TEMP2,ALMBDA,T1,T2,T3 DOUBLE PRECISION ALMBMN,TLMBDA,PLMBDA,PFPLS,DISC,A,B LOGICAL MXTAKE EXTERNAL FCN,DSCAL,DDOT INTRINSIC ABS,MAX,SQRT C MXTAKE=.FALSE. IRETCD=2 ALPHA=1D-4 C$ WRITE(IPR,954) C$ WRITE(IPR,955) (P(I),I=1,N) TMP=0D0 DO 5 I=1,N TMP=TMP+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 DSCAL(N,SCL,P(1),1) SLN=STEPMX C$ WRITE(IPR,954) C$ WRITE(IPR,955) (P(I),I=1,N) 10 CONTINUE SLP=DDOT(N,G(1),1,P(1),1) RLN=0D0 DO 15 I=1,N TEMP=1D0 TEMP1=ABS(X(I)) TEMP2=MAX(TEMP1,TEMP) TEMP1=ABS(P(I)) RLN=MAX(RLN,TEMP1/TEMP2) 15 CONTINUE ALMBMN=STEPTL/RLN ALMBDA=1.0D0 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) THEN RETURN END IF DO 105 I=1,N XPLS(I)=X(I) + ALMBDA*P(I) 105 CONTINUE DO 101 K=1,N W2(K)=XPLS(K)*TYPX(K) 101 CONTINUE CALL FCN(N,W2,FPLS) NFCNT=NFCNT+1 C$ WRITE(IPR,956) ALMBDA C$ WRITE(IPR,951) C$ WRITE(IPR,955) (XPLS(I),I=1,N) C$ WRITE(IPR,953) FPLS IF(FPLS.GT. F+SLP*ALPHA*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.1D0 .AND. SLN.GT. .99D0*STEPMX) MXTAKE=.TRUE. GO TO 100 C C SOLUTION NOT (YET) FOUND C C ELSE 130 IF(ALMBDA .GE. ALMBMN) GO TO 140 C IF(ALMBDA .LT. ALMBMN) 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.1D0) GO TO 150 C IF(ALMBDA.EQ.1.0) C THEN C C FIRST BACKTRACK: QUADRATIC FIT C TLMBDA=-SLP/(2D0*(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=1D0/(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.0D0,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.0D0,A)*SQRT(DISC))/(3.0*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*.1 GO TO 190 C ELSE 180 ALMBDA=TLMBDA C ENDIF C ENDIF C ENDIF 190 GO TO 100 956 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 ALMBMN=,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 C ---------------------- C | Z H Z | C ---------------------- SUBROUTINE ZHZ(NR,N,Y,H,U,T) C C PURPOSE: C COMPUTE QTHQ(N,N) AND ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND C FIRST N-1 COLUMNS OF QTHQ C C--------------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C Y(N) --> FIRST BASIS IN Q C H(N,N) <--> ON INPUT : HESSIAN C ON OUTPUT: QTHQ (ZTHZ) C U(N) <-- VECTOR TO FORM Q AND Z C T(N) --> WORKSPACE C INTEGER NR,N,I,J DOUBLE PRECISION T(N),Y(N),H(NR,N),U(N),S,SGN,YNORM DOUBLE PRECISION DDOT,D,TEMP1,TEMP2 EXTERNAL DDOT,DCOPY INTRINSIC SQRT,ABS C C U=Y+SGN(Y(N))||Y||E(N) IF( Y(N) .NE. 0D0)THEN SGN=Y(N)/ABS(Y(N)) ELSE SGN=1D0 END IF YNORM=DDOT(N,Y(1),1,Y(1),1) YNORM=SQRT(YNORM) U(N)=Y(N)+SGN*YNORM CALL DCOPY(N-1,Y(1),1,U(1),1) C C D=UTU/2 D=DDOT(N,U(1),1,U(1),1) D=D/2D0 C C T=2HU/UTU DO 40 I=1,N T(I)=0D0 DO 30 J=1,N T(I)=H(I,J)*U(J)+T(I) 30 CONTINUE T(I)=T(I)/D 40 CONTINUE C C S=4UHU/(UTU)**2 S = DDOT(N,U(1),1,T(1),1) S=S/D C C COMPUTE QTHQ (ZTHZ) DO 70 I=1,N TEMP1=U(I) DO 60 J=1,N TEMP2=U(J) H(I,J)=H(I,J)-U(I)*T(J)-T(I)*U(J)+U(I)*U(J)*S 60 CONTINUE 70 CONTINUE RETURN END C ---------------------- C | S O L V E W | C ---------------------- SUBROUTINE SOLVEW(NR,N,AL,U,W,B) C C PURPOSE: C SOLVE L*W=ZT*V C C---------------------------------------------------------------------- C C PARAMETERS: C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C AL(N-1,N-1) --> LOWER TRIAGULAR MATRIX C U(N) --> VECTOR TO FORM Z C W(N) --> ON INPUT : VECTOR V IN SYSTEM OF LINEAR EQUATIONS C ON OUTPUT: SOLUTION OF SYSTEM OF LINEAR EQUATIONS C B(N) --> WORKSPACE TO STORE ZT*V INTEGER NR,N,I,J DOUBLE PRECISION B(N),AL(NR,1),W(N),U(N),D,DDOT EXTERNAL DDOT C C FORM ZT*V (STORED IN B) D = DDOT(N,U(1),1,U(1),1) D=D/2D0 DO 20 I=1,N-1 B(I)=0D0 DO 15 J=1,N B(I)=B(I)+U(J)*U(I)*W(J)/D 15 CONTINUE B(I)=W(I)-B(I) 20 CONTINUE C C SOLVE LW=ZT*V CALL FORSLV(NR,N-1,AL,W,B) RETURN END C ---------------------- C | D S T A R | C ---------------------- SUBROUTINE DSTAR(NR,N,U,S,W1,W2,W3,SIGMA,AL,D) C C PURPOSE: C COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) C C------------------------------------------------------------------------ C C PARAMETERS: C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C U(N) --> VECTOR TO FORM Z C S(N) --> PREVIOUS STEP C W1(N) --> L**-1*ZT*A, WHERE A IS DESCRIBED IN SUBROUTINE SLVMDL C W2(N) --> L**-1*ZT*SH, WHERE H IS CURRENT HESSIAN C W3(N) --> L**-1*ZT*G, WHERE G IS CURRENT GRADIENT C SIGMA --> SOLUTION FOR REDUCED ONE VARIABLE MODEL C AL(N-1,N-1) --> LOWER TRIANGULAR MATRIX L C D(N) --> TENSOR STEP C INTEGER NR,N,I DOUBLE PRECISION U(N),S(N),W1(N-1),W2(N-1),W3(N-1),AL(NR,1),D(N) DOUBLE PRECISION SIGMA,DDOT,UTU,TEMP,UTT EXTERNAL DDOT IF (N.EQ.1) THEN D(1)=SIGMA * S(1) ELSE C C COMPUTE T(SIGMA)=-(ZTHZ)*ZT*(G+SIGMA*SH+SIGMA**2*A/2) (STORED IN D) DO 10 I=1,N-1 W2(I)=W3(I)+SIGMA*W2(I)+0.5D0*W1(I)*SIGMA*SIGMA 10 CONTINUE CALL BAKSLV(NR,N-1,AL,D,W2) D(N)=0D0 C C COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) UTU=DDOT(N,U(1),1,U(1),1) UTT=DDOT(N,U(1),1,D(1),1) TEMP=UTT/UTU DO 50 I=1,N D(I)=SIGMA*S(I)-(D(I)-2D0*U(I)*TEMP) 50 CONTINUE END IF RETURN END C ---------------------- C | M K M D L | C ---------------------- SUBROUTINE MKMDL(NR,N,F,FP,G,GP,S,H,ALPHA,BETA,SH,A) C C PURPOSE: C FORM TENSOR MODEL C C----------------------------------------------------------------------- C C PARAMETERS: C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C F --> CURRENT FUNCTION VALUE C FP --> PREVIOUS FUNCTION VALUE C G(N) --> CURRENT GRADIENT C GP(N) --> PREVIOUS GRADIENT C S(N) --> STEP TO PREVIOUS POINT C H(N,N) --> HESSIAN C ALPHA <-- SCALAR TO FORM 3RD ORDER TERM OF TENSOR MODEL C BETA <-- SCALAR TO FORM 4TH ORDER TERM OF TENSOR MODEL C SH(N) <-- SH C A(N) <-- A=2*(GP-G-SH-S*BETA/(6*STS)) C INTEGER NR,N,I,J DOUBLE PRECISION F,FP,G(N),GP(N),S(N),H(NR,N),SH(N),A(N) DOUBLE PRECISION ALPHA,BETA,DDOT,GS,GPS,SHS,B1,B2,STS EXTERNAL DDOT C C COMPUTE SH DO 20 I=1,N SH(I)=0D0 DO 10 J=1,N SH(I)=SH(I)+S(J)*H(J,I) 10 CONTINUE 20 CONTINUE GS=DDOT(N,G(1),1,S(1),1) GPS=DDOT(N,GP(1),1,S(1),1) SHS=DDOT(N,SH(1),1,S(1),1) B1=GPS-GS-SHS B2=FP-F-GS-.5D0*SHS ALPHA=24D0*B2-6D0*B1 BETA=24D0*B1-72D0*B2 C C COMPUTE A STS=DDOT(N,S(1),1,S(1),1) DO 50 I=1,N A(I)=2D0*(GP(I)-G(I)-SH(I)-S(I)*BETA/(6D0*STS)) 50 CONTINUE RETURN END C ---------------------- C | S I G M A | C ---------------------- SUBROUTINE SIGMA(SGSTAR,A,B,C,D) C C PURPOSE: C COMPUTE DESIRABLE ROOT OF REDUCED ONE VARIABLE EQUATION C C------------------------------------------------------------------------- C C PARAMETERS: C SGSTAR --> DESIRABLE ROOT C A --> COEFFICIENT OF 3RD ORDER TERM C B --> COEFFICIENT OF 2ND ORDER TERM C C --> COEFFICIENT OF 1ST ORDER TERM C D --> COEFFICIENT OF CONSTANT TERM C DOUBLE PRECISION SGSTAR,A,B,C,D DOUBLE PRECISION S1,S2,S3 C C COMPUTE ALL THREE ROOTS CALL ROOTS(S1,S2,S3,A,B,C,D) C C SORT ROOTS CALL SORTRT(S1,S2,S3) C C CHOOSE DESIRABLE ROOT IF( A .GT. 0D0 )THEN SGSTAR=S3 IF( S2 .GE. 0D0 )THEN SGSTAR=S1 END IF ELSE C IF A < 0 THEN SGSTAR=S2 IF( (S1.GT.0D0) .OR. (S3.LT.0D0) )THEN IF( S1 .GT. 0D0)THEN SGSTAR=S1 ELSE SGSTAR=S3 END IF A=0D0 END IF END IF RETURN END C ---------------------- C | R O O T S | C ---------------------- SUBROUTINE ROOTS(S1,S2,S3,A,B,C,D) C C PURPOSE: C COMPUTE ROOT(S) OF 3RD ORDER EQUATION C C--------------------------------------------------------------------------- C C PARAMETERS: C S1 <-- ROOT (IF THREE ROOTS ARE C S2 <-- ROOT EQUAL, THEN S1=S2=S3) C S3 <-- ROOT C A --> COEFFICIENT OF 3RD ORDER TERM C B --> COEFFICIENT OF 2ND ORDER TERM C C --> COEFFICIENT OF 1ST ORDER TERM C D --> COEFFICIENT OF CONSTANT TERM DOUBLE PRECISION S1,S2,S3,A,B,C,D DOUBLE PRECISION PI,A1,A2,A3,Q,R,V,S,T,TEMP,THETA INTRINSIC ACOS, COS, SQRT C C SET VALUE OF PI PI=3.141592653589793D0 A1=B/A A2=C/A A3=D/A Q=(.3D01*A2-A1*A1)/.9D01 R=(.9D01*A1*A2-.27D02*A3-.2D01*A1*A1*A1)/.54D02 V=Q*Q*Q+R*R IF(V .GT. 0D0)THEN S=R+SQRT(V) T=R-SQRT(V) IF(T .LT. 0D0)THEN T=-(-T)**(1D0/3D0) ELSE T=T**(1D0/3D0) END IF IF(S .LT. 0)THEN S=-(-S)**(1D0/3D0) ELSE S=S**(1D0/3D0) END IF S1=S+T-A1/3D0 S3=S1 S2=S1 ELSE TEMP=R/SQRT(-Q**.3D01) THETA=ACOS(TEMP) THETA=THETA/.3D01 TEMP=.2D01*SQRT(-Q) S1=TEMP*COS(THETA)-A1/.3D01 S2=TEMP*COS(THETA+PI*.2D01/.3D01)-A1/.3D01 S3=TEMP*COS(THETA+PI*.4D01/.3D01)-A1/.3D01 END IF RETURN END C ---------------------- C | S O R T R T | C ---------------------- SUBROUTINE SORTRT(S1,S2,S3) C C PURPOSE: C SORT ROOTS INTO ASCENDING ORDER C C----------------------------------------------------------------------------- C C PARAMETERS: C S1 <--> ROOT C S2 <--> ROOT C S3 <--> ROOT C DOUBLE PRECISION S1,S2,S3 DOUBLE PRECISION T IF(S1 .GT. S2)THEN T=S1 S1=S2 S2=T END IF IF(S2 .GT. S3)THEN T=S2 S2=S3 S3=T END IF IF(S1 .GT. S2)THEN T=S1 S1=S2 S2=T END IF RETURN END C ---------------------- C | F S T O F D | C ---------------------- SUBROUTINE FSTOFD(NR,M,N,XPLS,FCN,FPLS,A,TYPX,RNOISE, Z FHAT,ICASE,NFCNT) 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 RNOISE --> RELATIVE NOISE IN FCN [F(X)] C FHAT(M) --> WORKSPACE C ICASE --> =1 OPTIMIZATION (GRADIENT) C =2 SYSTEMS C =3 OPTIMIZATION (HESSIAN) C NFCNT <--> NUMBER OF FUNCTION EVALUTIONS C C INTERNAL VARIABLES C ------------------ C STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION C INTEGER NR,M,N,ICASE,NFCNT DOUBLE PRECISION XPLS(N),FPLS(M),FHAT(M),TYPX(N),A(NR,N) DOUBLE PRECISION RNOISE INTEGER I,J,NM1,JP1 DOUBLE PRECISION XTMPJ,STEPSZ EXTERNAL FCN INTRINSIC ABS, MAX, SQRT 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 XTMPJ=XPLS(J) STEPSZ=SQRT(RNOISE)*MAX(ABS(XPLS(J)),1.D0) XPLS(J)=XTMPJ+STEPSZ CALL FCN(N,XPLS,FHAT) NFCNT = NFCNT + 1 XPLS(J)=XTMPJ DO 20 I=1,M A(I,J)=(FHAT(I)-FPLS(I))/STEPSZ A(I,J)=A(I,J)*TYPX(J) 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 C ---------------------- C | S N D O F D | C ---------------------- SUBROUTINE SNDOFD(NR,N,XPLS,FCN,FPLS,A,TYPX, Z RNOISE,STEPSZ,ANBR,NFCNT) 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 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 NFCNT <--> NUMBER OF FUNCTION EVALUATIONS C C INTEGER NR,N,NFCNT DOUBLE PRECISION XPLS(N),TYPX(N),STEPSZ(N),ANBR(N),A(NR,N) DOUBLE PRECISION FPLS,RNOISE,OV3,XTMPI,XTMPJ,FHAT INTEGER I,J,IP1 EXTERNAL FCN INTRINSIC ABS, MAX C C FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION C OF I-TH UNIT VECTOR. C OV3 = 1D0/3D0 DO 10 I=1,N XTMPI=XPLS(I) STEPSZ(I)=RNOISE**OV3 * MAX(ABS(XPLS(I)),1D0) XPLS(I)=XTMPI+STEPSZ(I) CALL FCN(N,XPLS,ANBR(I)) NFCNT = NFCNT + 1 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) NFCNT = NFCNT + 1 A(I,I)=((FPLS-ANBR(I))+(FHAT-ANBR(I)))/(STEPSZ(I)*STEPSZ(I)) A(I,I)=A(I,I)*(TYPX(I)*TYPX(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) NFCNT = NFCNT + 1 A(J,I)=((FPLS-ANBR(I))+(FHAT-ANBR(J)))/(STEPSZ(I)*STEPSZ(J)) A(J,I)=A(J,I)*(TYPX(I)*TYPX(J)) XPLS(J)=XTMPJ 20 CONTINUE 25 XPLS(I)=XTMPI 30 CONTINUE RETURN END C ---------------------- C | B A K S L V | C ---------------------- 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 INTEGER NR,N DOUBLE PRECISION A(NR,N),X(N),B(N) INTEGER I,IP1,J DOUBLE PRECISION SUM 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 C ---------------------- C | F O R S L V | C ---------------------- 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 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 THEN VECTORS B AND X MAY SHARE THE SAME STORAGE C INTEGER NR,N DOUBLE PRECISION A(NR,N),X(N),B(N) INTEGER I,J,IM1 DOUBLE PRECISION SUM 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 C ---------------------- C | C H O L D R | C ---------------------- SUBROUTINE CHOLDR(NR,N,H,G,EPS,PIVOT,E,DIAG,ADDMAX) C C PURPOSE: C DRIVER FOR CHOLESKY DECOMPOSITION C C---------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION C N --> DIMENSION OF PROBLEM C H(N,N) --> MATRIX C G(N) --> WORK SPACE C EPS --> MACHINE EPSILON C PIVOT(N) --> PIVOTING VECTOR C E(N) --> DIAGONAL MATRIX ADDED TO H FOR MAKING H P.D. C DIAG(N) --> DIAGONAL OF H C ADDMAX --> ADDMAX * I IS ADDED TO H INTEGER NR,N DOUBLE PRECISION H(NR,N),G(N),EPS,E(N),DIAG(N), Z ADDMAX,TEMP INTEGER PIVOT(N),I,J,K LOGICAL REDO DOUBLE PRECISION TAU1,TAU2 INTRINSIC SQRT REDO=.FALSE. C C SAVE DIAGONAL OF H DO 10 I=1,N DIAG(I)=H(I,I) 10 CONTINUE TAU1=EPS**(1D0/3D0) TAU2=TAU1 CALL MODCHL(NR,N,H,G,EPS,TAU1,TAU2,PIVOT,E) ADDMAX=E(N) DO 22 I=1,N IF(PIVOT(I) .NE. I)REDO=.TRUE. 22 CONTINUE IF((ADDMAX .GT. 0D0) .OR. REDO)THEN C******************************** C * C H IS NOT P.D. * C * C******************************** C C H=H+UI DO 30 I=2,N DO 32 J=1,I-1 H(I,J)=H(J,I) 32 CONTINUE 30 CONTINUE DO 34 I=1,N PIVOT(I)=I H(I,I)=DIAG(I)+ADDMAX 34 CONTINUE C******************************** C * C COMPUTE L * C * C******************************** DO 40 J=1,N C C COMPUTE L(J,J) TEMP=0D0 IF(J .GT. 1)THEN DO 41 I=1,J-1 TEMP=TEMP+H(J,I)*H(J,I) 41 CONTINUE END IF H(J,J)=H(J,J)-TEMP H(J,J)=SQRT(H(J,J)) C C COMPUTE L(I,J) DO 43 I=J+1,N TEMP=0D0 IF(J .GT. 1)THEN DO 45 K=1,J-1 TEMP=TEMP+H(I,K)*H(J,K) 45 CONTINUE END IF H(I,J)=H(J,I)-TEMP H(I,J)=H(I,J)/H(J,J) 43 CONTINUE 40 CONTINUE END IF RETURN END C ---------------------- C | M O D C H L | C ---------------------- C********************************************************************* C C SUBROUTINE NAME: MODCHL C C AUTHORS : ELIZABETH ESKOW AND ROBERT B. SCHNABEL C C DATE : DECEMBER, 1988 C C PURPOSE : PERFORM A MODIFIED CHOLESKY FACTORIZATION C OF THE FORM (PTRANSPOSE)AP + E = L(LTRANSPOSE), C WHERE L IS STORED IN THE LOWER TRIANGLE OF THE C ORIGINAL MATRIX A. C THE FACTORIZATION HAS 2 PHASES: C PHASE 1: PIVOT ON THE MAXIMUM DIAGONAL ELEMENT. C CHECK THAT THE NORMAL CHOLESKY UPDATE C WOULD RESULT IN A POSITIVE DIAGONAL C AT THE CURRENT ITERATION, AND C IF SO, DO THE NORMAL CHOLESKY UPDATE, C OTHERWISE SWITCH TO PHASE 2. C PHASE 2: PIVOT ON THE MINIMUM OF THE NEGATIVES C OF THE LOWER GERSCHGORIN BOUND C ESTIMATES. C COMPUTE THE AMOUNT TO ADD TO THE C PIVOT ELEMENT AND ADD THIS C TO THE PIVOT ELEMENT. C DO THE CHOLESKY UPDATE. C UPDATE THE ESTIMATES OF THE C GERSCHGORIN BOUNDS. C C INPUT : NDIM - LARGEST DIMENSION OF MATRIX THAT C WILL BE USED C C N - DIMENSION OF MATRIX A C C A - N*N SYMMETRIC MATRIX (ONLY LOWER TRIANGULAR C PORTION OF A, INCLUDING THE MAIN DIAGONAL, IS USED) C C G - N*1 WORK ARRAY C C MCHEPS - MACHINE PRECISION C C TAU1 - TOLERANCE USED FOR DETERMINING WHEN TO SWITCH TO C PHASE 2 C C TAU2 - TOLERANCE USED FOR DETERMINING THE MAXIMUM C CONDITION NUMBER OF THE FINAL 2X2 SUBMATRIX. C C C OUTPUT : L - STORED IN THE MATRIX A (IN LOWER TRIANGULAR C PORTION OF A, INCLUDING THE MAIN DIAGONAL) C C P - A RECORD OF HOW THE ROWS AND COLUMNS C OF THE MATRIX WERE PERMUTED WHILE C PERFORMING THE DECOMPOSITION C C E - N*1 ARRAY, THE ITH ELEMENT IS THE C AMOUNT ADDED TO THE DIAGONAL OF A C AT THE ITH ITERATION C C C************************************************************************ SUBROUTINE MODCHL(NDIM,N,A,G,MCHEPS,TAU1,TAU2,P,E) INTEGER N,NDIM DOUBLE PRECISION A(NDIM,N),G(N),MCHEPS,TAU1,TAU2 INTEGER P(N) DOUBLE PRECISION E(N) C C J - CURRENT ITERATION NUMBER C IMING - INDEX OF THE ROW WITH THE MIN. OF THE C NEG. LOWER GERSCH. BOUNDS C IMAXD - INDEX OF THE ROW WITH THE MAXIMUM DIAG. C ELEMENT C I,ITEMP,JPL,K - TEMPORARY INTEGER VARIABLES C DELTA - AMOUNT TO ADD TO AJJ AT THE JTH ITERATION C GAMMA - THE MAXIMUM DIAGONAL ELEMENT OF THE ORIGINAL C MATRIX A. C NORMJ - THE 1 NORM OF A(COLJ), ROWS J+1 --> N. C MING - THE MINIMUM OF THE NEG. LOWER GERSCH. BOUNDS C MAXD - THE MAXIMUM DIAGONAL ELEMENT C TAUGAM - TAU1 * GAMMA C PHASE1 - LOGICAL, TRUE IF IN PHASE1, OTHERWISE FALSE C DELTA1,TEMP,JDMIN,TDMIN,TEMPJJ - TEMPORARY DOUBLE PRECISION VARS. C INTEGER J,IMING,I,IMAXD,ITEMP,JP1,K DOUBLE PRECISION DELTA,GAMMA DOUBLE PRECISION NORMJ, MING,MAXD DOUBLE PRECISION DELTA1,TEMP,JDMIN,TDMIN,TAUGAM,TEMPJJ LOGICAL PHASE1 INTRINSIC ABS, MAX, SQRT, MIN CALL INIT(N, NDIM, A, PHASE1, DELTA, P, G, E, * MING,TAU1,GAMMA,TAUGAM) C CHECK FOR N=1 IF (N.EQ.1) THEN DELTA = (TAU2 * ABS(A(1,1))) - A(1,1) IF (DELTA .GT. 0) E(1)=DELTA IF (A(1,1) .EQ. 0) E(1) = TAU2 A(1,1)=SQRT(A(1,1)+E(1)) ENDIF C C DO 200 J = 1, N-1 C C PHASE 1 C IF ( PHASE1 ) THEN C C FIND INDEX OF MAXIMUM DIAGONAL ELEMENT A(I,I) WHERE I>=J C MAXD = A(J,J) IMAXD = J DO 20 I = J+1, N IF (MAXD .LT. A(I,I)) THEN MAXD = A(I,I) IMAXD = I END IF 20 CONTINUE C C PIVOT TO THE TOP THE ROW AND COLUMN WITH THE MAX DIAG C IF (IMAXD .NE. J) THEN C C SWAP ROW J WITH ROW OF MAX DIAG C DO 30 I = 1, J-1 TEMP = A(J,I) A(J,I) = A(IMAXD,I) A(IMAXD,I) = TEMP 30 CONTINUE C C SWAP COLJ AND ROW MAXDIAG BETWEEN J AND MAXDIAG C DO 35 I = J+1,IMAXD-1 TEMP = A(I,J) A(I,J) = A(IMAXD,I) A(IMAXD,I) = TEMP 35 CONTINUE C C SWAP COLUMN J WITH COLUMN OF MAX DIAG C DO 40 I = IMAXD+1, N TEMP = A(I,J) A(I,J) = A(I,IMAXD) A(I,IMAXD) = TEMP 40 CONTINUE C C SWAP DIAG ELEMENTS C TEMP = A(J,J) A(J,J) = A(IMAXD,IMAXD) A(IMAXD,IMAXD) = TEMP C C SWAP ELEMENTS OF THE PERMUTATION VECTOR C ITEMP = P(J) P(J) = P(IMAXD) P(IMAXD) = ITEMP END IF C CHECK TO SEE WHETHER THE NORMAL CHOLESKY UPDATE FOR THIS C ITERATION WOULD RESULT IN A POSITIVE DIAGONAL, C AND IF NOT THEN SWITCH TO PHASE 2. JP1 = J+1 TEMPJJ=A(J,J) IF (TEMPJJ.GT.0) THEN JDMIN=A(JP1,JP1) DO 60 I = JP1, N TEMP = A(I,J) * A(I,J) / TEMPJJ TDMIN = A(I,I) - TEMP JDMIN = MIN(JDMIN, TDMIN) 60 CONTINUE IF (JDMIN .LT. TAUGAM) PHASE1 = .FALSE. ELSE PHASE1 = .FALSE. END IF IF (PHASE1) THEN C C DO THE NORMAL CHOLESKY UPDATE IF STILL IN PHASE 1 C A(J,J) = SQRT(A(J,J)) TEMPJJ = A(J,J) DO 70 I = JP1, N A(I,J) = A(I,J) / TEMPJJ 70 CONTINUE DO 80 I=JP1,N TEMP=A(I,J) DO 75 K = JP1, I A(I,K) = A(I,K) - (TEMP * A(K,J)) 75 CONTINUE 80 CONTINUE IF (J .EQ. N-1) A(N,N)=SQRT(A(N,N)) ELSE C C CALCULATE THE NEGATIVES OF THE LOWER GERSCHGORIN BOUNDS C CALL GERSCH(NDIM,N,A,J,G) END IF END IF C C PHASE 2 C IF (.NOT. PHASE1) THEN IF (J .NE. N-1) THEN C C FIND THE MINIMUM NEGATIVE GERSHGORIN BOUND C IMING=J MING = G(J) DO 90 I = J+1,N IF (MING .GT. G(I)) THEN MING = G(I) IMING = I END IF 90 CONTINUE C C PIVOT TO THE TOP THE ROW AND COLUMN WITH THE C MINIMUM NEGATIVE GERSCHGORIN BOUND C IF (IMING .NE. J) THEN C C SWAP ROW J WITH ROW OF MIN GERSCH BOUND C DO 100 I = 1, J-1 TEMP = A(J,I) A(J,I) = A(IMING,I) A(IMING,I) = TEMP 100 CONTINUE C C SWAP COLJ WITH ROW IMING FROM J TO IMING C DO 105 I = J+1,IMING-1 TEMP = A(I,J) A(I,J) = A(IMING,I) A(IMING,I) = TEMP 105 CONTINUE C C SWAP COLUMN J WITH COLUMN OF MIN GERSCH BOUND C DO 110 I = IMING+1, N TEMP = A(I,J) A(I,J) = A(I,IMING) A(I,IMING) = TEMP 110 CONTINUE C C SWAP DIAGONAL ELEMENTS C TEMP = A(J,J) A(J,J) = A(IMING,IMING) A(IMING,IMING) = TEMP C C SWAP ELEMENTS OF THE PERMUTATION VECTOR C ITEMP = P(J) P(J) = P(IMING) P(IMING) = ITEMP C C SWAP ELEMENTS OF THE NEGATIVE GERSCHGORIN BOUNDS VECTOR C TEMP = G(J) G(J) = G(IMING) G(IMING) = TEMP END IF C C CALCULATE DELTA AND ADD TO THE DIAGONAL. C DELTA=MAX{0,-A(J,J) + MAX{NORMJ,TAUGAM},DELTA_PREVIOUS} C WHERE NORMJ=SUM OF |A(I,J)|,FOR I=1,N, C DELTA_PREVIOUS IS THE DELTA COMPUTED AT THE PREVIOUS ITERATION, C AND TAUGAM IS TAU1*GAMMA. C NORMJ = 0.0 DO 140 I = J+1, N NORMJ = NORMJ + ABS(A(I,J)) 140 CONTINUE TEMP = MAX(NORMJ,TAUGAM) DELTA1 = TEMP - A(J,J) TEMP = 0.0 DELTA1 = MAX(TEMP, DELTA1) DELTA = MAX(DELTA1,DELTA) E(J) = DELTA A(J,J) = A(J,J) + E(J) C C UPDATE THE GERSCHGORIN BOUND ESTIMATES C (NOTE: G(I) IS THE NEGATIVE OF THE C GERSCHGORIN LOWER BOUND.) C IF (A(J,J) .NE. NORMJ) THEN TEMP = (NORMJ/A(J,J)) - 1.0 DO 150 I = J+1, N G(I) = G(I) + ABS(A(I,J)) * TEMP 150 CONTINUE END IF C C DO THE CHOLESKY UPDATE C A(J,J) = SQRT(A(J,J)) TEMPJJ = A(J,J) DO 160 I = J+1, N A(I,J) = A(I,J) / TEMPJJ 160 CONTINUE DO 180 I = J+1, N TEMP = A(I,J) DO 170 K = J+1, I A(I,K) = A(I,K) - (TEMP * A(K,J)) 170 CONTINUE 180 CONTINUE ELSE CALL FIN2X2(NDIM, N, A, E, J, TAU2, DELTA,GAMMA) END IF END IF 200 CONTINUE RETURN END C************************************************************************ C SUBROUTINE NAME : INIT C C PURPOSE : SET UP FOR START OF CHOLESKY FACTORIZATION C C INPUT : N, NDIM, A, TAU1 C C OUTPUT : PHASE1 - BOOLEAN VALUE SET TO TRUE IF IN PHASE ONE, C OTHERWISE FALSE. C DELTA - AMOUNT TO ADD TO AJJ AT ITERATION J C P,G,E - DESCRIBED ABOVE IN MODCHL C MING - THE MINIMUM NEGATIVE GERSCHGORIN BOUND C GAMMA - THE MAXIMUM DIAGONAL ELEMENT OF A C TAUGAM - TAU1 * GAMMA C C************************************************************************ SUBROUTINE INIT(N,NDIM,A,PHASE1,DELTA,P,G,E,MING, * TAU1,GAMMA,TAUGAM) INTEGER N,NDIM DOUBLE PRECISION A(NDIM,N) LOGICAL PHASE1 DOUBLE PRECISION DELTA,G(N),E(N) INTEGER P(N),I DOUBLE PRECISION MING,TAU1,GAMMA,TAUGAM INTRINSIC ABS, MAX PHASE1 = .TRUE. DELTA = 0.0 MING = 0.0 DO 10 I=1,N P(I)=I G(I)= 0.0 E(I) = 0.0 10 CONTINUE C C FIND THE MAXIMUM MAGNITUDE OF THE DIAGONAL ELEMENTS. C IF ANY DIAGONAL ELEMENT IS NEGATIVE, THEN PHASE1 IS FALSE. C GAMMA = 0.0 DO 20 I=1,N GAMMA=MAX(GAMMA,ABS(A(I,I))) IF (A(I,I) .LT. 0.0) PHASE1 = .FALSE. 20 CONTINUE TAUGAM = TAU1 * GAMMA C C IF NOT IN PHASE1, THEN CALCULATE THE INITIAL GERSCHGORIN BOUNDS C NEEDED FOR THE START OF PHASE2. C IF ( .NOT.PHASE1) CALL GERSCH(NDIM,N,A,1,G) RETURN END C************************************************************************ C C SUBROUTINE NAME : GERSCH C C PURPOSE : CALCULATE THE NEGATIVE OF THE GERSCHGORIN BOUNDS C CALLED ONCE AT THE START OF PHASE II. C C INPUT : NDIM, N, A, J C C OUTPUT : G - AN N VECTOR CONTAINING THE NEGATIVES OF THE C GERSCHGORIN BOUNDS. C C************************************************************************ SUBROUTINE GERSCH(NDIM, N, A, J, G) INTEGER NDIM, N, J DOUBLE PRECISION A(NDIM,N), G(N) INTEGER I, K DOUBLE PRECISION OFFROW INTRINSIC ABS DO 30 I = J, N OFFROW = 0.0 DO 10 K = J, I-1 OFFROW = OFFROW + ABS(A(I,K)) 10 CONTINUE DO 20 K = I+1, N OFFROW = OFFROW + ABS(A(K,I)) 20 CONTINUE G(I) = OFFROW - A(I,I) 30 CONTINUE RETURN END C************************************************************************ C C SUBROUTINE NAME : FIN2X2 C C PURPOSE : HANDLES FINAL 2X2 SUBMATRIX IN PHASE II. C FINDS EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX, C CALCULATES THE AMOUNT TO ADD TO THE DIAGONAL, C ADDS TO THE FINAL 2 DIAGONAL ELEMENTS, C AND DOES THE FINAL UPDATE. C C INPUT : NDIM, N, A, E, J, TAU2, C DELTA - AMOUNT ADDED TO THE DIAGONAL IN THE C PREVIOUS ITERATION C C OUTPUT : A - MATRIX WITH COMPLETE L FACTOR IN THE LOWER TRIANGLE, C E - N*1 VECTOR CONTAINING THE AMOUNT ADDED TO THE DIAGONAL C AT EACH ITERATION, C DELTA - AMOUNT ADDED TO DIAGONAL ELEMENTS N-1 AND N. C C************************************************************************ SUBROUTINE FIN2X2(NDIM, N, A, E, J, TAU2, DELTA,GAMMA) INTEGER NDIM, N, J DOUBLE PRECISION A(NDIM,N), E(N), TAU2, DELTA,GAMMA DOUBLE PRECISION T1, T2, T3,LMBD1,LMBD2,LMBDHI,LMBDLO DOUBLE PRECISION DELTA1, TEMP,T1A,T2A INTRINSIC SQRT, MAX, MIN C C FIND EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX C T1 = A(N-1,N-1) + A(N,N) T2 = A(N-1,N-1) - A(N,N) T1A=ABS(T2) T2A= 2.D0 * ABS(A(N,N-1)) IF (T1A .GE. T2A) THEN IF (T1A .GT. 0) T2A = T2A / T1A T3 = T1A * SQRT(1.D0 + (T2A**2)) ELSE T1A = T1A / T2A T3 = T2A * SQRT(1.D0 + (T1A**2)) ENDIF LMBD1 = (T1 - T3)/2. LMBD2 = (T1 + T3)/2. LMBDHI = MAX(LMBD1,LMBD2) LMBDLO = MIN(LMBD1,LMBD2) C C FIND DELTA SUCH THAT: C 1. THE L2 CONDITION NUMBER OF THE FINAL C 2X2 SUBMATRIX + DELTA*I <= TAU2 C 2. DELTA >= PREVIOUS DELTA, C 3. LMBDLO + DELTA >= TAU2 * GAMMA, C WHERE LMBDLO IS THE SMALLEST EIGENVALUE OF THE FINAL C 2X2 SUBMATRIX C DELTA1=(LMBDHI-LMBDLO)/(1.0-TAU2) DELTA1= MAX(DELTA1,GAMMA) DELTA1= TAU2 * DELTA1 - LMBDLO TEMP = 0.0 DELTA = MAX(DELTA, TEMP) DELTA = MAX(DELTA, DELTA1) IF (DELTA .GT. 0.0) THEN A(N-1,N-1) = A(N-1,N-1) + DELTA A(N,N) = A(N,N) + DELTA E(N-1) = DELTA E(N) = DELTA END IF C C FINAL UPDATE C A(N-1,N-1) = SQRT(A(N-1,N-1)) A(N,N-1) = A(N,N-1)/A(N-1,N-1) A(N,N) = A(N,N) - (A(N,N-1)**2) A(N,N) = SQRT(A(N,N)) RETURN END C ---------------------- C | S L V M D L | C ---------------------- SUBROUTINE SLVMDL(NR,N,H,U,T,E,DIAG,S,G, Z PIVOT,W1,W2,W3,ALPHA,BETA,NOMIN,EPS) C C PURPOSE: C COMPUTE TENSOR AND NEWTON STEPS C C---------------------------------------------------------------------------- C C PARAMETERS: C C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C H(N,N) --> HESSIAN C U(N) --> VECTOR TO FORM Q IN QR C T(N) --> WORKSPACE C E(N) --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION C DIAG(N) --> DIAGONAL OF HESSIAN C S(N) --> STEP TO PREVIOUS POINT (FOR TENSOR MODEL) C G(N) --> CURRENT GRADIENT C PIVOT(N) --> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION C W1(N) <--> ON INPUT: A=2*(GP-G-HS-S*BETA/(6*STS)) C ON OUTPUT: TENSOR STEP C W2(N) --> SH C W3(N) <-- NEWTON STEP C ALPHA --> SCALAR FOR 3RD ORDER TERM OF TENSOR MODEL C BETA --> SCALAR FOR 4TH ORDER TERM OF TENSOR MODEL C NOMIN <-- =.TRUE. IF TENSOR MODEL HAS NO MINIMIZER C EPS --> MACHINE EPSILON C INTEGER NR,N,PIVOT(N) DOUBLE PRECISION H(NR,N),U(N),T(N),E(N),DIAG(N),S(N), Z G(N),W1(N),W2(N),W3(N),ALPHA,BETA,EPS INTEGER I,J DOUBLE PRECISION W11,W22,W33,SG,ADDMAX,SHS,CA,W12 DOUBLE PRECISION W13,W23,SGSTAR,CB,CC,CD,R,R2,R1,TEMP DOUBLE PRECISION UU,SS LOGICAL NOMIN C C S O L V E M O D E L C NOMIN=.FALSE. C C COMPUTE QTHQ(N,N), ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND N-1 C COLUMNS OF QTHQ IF (N.GT.1) THEN CALL ZHZ(NR,N,S,H,U,T) C C IN CHOLESKY DECOMPOSITION WILL STORE H(1,1) ... H(N-1,N-1) C IN DIAG(1) ... DIAG(N-1), STORE H(N,N) IN DIAG(N) FIRST DIAG(N)=H(N,N) C C COLESKY DECOMPOSITION FOR FIRST N-1 ROWS AND N-1 COLUMNS OF ZTHZ C ZTHZ(N-1,N-1)=LLT CALL CHOLDR(NR,N-1,H,T,EPS,PIVOT,E,DIAG,ADDMAX) END IF C C ON INPUT: SH IS STORED IN W2 SHS=0D0 DO 100 I=1,N SHS=SHS+W2(I)*S(I) W3(I)=G(I) 100 CONTINUE C C COMPUTE W1,W2,W3 C W1=L**-1*ZT*A C W2=L**-1*ZT*SH C W3=L**-1*ZT*G C IF (N.GT.1) THEN CALL SOLVEW(NR,N,H,U,W1,T) CALL SOLVEW(NR,N,H,U,W2,T) CALL SOLVEW(NR,N,H,U,W3,T) END IF C C COMPUTE COEFFICIENTS CA,CB,CC AND CD FOR REDUCED ONE VARIABLE C 3RD ORDER EQUATION W11=0D0 DO 110 I=1,N-1 W11=W11+W1(I)*W1(I) 110 CONTINUE CA=BETA/6D0-W11/2D0 W12=0D0 DO 120 I=1,N-1 W12=W12+W1(I)*W2(I) 120 CONTINUE CB=ALPHA/2D0-3D0*W12/2D0 W13=0D0 DO 130 I=1,N-1 W13=W13+W1(I)*W3(I) 130 CONTINUE W22=0D0 DO 133 I=1,N-1 W22=W22+W2(I)*W2(I) 133 CONTINUE CC=SHS-W22-W13 SG=0D0 DO 140 I=1,N SG=SG+S(I)*G(I) 140 CONTINUE W23=0D0 DO 145 I=1,N-1 W23=W23+W2(I)*W3(I) 145 CONTINUE CD=SG-W23 W33=0D0 DO 147 I=1,N-1 W33=W33+W3(I)*W3(I) 147 CONTINUE C C COMPUTE DESIRABLE ROOT, SGSTAR, OF 3RD ORDER EQUATION IF(CA .NE. 0D0)THEN CALL SIGMA(SGSTAR,CA,CB,CC,CD) IF(CA .EQ. 0D0)THEN NOMIN=.TRUE. GO TO 200 END IF ELSE C 2ND ORDER ( CA=0 ) IF(CB .NE. 0D0)THEN R=CC*CC-4D0*CB*CD IF(R .LT. 0D0)THEN NOMIN=.TRUE. GO TO 200 ELSE R1=(-CC+SQRT(R))/(2D0*CB) R2=(-CC-SQRT(R))/(2D0*CB) IF(R2 .LT. R1)THEN TEMP=R1 R1=R2 R2=TEMP END IF IF(CB .GT. 0D0)THEN SGSTAR=R2 ELSE SGSTAR=R1 END IF IF( ( (R1 .GT. 0D0) .AND. (SGSTAR .EQ. R2) ) .OR. Z ( (R2 .LT. 0D0) .AND. (SGSTAR .EQ. R1) ) )THEN NOMIN=.TRUE. GO TO 200 END IF END IF ELSE C 1ST ORDER (CA=0,CB=0) IF(CC .GT. 0D0)THEN SGSTAR=-CD/CC ELSE NOMIN=.TRUE. GO TO 200 END IF END IF END IF C C FIND TENSOR STEP, W1 (FUNCTION OF SGSTAR) CALL DSTAR(NR,N,U,S,W1,W2,W3,SGSTAR,H,W1) C C COMPUTE DN 200 UU=0D0 SS=0D0 DO 202 I=1,N UU=UU+U(I)*U(I) SS=SS+S(I)*S(I) 202 CONTINUE UU=UU/2D0 SS=DSQRT(SS) IF (N.EQ.1) THEN CALL CHOLDR(NR,N,H,T,EPS,PIVOT,E,DIAG,ADDMAX) ELSE C C COMPUTE LAST ROW OF L(N,N) DO 220 I=1,N-1 TEMP=0D0 IF(I .GT. 1)THEN DO 210 J=1,I-1 TEMP=TEMP+H(N,J)*H(I,J) 210 CONTINUE END IF H(N,I)=(H(I,N)-TEMP)/H(I,I) 220 CONTINUE TEMP=0D0 DO 224 I=1,N-1 TEMP=TEMP+H(N,I)*H(N,I) 224 CONTINUE H(N,N)=DIAG(N)-TEMP+ADDMAX IF(H(N,N) .GT. 0D0)THEN H(N,N)=DSQRT(H(N,N)) ELSE C AFTER ADDING THE LAST COLUMN AND ROW C QTHQ IS NOT P.D., NEED TO REDO CHOLESKY DECOMPOSITION DO 232 I=2,N DO 230 J=1,I-1 H(I,J)=H(J,I) 230 CONTINUE H(I,I)=DIAG(I) 232 CONTINUE H(1,1)=DIAG(1) CALL CHOLDR(NR,N,H,T,EPS,PIVOT,E,DIAG,ADDMAX) END IF END IF C C SOLVE QTHQ*QT*W3=-QT*G, WHERE W3 IS NEWTON STEP C W2=-QT*G DO 302 I=1,N W2(I)=0D0 DO 300 J=1,N W2(I)=W2(I)+U(J)*U(I)*G(J)/UU 300 CONTINUE W2(I)=W2(I)-G(I) 302 CONTINUE CALL FORSLV(NR,N,H,W3,W2) CALL BAKSLV(NR,N,H,W2,W3) C W2=QT*W3 => W3=Q*W2 --- NEWTON STEP DO 312 I=1,N W3(I)=0D0 DO 310 J=1,N W3(I)=W3(I)+U(I)*U(J)*W2(J)/UU 310 CONTINUE W3(I)=W2(I)-W3(I) 312 CONTINUE RETURN END C ---------------------- C | O P T S T P | C ---------------------- SUBROUTINE OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX, + ITRMCD,GRADTL,STEPTL,FSCALE,ITNLIM,IRETCD,MXTAKE,IPR,MSG, + RGX,RSX) 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 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 <--> CONTROL OUTPUT ON INPUT AND CONTAIN STOPPING C CONDITION ON OUTPUT C C INTEGER N,ITNCNT,ICSCMX,ITRMCD,ITNLIM,IRETCD,IPR,MSG DOUBLE PRECISION XPLS(N),GPLS(N),X(N),FPLS,GRADTL,STEPTL,FSCALE INTEGER I,JTRMCD,IMSG DOUBLE PRECISION RGX,D,RELGRD,RELSTP,RSX LOGICAL MXTAKE INTRINSIC ABS, MAX C ITRMCD=0 IMSG=MSG RGX=0.D0 RSX=0.D0 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) C D=1D0 RGX=0.0 DO 100 I=1,N RELGRD=ABS(GPLS(I))*MAX(ABS(XPLS(I)),1.D0)/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.0 DO 120 I=1,N RELSTP=ABS(XPLS(I)-X(I))/MAX(ABS(XPLS(I)),1.D0) 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 C IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900) IF(IMSG .GE. 1) 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 C IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD IF (IMSG .GE. 1) 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 MSG=-ITRMCD RETURN C 900 FORMAT(48H OPTSTP STEP OF MAXIMUM LENGTH (STEPMX) TAKEN) 901 FORMAT(' OPTSTP RELATIVE GRADIENT CLOSE TO ZERO.'/ + ' OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.') 902 FORMAT(48H OPTSTP SUCCESSIVE ITERATES WITHIN TOLERANCE./ + 48H OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.) 903 FORMAT(52H OPTSTP 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(36H OPTSTP ITERATION LIMIT EXCEEDED./ + 28H OPTSTP ALGORITHM FAILED.) 905 FORMAT(39H OPTSTP 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 C SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*),DY(*) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 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 DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END C DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 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 DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + 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 60 DDOT = DTEMP RETURN END C SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DA,DX(*) INTEGER I,INCX,M,MP1,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END c*** wood.inp -3,-1,-3,-1 c*** samp.out __________________________________________ f(x)= Wood function __________________________________________ INITIAL APPROXIMATION TO THE SOLUTION X*: -3.0000000000000 -1.0000000000000 -3.0000000000000 -1.0000000000000 CALLING TENSRD - ALL INPUT PARAMETERS ARE SET TO THEIR DEFAULT VALUES. OUTPUT WILL BE WRITTEN TO THE STANDARD OUTPUT. INITIAL FUNCTION VALUE F= 0.1919200000000E+05 INITIAL GRADIENT G=-0.1200799943060E+05-0.2080000057570E+04-0.1080799948883E+05-0.1880000003362E+04 OPTSTP RELATIVE GRADIENT CLOSE TO ZERO. OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION. FINAL X= 0.9999996220733E+00 0.9999992737728E+00 0.1000000116495E+01 0.1000000265519E+01 GRADIENT G= 0.7442482819355E-07-0.5611498535672E-08-0.6145918548890E-07 0.4772662607369E-08 FUNCTION VALUE F(X)= 0.2560272516251E-11 AT ITERATION 33 RESULTS OF TENSRD: G= 0.7442482819355E-07-0.5611498535672E-08-0.6145918548890E-07 0.4772662607369E-08 XPLS= 0.99999962207329 0.99999927377285 1.0000001164954 1.0000002655194 FPLS= 2.5602725162505D-12 ITNNO= 33 MSG= -1 __________________________________________ f(x)= Wood function __________________________________________ INITIAL APPROXIMATION TO THE SOLUTION X*: -3.0000000000000 -1.0000000000000 -3.0000000000000 -1.0000000000000 CALLING TENSOR AFTER RESETTING INPUT PARAMETERS IPR, MSG AND METHOD. THE STANDARD METHOD IS USED IN THIS EXAMPLE. ADDITIONAL OUTPUT WILL BE WRITTEN TO FILE 'FORT8'. RESULTS OF TENSOR: G=-0.1640905875818E-06 0.1150812112878E-09 0.3063518809213E-06-0.5242027244489E-07 XPLS= 0.99999955541473 0.99999914071878 1.0000001835075 1.0000003988941 FPLS= 2.6899952376097D-12 ITNNO= 63 MSG= -1 c*** fort8 INITIAL FUNCTION VALUE F= 0.1919200000000E+05 INITIAL GRADIENT G=-0.1200799943060E+05-0.2080000057570E+04-0.1080799948883E+05-0.1880000003362E+04 ----------- ITERATION 1 ---------------- X=-0.2696343145685E+01 0.6172014725401E+01-0.2662886565912E+01 0.5860562423349E+01 F(X)= 0.1290481213184E+04 G(X)=-0.1191897729121E+04-0.1893647244555E+02-0.1186837501761E+04-0.2088316852776E+02 ----------- ITERATION 2 ---------------- X=-0.1896590009595E+01 0.2637182602354E+01-0.1879844351575E+01 0.2564967338913E+01 F(X)= 0.2958370305757E+03 G(X)=-0.7339858283615E+03-0.1279167608874E+03-0.6614212678077E+03-0.1103639764491E+03 ----------- ITERATION 3 ---------------- X=-0.1524535636069E+01 0.2004180871016E+01-0.1474894868019E+01 0.1809980863757E+01 F(X)= 0.6766798280681E+02 G(X)=-0.2002066803801E+03-0.2768352524620E+02-0.1989282815207E+03-0.2951572109524E+02 ----------- ITERATION 4 ---------------- X=-0.1215211143448E+01 0.1323413666211E+01-0.1183662277287E+01 0.1252349336508E+01 F(X)= 0.1733215643731E+02 G(X)=-0.7895903385440E+02-0.1913541382490E+02-0.6773411774559E+02-0.1526621791561E+02 ----------- ITERATION 5 ---------------- X=-0.1065831508517E+01 0.1103344764211E+01-0.1018583627212E+01 0.9990070147813E+00 F(X)= 0.8688244006022E+01 G(X)=-0.1805227673685E+02-0.4462501104270E+01-0.1815677411407E+02-0.4904835175276E+01 ----------- ITERATION 6 ---------------- X=-0.1000244200012E+01 0.1003564626144E+01-0.9576925727155E+00 0.9218464079588E+00 F(X)= 0.7892750205061E+01 G(X)=-0.2769708705120E+01-0.8601988759784E+00-0.2304838579320E+01-0.6672778563380E+00 ----------- ITERATION 7 ---------------- X=-0.9966593994294E+00 0.1003309174287E+01-0.9405563611904E+00 0.8956941935882E+00 F(X)= 0.7876512616767E+01 G(X)=-0.1495457772274E-01-0.2562993473623E-02-0.1402720650969E+00-0.5282599452162E-01 ----------- ITERATION 8 ---------------- X=-0.1007734502375E+01 0.1025337533482E+01-0.9287317852339E+00 0.8735612409246E+00 F(X)= 0.7876070070247E+01 G(X)=-0.6162745853460E-01-0.2992445050859E-01-0.1734865719497E+00-0.6904444420638E-01 ----------- ITERATION 9 ---------------- X=-0.1022382409417E+01 0.1054809391676E+01-0.9127183675228E+00 0.8439574548745E+00 F(X)= 0.7875255635463E+01 G(X)=-0.1418675842992E+00-0.7376886418672E-01-0.2430541584295E+00-0.1043556769800E+00 ----------- ITERATION 10 ---------------- X=-0.1044655483543E+01 0.1100293519237E+01-0.8875638901643E+00 0.7982532050856E+00 F(X)= 0.7873534623919E+01 G(X)=-0.3333673290206E+00-0.1709655976011E+00-0.4253847726666E+00-0.2024321299980E+00 ----------- ITERATION 11 ---------------- X=-0.1089951783105E+01 0.1194844179397E+01-0.8337651473323E+00 0.7031754529189E+00 F(X)= 0.7869234157333E+01 G(X)=-0.1193728866734E+01-0.5714114825626E+00-0.1262937406237E+01-0.6959341497953E+00 ----------- ITERATION 12 ---------------- X=-0.1112911528943E+01 0.1245364461880E+01-0.8027068658104E+00 0.6522622913650E+00 F(X)= 0.7864391715010E+01 G(X)=-0.1202073759010E+01-0.5703620258469E+00-0.1315578787092E+01-0.7397659787904E+00 ----------- ITERATION 13 ---------------- X=-0.1128613321687E+01 0.1280557741197E+01-0.7802757232921E+00 0.6167735935437E+00 F(X)= 0.7860134598756E+01 G(X)=-0.1192024972202E+01-0.5626697728306E+00-0.1329252388072E+01-0.7563169465346E+00 ----------- ITERATION 14 ---------------- X=-0.1148075791369E+01 0.1324721030608E+01-0.7513606848874E+00 0.5722133368515E+00 F(X)= 0.7854054071844E+01 G(X)=-0.1245461864328E+01-0.5822049050658E+00-0.1427933940260E+01-0.8311285947498E+00 ----------- ITERATION 15 ---------------- X=-0.1169183115974E+01 0.1373426437727E+01-0.7182232003292E+00 0.5230397706105E+00 F(X)= 0.7845957809356E+01 G(X)=-0.1327802843694E+01-0.6131379256001E+00-0.1576045780752E+01-0.9456130824600E+00 ----------- ITERATION 16 ---------------- X=-0.1190261672835E+01 0.1422972596696E+01-0.6829009214181E+00 0.4729931681798E+00 F(X)= 0.7835752418713E+01 G(X)=-0.1404968345385E+01-0.6407344828888E+00-0.1733513177925E+01-0.1075567482073E+01 ----------- ITERATION 17 ---------------- X=-0.1210980056626E+01 0.1472565305038E+01-0.6456511905259E+00 0.4228835361952E+00 F(X)= 0.7823064819746E+01 G(X)=-0.1470727215770E+01-0.6625602093614E+00-0.1892489447062E+01-0.1217702676275E+01 ----------- ITERATION 18 ---------------- X=-0.1231124910547E+01 0.1521639569001E+01-0.6065921587353E+00 0.3732906585503E+00 F(X)= 0.7807452284149E+01 G(X)=-0.1521795849396E+01-0.6775156682322E+00-0.2047807315378E+01-0.1370472044620E+01 ----------- ITERATION 19 ---------------- X=-0.1250493612657E+01 0.1569622711170E+01-0.5658618849171E+00 0.3248024238430E+00 F(X)= 0.7788437265437E+01 G(X)=-0.1555581963383E+01-0.6848405557371E+00-0.2194092399000E+01-0.1531963025467E+01 ----------- ITERATION 20 ---------------- X=-0.1268901199831E+01 0.1615956618861E+01-0.5236085333407E+00 0.2779913386268E+00 F(X)= 0.7765516141790E+01 G(X)=-0.1570393279718E+01-0.6841693397987E+00-0.2326121645020E+01-0.1700051059009E+01 ----------- ITERATION 21 ---------------- X=-0.1286185244078E+01 0.1660116884323E+01-0.4799795542923E+00 0.2333943739186E+00 F(X)= 0.7738167966887E+01 G(X)=-0.1565549953957E+01-0.6755440973855E+00-0.2439159248041E+01-0.1872595884989E+01 ----------- ITERATION 22 ---------------- X=-0.1302206254371E+01 0.1701621616361E+01-0.4351163147713E+00 0.1915038154912E+00 F(X)= 0.7705864576194E+01 G(X)=-0.1541341558386E+01-0.6593644431287E+00-0.2529125900440E+01-0.2047542345755E+01 ----------- ITERATION 23 ---------------- X=-0.1316846300606E+01 0.1740035393233E+01-0.3891504167230E+00 0.1527628420389E+00 F(X)= 0.7668080911694E+01 G(X)=-0.1498930193636E+01-0.6363320631518E+00-0.2592702993902E+01-0.2223023520501E+01 ----------- ITERATION 24 ---------------- X=-0.1330005758629E+01 0.1774967558347E+01-0.3422043628988E+00 0.1175673587816E+00 F(X)= 0.7624306668524E+01 G(X)=-0.1440175960822E+01-0.6073673629126E+00-0.2627303072248E+01-0.2397342693630E+01 ----------- ITERATION 25 ---------------- X=-0.1341602421722E+01 0.1806075577421E+01-0.2943850103556E+00 0.8626288522413E-01 F(X)= 0.7574051877334E+01 G(X)=-0.1367527670222E+01-0.5735581007652E+00-0.2631123196102E+01-0.2569126876474E+01 ----------- ITERATION 26 ---------------- X=-0.1351565851636E+01 0.1833054914512E+01-0.2457921188183E+00 0.5915500706487E-01 F(X)= 0.7516861795174E+01 G(X)=-0.1283820904231E+01-0.5360826097722E+00-0.2602964963524E+01-0.2737154916893E+01 ----------- ITERATION 27 ---------------- X=-0.1359835150229E+01 0.1855636722226E+01-0.1965207269757E+00 0.3651214118439E-01 F(X)= 0.7452329548211E+01 G(X)=-0.1192179116998E+01-0.4961740119998E+00-0.2542194626943E+01-0.2900330311379E+01 ----------- ITERATION 28 ---------------- X=-0.1366358912582E+01 0.1873590815221E+01-0.1466522456264E+00 0.1856290145042E-01 F(X)= 0.7380094725628E+01 G(X)=-0.1095909692619E+01-0.4550860592035E+00-0.2448731072063E+01-0.3057844388593E+01 ----------- ITERATION 29 ---------------- X=-0.1371090536467E+01 0.1886715474159E+01-0.9627053650711E-01 0.5506531139604E-02 F(X)= 0.7299868201853E+01 G(X)=-0.9984049867683E+00-0.4140686091823E+00-0.2322904148709E+01-0.3208865839344E+01 ----------- ITERATION 30 ---------------- X=-0.1373989294353E+01 0.1894842380227E+01-0.4545096902266E-01-0.2491501099488E-02 F(X)= 0.7211427138664E+01 G(X)=-0.9030844920595E+00-0.3743491391052E+00-0.2165469837631E+01-0.3352758444434E+01 ----------- ITERATION 31 ---------------- X=-0.1375018090829E+01 0.1897832244674E+01 0.5732204010298E-02-0.5295234213261E-02 F(X)= 0.7114629687205E+01 G(X)=-0.8133295193951E+00-0.3371287063836E+00-0.1977540483218E+01-0.3488938675253E+01 ----------- ITERATION 32 ---------------- X=-0.1374143159008E+01 0.1895575460495E+01 0.5720800835193E-01-0.2798248174670E-02 F(X)= 0.7009416973882E+01 G(X)=-0.7324360827179E+00-0.3035666380786E+00-0.1760552212966E+01-0.3616908161970E+01 ----------- ITERATION 33 ---------------- X=-0.1371333298053E+01 0.1887991320539E+01 0.1089048003545E+00 0.5077081778205E-02 F(X)= 0.6895823429310E+01 G(X)=-0.6635722430615E+00-0.2747813455817E+00-0.1516250894142E+01-0.3736182690431E+01 ----------- ITERATION 34 ---------------- X=-0.1366385379729E+01 0.1874555349309E+01 0.1626449942376E+00 0.1886485249949E-01 F(X)= 0.6769771893514E+01 G(X)=-0.6082531445074E+00-0.2511826582973E+00-0.1230383778677E+01-0.3868668288628E+01 ----------- ITERATION 35 ---------------- X=-0.1358748554499E+01 0.1853825502699E+01 0.2205667539493E+00 0.4002795595178E-01 F(X)= 0.6625376493808E+01 G(X)=-0.5717233289456E+00-0.2345912029583E+00-0.8742651663715E+00-0.4037599933322E+01 ----------- ITERATION 36 ---------------- X=-0.1347414898287E+01 0.1823183565384E+01 0.2850585034920E+00 0.7121762029085E-01 F(X)= 0.6454883514836E+01 G(X)=-0.5681211780236E+00-0.2302453072641E+00-0.3994915051491E+00-0.4269697707356E+01 ----------- ITERATION 37 ---------------- X=-0.1330406704622E+01 0.1777548144048E+01 0.3602221982708E+00 0.1175798331496E+00 F(X)= 0.6244718741784E+01 G(X)=-0.6343640035530E+00-0.2522116979829E+00 0.2999741499022E+00-0.4621866716269E+01 ----------- ITERATION 38 ---------------- X=-0.1303299982041E+01 0.1705726922282E+01 0.4542679040664E+00 0.1904558103964E+00 F(X)= 0.5968254564877E+01 G(X)=-0.8863912822887E+00-0.3460694111202E+00 0.1509343146729E+01-0.5242029721105E+01 ----------- ITERATION 39 ---------------- X=-0.1253716162569E+01 0.1577146569337E+01 0.5874350400206E+00 0.3208902621295E+00 F(X)= 0.5566745916050E+01 G(X)=-0.1828289682943E+01-0.7195360113415E+00 0.4290422410101E+01-0.6644651057617E+01 ----------- ITERATION 40 ---------------- X=-0.1101712237320E+01 0.1199703801791E+01 0.8842784661716E+00 0.7049761186993E+00 F(X)= 0.5098931875022E+01 G(X)=-0.1040210412639E+02-0.4620662319731E+01 0.2427194321962E+02-0.1586035568094E+02 ----------- ITERATION 41 ---------------- X=-0.9843404511569E+00 0.9626637058192E+00 0.9814818416140E+00 0.9622072795127E+00 F(X)= 0.3998424212016E+01 G(X)=-0.6434409078568E+01-0.2754969110606E+01 0.3514033101767E+00-0.1700547069561E+01 ----------- ITERATION 42 ---------------- X=-0.9224273654140E+00 0.8426215927111E+00 0.1045141122290E+01 0.1087758158060E+01 F(X)= 0.3760923914087E+01 G(X)=-0.6889094635413E+01-0.3091559175172E+01 0.1806679049754E+01-0.2164499541908E+01 ----------- ITERATION 43 ---------------- X=-0.8220130237414E+00 0.6616265176119E+00 0.1134579377661E+01 0.1277043721858E+01 F(X)= 0.3442564400084E+01 G(X)=-0.8273230956239E+01-0.4165453871487E+01 0.4446232970179E+01-0.2944303345622E+01 ----------- ITERATION 44 ---------------- X=-0.5809024125060E+00 0.2848118988724E+00 0.1303728919462E+01 0.1674118977783E+01 F(X)= 0.3137376548868E+01 G(X)=-0.1539228552088E+02-0.1162638320384E+02 0.1261801082163E+02-0.5149736893622E+01 ----------- ITERATION 45 ---------------- X=-0.4885893648635E+00 0.2346511832659E+00 0.1316966410812E+01 0.1735836590727E+01 F(X)= 0.2552272260416E+01 G(X)=-0.3772283383262E+01-0.1704154970554E+01-0.4689023492564E-01-0.3151054465950E-01 ----------- ITERATION 46 ---------------- X=-0.4407669906764E+00 0.1889855647924E+00 0.1333855759465E+01 0.1780365948333E+01 F(X)= 0.2452838494288E+01 G(X)=-0.3814190039703E+01-0.1989237400696E+01 0.9402845550183E-01-0.7963101219392E-01 ----------- ITERATION 47 ---------------- X=-0.3839374961872E+00 0.1404152609539E+00 0.1351581924025E+01 0.1827739178100E+01 F(X)= 0.2338714175393E+01 G(X)=-0.3841783168525E+01-0.2372920661942E+01 0.2334184587396E+00-0.1256541413196E+00 ----------- ITERATION 48 ---------------- X=-0.3121786492297E+00 0.8777124997029E-01 0.1370533240494E+01 0.1879123379751E+01 F(X)= 0.2200392338085E+01 G(X)=-0.3833643427519E+01-0.2957226148317E+01 0.3651230713308E+00-0.1666681078916E+00 ----------- ITERATION 49 ---------------- X=-0.2108726860497E+00 0.2941617702453E-01 0.1391239964958E+01 0.1936188485323E+01 F(X)= 0.2017380273608E+01 G(X)=-0.3691292106408E+01-0.4079480300642E+01 0.4620464092225E+00-0.1913740775218E+00 ----------- ITERATION 50 ---------------- X= 0.3190783272781E-01-0.5228572138646E-01 0.1419775417400E+01 0.2016613142533E+01 F(X)= 0.1838368430719E+01 G(X)=-0.1255860071525E+01-0.1178799445180E+02 0.4046687481961E+00-0.1465023131455E+00 ----------- ITERATION 51 ---------------- X= 0.1178310613226E+00 0.1179934749586E-01 0.1397704116897E+01 0.1954657587299E+01 F(X)= 0.1325651560520E+01 G(X)=-0.1666075403460E+01-0.1476391643433E+01 0.2516150494056E+00-0.8774139500156E-01 ----------- ITERATION 52 ---------------- X= 0.1874648975354E+00 0.2886477970278E-01 0.1391766887314E+01 0.1938091189018E+01 F(X)= 0.1193174153393E+01 G(X)=-0.1154284749958E+01-0.2298383926031E+01 0.2443893674169E+00-0.8532757149964E-01 ----------- ITERATION 53 ---------------- X= 0.2420954003406E+00 0.5069249046199E-01 0.1384141685979E+01 0.1916902718889E+01 F(X)= 0.1087128048236E+01 G(X)=-0.7490735634423E+00-0.2604872834691E+01 0.2428597108188E+00-0.8503562499704E-01 ----------- ITERATION 54 ---------------- X= 0.3159670901037E+00 0.8886604333971E-01 0.1370739289565E+01 0.1879853590247E+01 F(X)= 0.9480065441772E+00 G(X)= 0.1829284670467E-01-0.3177632959248E+01 0.2838722378051E+00-0.1004735961367E+00 ----------- ITERATION 55 ---------------- X= 0.4325149249586E+00 0.1684705245474E+00 0.1342521417616E+01 0.1802627505599E+01 F(X)= 0.7493575603471E+00 G(X)= 0.2082707359336E+01-0.4624594437267E+01 0.5575985539654E+00-0.2037275386769E+00 ----------- ITERATION 56 ---------------- X= 0.6183338774436E+00 0.3513954543615E+00 0.1276218834503E+01 0.1625490152295E+01 F(X)= 0.4863242725035E+00 G(X)= 0.6889501807836E+01-0.6905369258433E+01 0.2043046791477E+01-0.7914488021691E+00 ----------- ITERATION 57 ---------------- X= 0.7168898446265E+00 0.5068648759514E+00 0.1216217220437E+01 0.1476369027545E+01 F(X)= 0.2294125091616E+00 G(X)= 0.1460053432197E+01-0.1942453971873E+01 0.1665100756281E+01-0.6481703558950E+00 ----------- ITERATION 58 ---------------- X= 0.8990421298967E+00 0.7761976557296E+00 0.1108425434069E+01 0.1217227109556E+01 F(X)= 0.1563934933716E+00 G(X)= 0.1133427810266E+02-0.6635526229401E+01 0.4757797307166E+01-0.2091664957464E+01 ----------- ITERATION 59 ---------------- X= 0.9252550978672E+00 0.8559306593130E+00 0.1069258226315E+01 0.1141721909416E+01 F(X)= 0.1883692279475E-01 G(X)=-0.8791734906131E-01-0.1373707579593E+00 0.7510531106887E+00-0.2762108848727E+00 ----------- ITERATION 60 ---------------- X= 0.9967833282706E+00 0.9884597968375E+00 0.1007781411738E+01 0.1011789658507E+01 F(X)= 0.4067291887128E-02 G(X)= 0.2033877813148E+01-0.1023114720527E+01 0.1406451472003E+01-0.6804104798721E+00 ----------- ITERATION 61 ---------------- X= 0.9983854337581E+00 0.9967711196050E+00 0.1001670421173E+01 0.1003291955310E+01 F(X)= 0.9929949105813E-05 G(X)=-0.2276119169655E-02-0.5101342463444E-03 0.2198723560537E-01-0.6733080553696E-02 ----------- ITERATION 62 ---------------- X= 0.9999896131831E+00 0.9999767166964E+00 0.1000012954723E+01 0.1000023097305E+01 F(X)= 0.1832888383061E-08 G(X)= 0.9958075229858E-03-0.5114699822275E-03 0.1049769896307E-02-0.4974938710699E-03 OPTSTP RELATIVE GRADIENT CLOSE TO ZERO. OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION. FINAL X= 0.9999995554147E+00 0.9999991407188E+00 0.1000000183507E+01 0.1000000398894E+01 GRADIENT G=-0.1640905875818E-06 0.1150812112878E-09 0.3063518809213E-06-0.5242027244489E-07 FUNCTION VALUE F(X)= 0.2689995237610E-11 AT ITERATION 63