# to unbundle, sh this file (in an empty directory) echo README 1>&2 sed >README <<'//GO.SYSIN DD README' 's/^-//' -tn.f Subroutines TN and TNBC, which solve unconstrained and simply - bounded optimization problems by a truncated Newton algorithm. - You must adapt subroutine MCHPR1 to your machine. Comments in - TN and TNBC explain their usage. - -blas.f Standard BLAS (basic linear algebra subroutines) used by tn.f, - plus one nonstandard one (DXPY). - - - Sample calling programs: - -mainc.f Customized, no bounds. -maincb.f Customized, bounds on variables. -mainfc.f Finite-differencing, customized, no bounds. -mains.f Easy to use, no bounds. -mainsb.f Easy to use, bounds on variables. - -Output from running these sample programs on an SGI 4D/240S computer - (IEEE arithmetic) appears in the corresponding .out files: - mainc.out is from mainc.f, maincb.out from maincb.f, etc. - -All the Fortran files were written and supplied by - - Stephen G. Nash - George Mason University - Dept. of Operational Research & Applied Statistics - Fairfax, VA 22030 - - SNASH%GMUVAX.BITNET@CUNYVM.CUNY.EDU //GO.SYSIN DD README echo tn.f 1>&2 sed >tn.f <<'//GO.SYSIN DD tn.f' 's/^-//' -C%% TRUNCATED-NEWTON METHOD: SUBROUTINES -C FOR OTHER MACHINES, MODIFY ROUTINE MCHPR1 (MACHINE EPSILON) -C WRITTEN BY: STEPHEN G. NASH -C OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. -C GEORGE MASON UNIVERSITY -C FAIRFAX, VA 22030 -C****************************************************************** - SUBROUTINE TN (IERROR, N, X, F, G, W, LW, SFUN) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER IERROR, N, LW - DOUBLE PRECISION X(N), G(N), F, W(LW) -C -C THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM -C -C MINIMIZE F(X) -C X -C -C WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS -C A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA -C THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), -C PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES -C NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A -C GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW. -C IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS -C ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. -C -C SUBROUTINE PARAMETERS: -C -C IERROR - (INTEGER) ERROR CODE -C ( 0 => NORMAL RETURN) -C ( 2 => MORE THAN MAXFUN EVALUATIONS) -C ( 3 => LINE SEARCH FAILED TO FIND -C ( LOWER POINT (MAY NOT BE SERIOUS) -C (-1 => ERROR IN INPUT PARAMETERS) -C N - (INTEGER) NUMBER OF VARIABLES -C X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL -C ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. -C G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL -C VALUE OF THE GRADIENT -C F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE -C OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE -C OF THE OBJECTIVE FUNCTION AT THE SOLUTION -C W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N -C LW - (INTEGER) THE DECLARED DIMENSION OF W -C SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION -C AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE -C THE CALLING SEQUENCE -C SUBROUTINE SFUN (N, X, F, G) -C INTEGER N -C DOUBLE PRECISION X(N), G(N), F -C -C THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE -C LMQN. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE -C OF THIS ALGORITHM SHOULD CALL LMQN DIRECTLY. -C -C---------------------------------------------------------------------- -C THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON -C ALGORITHM. THE PARAMETERS ARE: -C -C ETA - SEVERITY OF THE LINESEARCH -C MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS -C XTOL - DESIRED ACCURACY FOR THE SOLUTION X* -C STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH -C ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES -C MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT -C 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. -C MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP -C - DOUBLE PRECISION ETA, ACCRCY, XTOL, STEPMX, DSQRT, MCHPR1 - EXTERNAL SFUN -C -C SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE -C - MAXIT = N/2 - IF (MAXIT .GT. 50) MAXIT = 50 - IF (MAXIT .LE. 0) MAXIT = 1 - MSGLVL = 1 - MAXFUN = 150*N - ETA = .25D0 - STEPMX = 1.D1 - ACCRCY = 1.D2*MCHPR1() - XTOL = DSQRT(ACCRCY) -C -C MINIMIZE THE FUNCTION -C - CALL LMQN (IERROR, N, X, F, G, W, LW, SFUN, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) -C -C PRINT THE RESULTS -C - IF (IERROR .NE. 0) WRITE(*,800) IERROR - WRITE(*,810) F - IF (MSGLVL .LT. 1) RETURN - WRITE(*,820) - NMAX = 10 - IF (N .LT. NMAX) NMAX = N - WRITE(*,830) (I,X(I),I=1,NMAX) - RETURN -800 FORMAT(//,' ERROR CODE =', I3) -810 FORMAT(//,' OPTIMAL FUNCTION VALUE = ', 1PD22.15) -820 FORMAT(10X, 'CURRENT SOLUTION IS (AT MOST 10 COMPONENTS)', /, - * 14X, 'I', 11X, 'X(I)') -830 FORMAT(10X, I5, 2X, 1PD22.15) - END -C -C - SUBROUTINE TNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER IERROR, N, LW, IPIVOT(N) - DOUBLE PRECISION X(N), G(N), F, W(LW), LOW(N), UP(N) -C -C THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM -C -C MINIMIZE F(X) -C X -C SUBJECT TO LOW <= X <= UP -C -C WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS -C A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA -C THE LANCZOS ALGORITHM" BY S.G. NASH (TECHNICAL REPORT 378, MATH. -C THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), -C PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES -C NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A -C GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW. -C IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS -C ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. -C -C SUBROUTINE PARAMETERS: -C -C IERROR - (INTEGER) ERROR CODE -C ( 0 => NORMAL RETURN -C ( 2 => MORE THAN MAXFUN EVALUATIONS -C ( 3 => LINE SEARCH FAILED TO FIND LOWER -C ( POINT (MAY NOT BE SERIOUS) -C (-1 => ERROR IN INPUT PARAMETERS -C N - (INTEGER) NUMBER OF VARIABLES -C X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL -C ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. -C G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL -C VALUE OF THE GRADIENT -C F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE -C OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE -C OF THE OBJECTIVE FUNCTION AT THE SOLUTION -C W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N -C LW - (INTEGER) THE DECLARED DIMENSION OF W -C SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION -C AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE -C THE CALLING SEQUENCE -C SUBROUTINE SFUN (N, X, F, G) -C INTEGER N -C DOUBLE PRECISION X(N), G(N), F -C LOW, UP - (REAL*8) VECTORS OF LENGTH AT LEAST N CONTAINING -C THE LOWER AND UPPER BOUNDS ON THE VARIABLES. IF -C THERE ARE NO BOUNDS ON A PARTICULAR VARIABLE, SET -C THE BOUNDS TO -1.D38 AND 1.D38, RESPECTIVELY. -C IPIVOT - (INTEGER) WORK VECTOR OF LENGTH AT LEAST N, USED -C TO RECORD WHICH VARIABLES ARE AT THEIR BOUNDS. -C -C THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE -C LMQNBC. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE -C OF THIS ALGORITHM SHOULD CALL LMQBC DIRECTLY. -C -C---------------------------------------------------------------------- -C THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON -C ALGORITHM. THE PARAMETERS ARE: -C -C ETA - SEVERITY OF THE LINESEARCH -C MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS -C XTOL - DESIRED ACCURACY FOR THE SOLUTION X* -C STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH -C ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES -C MSGLVL - CONTROLS QUANTITY OF PRINTED OUTPUT -C 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. -C MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP -C - DOUBLE PRECISION ETA, ACCRCY, XTOL, STEPMX, DSQRT, MCHPR1 - EXTERNAL SFUN -C -C SET PARAMETERS FOR THE OPTIMIZATION ROUTINE -C - MAXIT = N/2 - IF (MAXIT .GT. 50) MAXIT = 50 - IF (MAXIT .LE. 0) MAXIT = 1 - MSGLVL = 1 - MAXFUN = 150*N - ETA = .25D0 - STEPMX = 1.D1 - ACCRCY = 1.D2*MCHPR1() - XTOL = DSQRT(ACCRCY) -C -C MINIMIZE FUNCTION -C - CALL LMQNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) -C -C PRINT RESULTS -C - IF (IERROR .NE. 0) WRITE(*,800) IERROR - WRITE(*,810) F - IF (MSGLVL .LT. 1) RETURN - WRITE(*,820) - NMAX = 10 - IF (N .LT. NMAX) NMAX = N - WRITE(*,830) (I,X(I),I=1,NMAX) - RETURN -800 FORMAT(//,' ERROR CODE =', I3) -810 FORMAT(//,' OPTIMAL FUNCTION VALUE = ', 1PD22.15) -820 FORMAT(10X, 'CURRENT SOLUTION IS (AT MOST 10 COMPONENTS)', /, - * 14X, 'I', 11X, 'X(I)') -830 FORMAT(10X, I5, 2X, 1PD22.15) - END -C -C - SUBROUTINE LMQN (IFAIL, N, X, F, G, W, LW, SFUN, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER MSGLVL, N, MAXFUN, IFAIL, LW - DOUBLE PRECISION X(N), G(N), W(LW), ETA, XTOL, STEPMX, F, ACCRCY -C -C THIS ROUTINE IS A TRUNCATED-NEWTON METHOD. -C THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY -C QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED -C IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3). -C FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TN. -C - INTEGER I, ICYCLE, IOLDG, IPK, IYK, LOLDG, LPK, LSR, - * LWTEST, LYK, LYR, NFTOTL, NITER, NM1, NUMF, NWHY - DOUBLE PRECISION ABSTOL, ALPHA, DIFNEW, DIFOLD, EPSMCH, - * EPSRED, FKEEP, FM, FNEW, FOLD, FSTOP, FTEST, GNORM, GSK, - * GTG, GTPNEW, OLDF, OLDGTP, ONE, PE, PEPS, PNORM, RELTOL, - * RTEPS, RTLEPS, RTOL, RTOLSQ, SMALL, SPE, TINY, - * TNYTOL, TOLEPS, XNORM, YKSK, YRSR, ZERO - LOGICAL LRESET, UPD1 -C -C THE FOLLOWING IMSL AND STANDARD FUNCTIONS ARE USED -C - DOUBLE PRECISION DABS, DDOT, DSQRT, STEP1, DNRM2 - EXTERNAL SFUN - COMMON /SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, - * LOLDG,LHG,LHYK,LPK,LEMAT,LWTEST -C -C INITIALIZE PARAMETERS AND CONSTANTS -C - IF (MSGLVL .GE. -2) WRITE(*,800) - CALL SETPAR(N) - UPD1 = .TRUE. - IRESET = 0 - NFEVAL = 0 - NMODIF = 0 - NLINCG = 0 - FSTOP = F - ZERO = 0.D0 - ONE = 1.D0 - NM1 = N - 1 -C -C WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) -C - LHYR = LOLDG -C -C CHECK PARAMETERS AND SET CONSTANTS -C - CALL CHKUCP(LWTEST,MAXFUN,NWHY,N,ALPHA,EPSMCH, - * ETA,PEPS,RTEPS,RTOL,RTOLSQ,STEPMX,FTEST, - * XTOL,XNORM,X,LW,SMALL,TINY,ACCRCY) - IF (NWHY .LT. 0) GO TO 120 - CALL SETUCR(SMALL,NFTOTL,NITER,N,F,FNEW, - * FM,GTG,OLDF,SFUN,G,X) - FOLD = FNEW - IF (MSGLVL .GE. 1) WRITE(*,810) NITER,NFTOTL,NLINCG,FNEW,GTG -C -C CHECK FOR SMALL GRADIENT AT THE STARTING POINT. -C - FTEST = ONE + DABS(FNEW) - IF (GTG .LT. 1.D-4*EPSMCH*FTEST*FTEST) GO TO 90 -C -C SET INITIAL VALUES TO OTHER PARAMETERS -C - ICYCLE = NM1 - TOLEPS = RTOL + RTEPS - RTLEPS = RTOLSQ + EPSMCH - GNORM = DSQRT(GTG) - DIFNEW = ZERO - EPSRED = 5.0D-2 - FKEEP = FNEW -C -C SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. -C - IDIAGB = LDIAGB - DO 10 I = 1,N - W(IDIAGB) = ONE - IDIAGB = IDIAGB + 1 -10 CONTINUE -C -C ..................START OF MAIN ITERATIVE LOOP.......... -C -C COMPUTE THE NEW SEARCH DIRECTION -C - MODET = MSGLVL - 3 - CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), - * W(LDIAGB),W(LEMAT),X,G,W(LZK), - * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, - * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,SFUN,.FALSE.,IPIVOT, - * ACCRCY,GTPNEW,GNORM,XNORM) -20 CONTINUE - CALL DCOPY(N,G,1,W(LOLDG),1) - PNORM = DNRM2(N,W(LPK),1) - OLDF = FNEW - OLDGTP = GTPNEW -C -C PREPARE TO COMPUTE THE STEP LENGTH -C - PE = PNORM + EPSMCH -C -C COMPUTE THE ABSOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH -C - RELTOL = RTEPS*(XNORM + ONE)/PE - ABSTOL = - EPSMCH*FTEST/(OLDGTP - EPSMCH) -C -C COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN -C THE LINEAR SEARCH -C - TNYTOL = EPSMCH*(XNORM + ONE)/PE - SPE = STEPMX/PE -C -C SET THE INITIAL STEP LENGTH. -C - ALPHA = STEP1(FNEW,FM,OLDGTP,SPE) -C -C PERFORM THE LINEAR SEARCH -C - CALL LINDER(N,SFUN,SMALL,EPSMCH,RELTOL,ABSTOL,TNYTOL, - * ETA,ZERO,SPE,W(LPK),OLDGTP,X,FNEW,ALPHA,G,NUMF, - * NWHY,W,LW) -C - FOLD = FNEW - NITER = NITER + 1 - NFTOTL = NFTOTL + NUMF - GTG = DDOT(N,G,1,G,1) - IF (MSGLVL .GE. 1) WRITE(*,810) NITER,NFTOTL,NLINCG,FNEW,GTG - IF (NWHY .LT. 0) GO TO 120 - IF (NWHY .EQ. 0 .OR. NWHY .EQ. 2) GO TO 30 -C -C THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT -C - NWHY = 3 - GO TO 100 -30 IF (NWHY .LE. 1) GO TO 40 - CALL SFUN(N,X,FNEW,G) - NFTOTL = NFTOTL + 1 -C -C TERMINATE IF MORE THAN MAXFUN EVALUTATIONS HAVE BEEN MADE -C -40 NWHY = 2 - IF (NFTOTL .GT. MAXFUN) GO TO 110 - NWHY = 0 -C -C SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS -C - DIFOLD = DIFNEW - DIFNEW = OLDF - FNEW -C -C IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE -C PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. -C - IF (ICYCLE .NE. 1) GO TO 50 - IF (DIFNEW .GT. 2.0D0 *DIFOLD) EPSRED = EPSRED + EPSRED - IF (DIFNEW .LT. 5.0D-1*DIFOLD) EPSRED = 5.0D-1*EPSRED -50 CONTINUE - GNORM = DSQRT(GTG) - FTEST = ONE + DABS(FNEW) - XNORM = DNRM2(N,X,1) -C -C TEST FOR CONVERGENCE -C - IF ((ALPHA*PNORM .LT. TOLEPS*(ONE + XNORM) - * .AND. DABS(DIFNEW) .LT. RTLEPS*FTEST - * .AND. GTG .LT. PEPS*FTEST*FTEST) - * .OR. GTG .LT. 1.D-4*ACCRCY*FTEST*FTEST) GO TO 90 -C -C COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE -C IN THE GRADIENTS -C - ISK = LSK - IPK = LPK - IYK = LYK - IOLDG = LOLDG - DO 60 I = 1,N - W(IYK) = G(I) - W(IOLDG) - W(ISK) = ALPHA*W(IPK) - IPK = IPK + 1 - ISK = ISK + 1 - IYK = IYK + 1 - IOLDG = IOLDG + 1 -60 CONTINUE -C -C SET UP PARAMETERS USED IN UPDATING THE DIRECTION OF SEARCH. -C - YKSK = DDOT(N,W(LYK),1,W(LSK),1) - LRESET = .FALSE. - IF (ICYCLE .EQ. NM1 .OR. DIFNEW .LT. - * EPSRED*(FKEEP-FNEW)) LRESET = .TRUE. - IF (LRESET) GO TO 70 - YRSR = DDOT(N,W(LYR),1,W(LSR),1) - IF (YRSR .LE. ZERO) LRESET = .TRUE. -70 CONTINUE - UPD1 = .FALSE. -C -C COMPUTE THE NEW SEARCH DIRECTION -C - MODET = MSGLVL - 3 - CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), - * W(LDIAGB),W(LEMAT),X,G,W(LZK), - * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, - * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,SFUN,.FALSE.,IPIVOT, - * ACCRCY,GTPNEW,GNORM,XNORM) - IF (LRESET) GO TO 80 -C -C STORE THE ACCUMULATED CHANGE IN THE POINT AND GRADIENT AS AN -C "AVERAGE" DIRECTION FOR PRECONDITIONING. -C - CALL DXPY(N,W(LSK),1,W(LSR),1) - CALL DXPY(N,W(LYK),1,W(LYR),1) - ICYCLE = ICYCLE + 1 - GOTO 20 -C -C RESET -C -80 IRESET = IRESET + 1 -C -C INITIALIZE THE SUM OF ALL THE CHANGES IN X. -C - CALL DCOPY(N,W(LSK),1,W(LSR),1) - CALL DCOPY(N,W(LYK),1,W(LYR),1) - FKEEP = FNEW - ICYCLE = 1 - GO TO 20 -C -C ...............END OF MAIN ITERATION....................... -C -90 IFAIL = 0 - F = FNEW - RETURN -100 OLDF = FNEW -C -C LOCAL SEARCH HERE COULD BE INSTALLED HERE -C -110 F = OLDF -C -C SET IFAIL -C -120 IFAIL = NWHY - RETURN -800 FORMAT(//' NIT NF CG', 9X, 'F', 21X, 'GTG',//) -810 FORMAT(' ',I3,1X,I4,1X,I4,1X,1PD22.15,2X,1PD15.8) - END -C -C - SUBROUTINE LMQNBC (IFAIL, N, X, F, G, W, LW, SFUN, LOW, UP, - * IPIVOT, MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER MSGLVL,N,MAXFUN,IFAIL,LW - INTEGER IPIVOT(N) - DOUBLE PRECISION ETA,XTOL,STEPMX,F,ACCRCY - DOUBLE PRECISION X(N),G(N),W(LW),LOW(N),UP(N) -C -C THIS ROUTINE IS A BOUNDS-CONSTRAINED TRUNCATED-NEWTON METHOD. -C THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY -C QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED -C IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3). -C FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TNBC. -C - INTEGER I, ICYCLE, IOLDG, IPK, IYK, LOLDG, LPK, LSR, - * LWTEST, LYK, LYR, NFTOTL, NITER, NM1, NUMF, NWHY - DOUBLE PRECISION ABSTOL, ALPHA, DIFNEW, DIFOLD, EPSMCH, EPSRED, - * FKEEP, FLAST, FM, FNEW, FOLD, FSTOP, FTEST, GNORM, GSK, - * GTG, GTPNEW, OLDF, OLDGTP, ONE, PE, PEPS, PNORM, RELTOL, - * RTEPS, RTLEPS, RTOL, RTOLSQ, SMALL, SPE, TINY, - * TNYTOL, TOLEPS, XNORM, YKSK, YRSR, ZERO - LOGICAL CONV, LRESET, UPD1, NEWCON -C -C THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE USED -C - DOUBLE PRECISION DABS, DDOT, DNRM2, DSQRT, STEP1 - EXTERNAL SFUN - COMMON/SUBSCR/ LGV, LZ1, LZK, LV, LSK, LYK, LDIAGB, LSR, LYR, - * LOLDG, LHG, LHYK, LPK, LEMAT, LWTEST -C -C CHECK THAT INITIAL X IS FEASIBLE AND THAT THE BOUNDS ARE CONSISTENT -C - CALL CRASH(N,X,IPIVOT,LOW,UP,IER) - IF (IER .NE. 0) WRITE(*,800) - IF (IER .NE. 0) RETURN - IF (MSGLVL .GE. 1) WRITE(*,810) -C -C INITIALIZE VARIABLES -C - CALL SETPAR(N) - UPD1 = .TRUE. - IRESET = 0 - NFEVAL = 0 - NMODIF = 0 - NLINCG = 0 - FSTOP = F - CONV = .FALSE. - ZERO = 0.D0 - ONE = 1.D0 - NM1 = N - 1 -C -C WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) -C - LHYR = LOLDG -C -C CHECK PARAMETERS AND SET CONSTANTS -C - CALL CHKUCP(LWTEST,MAXFUN,NWHY,N,ALPHA,EPSMCH, - * ETA,PEPS,RTEPS,RTOL,RTOLSQ,STEPMX,FTEST, - * XTOL,XNORM,X,LW,SMALL,TINY,ACCRCY) - IF (NWHY .LT. 0) GO TO 160 - CALL SETUCR(SMALL,NFTOTL,NITER,N,F,FNEW, - * FM,GTG,OLDF,SFUN,G,X) - FOLD = FNEW - FLAST = FNEW -C -C TEST THE LAGRANGE MULTIPLIERS TO SEE IF THEY ARE NON-NEGATIVE. -C BECAUSE THE CONSTRAINTS ARE ONLY LOWER BOUNDS, THE COMPONENTS -C OF THE GRADIENT CORRESPONDING TO THE ACTIVE CONSTRAINTS ARE THE -C LAGRANGE MULTIPLIERS. AFTERWORDS, THE PROJECTED GRADIENT IS FORMED. -C - DO 10 I = 1,N - IF (IPIVOT(I) .EQ. 2) GO TO 10 - IF (-IPIVOT(I)*G(I) .GE. 0.D0) GO TO 10 - IPIVOT(I) = 0 -10 CONTINUE - CALL ZTIME(N,G,IPIVOT) - GTG = DDOT(N,G,1,G,1) - IF (MSGLVL .GE. 1) - * CALL MONIT(N,X,FNEW,G,NITER,NFTOTL,NFEVAL,LRESET,IPIVOT) -C -C CHECK IF THE INITIAL POINT IS A LOCAL MINIMUM. -C - FTEST = ONE + DABS(FNEW) - IF (GTG .LT. 1.D-4*EPSMCH*FTEST*FTEST) GO TO 130 -C -C SET INITIAL VALUES TO OTHER PARAMETERS -C - ICYCLE = NM1 - TOLEPS = RTOL + RTEPS - RTLEPS = RTOLSQ + EPSMCH - GNORM = DSQRT(GTG) - DIFNEW = ZERO - EPSRED = 5.0D-2 - FKEEP = FNEW -C -C SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. -C - IDIAGB = LDIAGB - DO 15 I = 1,N - W(IDIAGB) = ONE - IDIAGB = IDIAGB + 1 -15 CONTINUE -C -C ..................START OF MAIN ITERATIVE LOOP.......... -C -C COMPUTE THE NEW SEARCH DIRECTION -C - MODET = MSGLVL - 3 - CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), - * W(LDIAGB),W(LEMAT),X,G,W(LZK), - * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, - * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,SFUN,.TRUE.,IPIVOT, - * ACCRCY,GTPNEW,GNORM,XNORM) -20 CONTINUE - CALL DCOPY(N,G,1,W(LOLDG),1) - PNORM = DNRM2(N,W(LPK),1) - OLDF = FNEW - OLDGTP = GTPNEW -C -C PREPARE TO COMPUTE THE STEP LENGTH -C - PE = PNORM + EPSMCH -C -C COMPUTE THE ABSOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH -C - RELTOL = RTEPS*(XNORM + ONE)/PE - ABSTOL = - EPSMCH*FTEST/(OLDGTP - EPSMCH) -C -C COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN -C THE LINEAR SEARCH -C - TNYTOL = EPSMCH*(XNORM + ONE)/PE - CALL STPMAX(STEPMX,PE,SPE,N,X,W(LPK),IPIVOT,LOW,UP) -C -C SET THE INITIAL STEP LENGTH. -C - ALPHA = STEP1(FNEW,FM,OLDGTP,SPE) -C -C PERFORM THE LINEAR SEARCH -C - CALL LINDER(N,SFUN,SMALL,EPSMCH,RELTOL,ABSTOL,TNYTOL, - * ETA,ZERO,SPE,W(LPK),OLDGTP,X,FNEW,ALPHA,G,NUMF, - * NWHY,W,LW) - NEWCON = .FALSE. - IF (DABS(ALPHA-SPE) .GT. 1.D1*EPSMCH) GO TO 30 - NEWCON = .TRUE. - NWHY = 0 - CALL MODZ(N,X,W(LPK),IPIVOT,EPSMCH,LOW,UP,FLAST,FNEW) - FLAST = FNEW -C -30 IF (MSGLVL .GE. 3) WRITE(*,820) ALPHA,PNORM - FOLD = FNEW - NITER = NITER + 1 - NFTOTL = NFTOTL + NUMF -C -C IF REQUIRED, PRINT THE DETAILS OF THIS ITERATION -C - IF (MSGLVL .GE. 1) - * CALL MONIT(N,X,FNEW,G,NITER,NFTOTL,NFEVAL,LRESET,IPIVOT) - IF (NWHY .LT. 0) GO TO 160 - IF (NWHY .EQ. 0 .OR. NWHY .EQ. 2) GO TO 40 -C -C THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT -C - NWHY = 3 - GO TO 140 -40 IF (NWHY .LE. 1) GO TO 50 - CALL SFUN(N,X,FNEW,G) - NFTOTL = NFTOTL + 1 -C -C TERMINATE IF MORE THAN MAXFUN EVALUATIONS HAVE BEEN MADE -C -50 NWHY = 2 - IF (NFTOTL .GT. MAXFUN) GO TO 150 - NWHY = 0 -C -C SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS -C - DIFOLD = DIFNEW - DIFNEW = OLDF - FNEW -C -C IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE -C PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. -C - IF (ICYCLE .NE. 1) GO TO 60 - IF (DIFNEW .GT. 2.D0*DIFOLD) EPSRED = EPSRED + EPSRED - IF (DIFNEW .LT. 5.0D-1*DIFOLD) EPSRED = 5.0D-1*EPSRED -60 CALL DCOPY(N,G,1,W(LGV),1) - CALL ZTIME(N,W(LGV),IPIVOT) - GTG = DDOT(N,W(LGV),1,W(LGV),1) - GNORM = DSQRT(GTG) - FTEST = ONE + DABS(FNEW) - XNORM = DNRM2(N,X,1) -C -C TEST FOR CONVERGENCE -C - CALL CNVTST(CONV,ALPHA,PNORM,TOLEPS,XNORM,DIFNEW,RTLEPS, - * FTEST,GTG,PEPS,EPSMCH,GTPNEW,FNEW,FLAST,G,IPIVOT,N,ACCRCY) - IF (CONV) GO TO 130 - CALL ZTIME(N,G,IPIVOT) -C -C COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE -C IN THE GRADIENTS -C - IF (NEWCON) GO TO 90 - ISK = LSK - IPK = LPK - IYK = LYK - IOLDG = LOLDG - DO 70 I = 1,N - W(IYK) = G(I) - W(IOLDG) - W(ISK) = ALPHA*W(IPK) - IPK = IPK + 1 - ISK = ISK + 1 - IYK = IYK + 1 - IOLDG = IOLDG + 1 -70 CONTINUE -C -C SET UP PARAMETERS USED IN UPDATING THE PRECONDITIONING STRATEGY. -C - YKSK = DDOT(N,W(LYK),1,W(LSK),1) - LRESET = .FALSE. - IF (ICYCLE .EQ. NM1 .OR. DIFNEW .LT. - * EPSRED*(FKEEP-FNEW)) LRESET = .TRUE. - IF (LRESET) GO TO 80 - YRSR = DDOT(N,W(LYR),1,W(LSR),1) - IF (YRSR .LE. ZERO) LRESET = .TRUE. -80 CONTINUE - UPD1 = .FALSE. -C -C COMPUTE THE NEW SEARCH DIRECTION -C -90 IF (UPD1 .AND. MSGLVL .GE. 3) WRITE(*,830) - IF (NEWCON .AND. MSGLVL .GE. 3) WRITE(*,840) - MODET = MSGLVL - 3 - CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), - * W(LDIAGB),W(LEMAT),X,G,W(LZK), - * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, - * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,SFUN,.TRUE.,IPIVOT, - * ACCRCY,GTPNEW,GNORM,XNORM) - IF (NEWCON) GO TO 20 - IF (LRESET) GO TO 110 -C -C COMPUTE THE ACCUMULATED STEP AND ITS CORRESPONDING -C GRADIENT DIFFERENCE. -C - CALL DXPY(N,W(LSK),1,W(LSR),1) - CALL DXPY(N,W(LYK),1,W(LYR),1) - ICYCLE = ICYCLE + 1 - GOTO 20 -C -C RESET -C -110 IRESET = IRESET + 1 -C -C INITIALIZE THE SUM OF ALL THE CHANGES IN X. -C - CALL DCOPY(N,W(LSK),1,W(LSR),1) - CALL DCOPY(N,W(LYK),1,W(LYR),1) - FKEEP = FNEW - ICYCLE = 1 - GO TO 20 -C -C ...............END OF MAIN ITERATION....................... -C -130 IFAIL = 0 - F = FNEW - RETURN -140 OLDF = FNEW -C -C LOCAL SEARCH COULD BE INSTALLED HERE -C -150 F = OLDF - IF (MSGLVL .GE. 1) CALL MONIT(N,X, - * F,G,NITER,NFTOTL,NFEVAL,IRESET,IPIVOT) -C -C SET IFAIL -C -160 IFAIL = NWHY - RETURN -800 FORMAT(' THERE IS NO FEASIBLE POINT; TERMINATING ALGORITHM') -810 FORMAT(//' NIT NF CG', 9X, 'F', 21X, 'GTG',//) -820 FORMAT(' LINESEARCH RESULTS: ALPHA,PNORM',2(1PD12.4)) -830 FORMAT(' UPD1 IS TRUE - TRIVIAL PRECONDITIONING') -840 FORMAT(' NEWCON IS TRUE - CONSTRAINT ADDED IN LINESEARCH') - END -C -C - SUBROUTINE MONIT(N,X,F,G,NITER,NFTOTL,NFEVAL,IRESET,IPIVOT) -C -C PRINT RESULTS OF CURRENT ITERATION -C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION X(N),F,G(N),GTG - INTEGER IPIVOT(N) -C - GTG = 0.D0 - DO 10 I = 1,N - IF (IPIVOT(I) .NE. 0) GO TO 10 - GTG = GTG + G(I)*G(I) -10 CONTINUE - WRITE(*,800) NITER,NFTOTL,NFEVAL,F,GTG - RETURN -800 FORMAT(' ',I4,1X,I4,1X,I4,1X,1PD22.15,2X,1PD15.8) - END -C -C - SUBROUTINE ZTIME(N,X,IPIVOT) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION X(N) - INTEGER IPIVOT(N) -C -C THIS ROUTINE MULTIPLIES THE VECTOR X BY THE CONSTRAINT MATRIX Z -C - DO 10 I = 1,N - IF (IPIVOT(I) .NE. 0) X(I) = 0.D0 -10 CONTINUE - RETURN - END -C -C - SUBROUTINE STPMAX(STEPMX,PE,SPE,N,X,P,IPIVOT,LOW,UP) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION LOW(N),UP(N),X(N),P(N),STEPMX,PE,SPE,T - INTEGER IPIVOT(N) -C -C COMPUTE THE MAXIMUM ALLOWABLE STEP LENGTH -C - SPE = STEPMX / PE -C SPE IS THE STANDARD (UNCONSTRAINED) MAX STEP - DO 10 I = 1,N - IF (IPIVOT(I) .NE. 0) GO TO 10 - IF (P(I) .EQ. 0.D0) GO TO 10 - IF (P(I) .GT. 0.D0) GO TO 5 - T = LOW(I) - X(I) - IF (T .GT. SPE*P(I)) SPE = T / P(I) - GO TO 10 -5 T = UP(I) - X(I) - IF (T .LT. SPE*P(I)) SPE = T / P(I) -10 CONTINUE - RETURN - END -C -C - SUBROUTINE MODZ(N,X,P,IPIVOT,EPSMCH,LOW,UP,FLAST,FNEW) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION X(N), P(N), EPSMCH, DABS, TOL, LOW(N), UP(N), - * FLAST, FNEW - INTEGER IPIVOT(N) -C -C UPDATE THE CONSTRAINT MATRIX IF A NEW CONSTRAINT IS ENCOUNTERED -C - DO 10 I = 1,N - IF (IPIVOT(I) .NE. 0) GO TO 10 - IF (P(I) .EQ. 0.D0) GO TO 10 - IF (P(I) .GT. 0.D0) GO TO 5 - TOL = 1.D1 * EPSMCH * (DABS(LOW(I)) + 1.D0) - IF (X(I)-LOW(I) .GT. TOL) GO TO 10 - FLAST = FNEW - IPIVOT(I) = -1 - X(I) = LOW(I) - GO TO 10 -5 TOL = 1.D1 * EPSMCH * (DABS(UP(I)) + 1.D0) - IF (UP(I)-X(I) .GT. TOL) GO TO 10 - FLAST = FNEW - IPIVOT(I) = 1 - X(I) = UP(I) -10 CONTINUE - RETURN - END -C -C - SUBROUTINE CNVTST(CONV,ALPHA,PNORM,TOLEPS,XNORM,DIFNEW,RTLEPS, - * FTEST,GTG,PEPS,EPSMCH,GTPNEW,FNEW,FLAST,G,IPIVOT,N,ACCRCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - LOGICAL CONV,LTEST - INTEGER IPIVOT(N) - DOUBLE PRECISION G(N), ALPHA, PNORM, TOLEPS, XNORM, DIFNEW, - * RTLEPS, FTEST, GTG, PEPS, EPSMCH, GTPNEW, FNEW, FLAST, ONE, - * CMAX, T, ACCRCY -C -C TEST FOR CONVERGENCE -C - IMAX = 0 - CMAX = 0.D0 - LTEST = FLAST - FNEW .LE. -5.D-1*GTPNEW - DO 10 I = 1,N - IF (IPIVOT(I) .EQ. 0 .OR. IPIVOT(I) .EQ. 2) GO TO 10 - T = -IPIVOT(I)*G(I) - IF (T .GE. 0.D0) GO TO 10 - CONV = .FALSE. - IF (LTEST) GO TO 10 - IF (CMAX .LE. T) GO TO 10 - CMAX = T - IMAX = I -10 CONTINUE - IF (IMAX .EQ. 0) GO TO 15 - IPIVOT(IMAX) = 0 - FLAST = FNEW - RETURN -15 CONTINUE - CONV = .FALSE. - ONE = 1.D0 - IF ((ALPHA*PNORM .GE. TOLEPS*(ONE + XNORM) - * .OR. DABS(DIFNEW) .GE. RTLEPS*FTEST - * .OR. GTG .GE. PEPS*FTEST*FTEST) - * .AND. GTG .GE. 1.D-4*ACCRCY*FTEST*FTEST) RETURN - CONV = .TRUE. -C -C FOR DETAILS, SEE GILL, MURRAY, AND WRIGHT (1981, P. 308) AND -C FLETCHER (1981, P. 116). THE MULTIPLIER TESTS (HERE, TESTING -C THE SIGN OF THE COMPONENTS OF THE GRADIENT) MAY STILL NEED TO -C MODIFIED TO INCORPORATE TOLERANCES FOR ZERO. -C - RETURN - END -C -C - SUBROUTINE CRASH(N,X,IPIVOT,LOW,UP,IER) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION X(N),LOW(N),UP(N) - INTEGER IPIVOT(N) -C -C THIS INITIALIZES THE CONSTRAINT INFORMATION, AND ENSURES THAT THE -C INITIAL POINT SATISFIES LOW <= X <= UP. -C THE CONSTRAINTS ARE CHECKED FOR CONSISTENCY. -C - IER = 0 - DO 30 I = 1,N - IF (X(I) .LT. LOW(I)) X(I) = LOW(I) - IF (X(I) .GT. UP(I)) X(I) = UP(I) - IPIVOT(I) = 0 - IF (X(I) .EQ. LOW(I)) IPIVOT(I) = -1 - IF (X(I) .EQ. UP(I)) IPIVOT(I) = 1 - IF (UP(I) .EQ. LOW(I)) IPIVOT(I) = 2 - IF (LOW(I) .GT. UP(I)) IER = -I -30 CONTINUE - RETURN - END -C -C THE VECTORS SK AND YK, ALTHOUGH NOT IN THE CALL, -C ARE USED (VIA THEIR POSITION IN W) BY THE ROUTINE MSOLVE. -C - SUBROUTINE MODLNP(MODET,ZSOL,GV,R,V,DIAGB,EMAT, - * X,G,ZK,N,W,LW,NITER,MAXIT,NFEVAL,NMODIF,NLINCG, - * UPD1,YKSK,GSK,YRSR,LRESET,SFUN,BOUNDS,IPIVOT,ACCRCY, - * GTP,GNORM,XNORM) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER MODET,N,NITER,IPIVOT(1) - DOUBLE PRECISION ZSOL(N),G(N),GV(N),R(N),V(N),DIAGB(N),W(LW) - DOUBLE PRECISION EMAT(N),ZK(N),X(N),ACCRCY - DOUBLE PRECISION ALPHA,BETA,DELTA,GSK,GTP,PR, - * QOLD,QNEW,QTEST,RHSNRM,RNORM,RZ,RZOLD,TOL,VGV,YKSK,YRSR - DOUBLE PRECISION GNORM,XNORM - DOUBLE PRECISION DDOT,DNRM2 - LOGICAL FIRST,UPD1,LRESET,BOUNDS - EXTERNAL SFUN -C -C THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT -C ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH -C DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM. WHEN THE VALUE OF THE -C QUADRATIC MODEL IS SUFFICIENTLY REDUCED, -C THE ITERATION IS TERMINATED. -C -C PARAMETERS -C -C MODET - INTEGER WHICH CONTROLS AMOUNT OF OUTPUT -C ZSOL - COMPUTED SEARCH DIRECTION -C G - CURRENT GRADIENT -C GV,GZ1,V - SCRATCH VECTORS -C R - RESIDUAL -C DIAGB,EMAT - DIAGONAL PRECONDITONING MATRIX -C NITER - NONLINEAR ITERATION # -C FEVAL - VALUE OF QUADRATIC FUNCTION -C -C ************************************************************* -C INITIALIZATION -C ************************************************************* -C -C GENERAL INITIALIZATION -C - IF (MODET .GT. 0) WRITE(*,800) - IF (MAXIT .EQ. 0) RETURN - FIRST = .TRUE. - RHSNRM = GNORM - TOL = 1.D-12 - QOLD = 0.D0 -C -C INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM -C - CALL INITPC(DIAGB,EMAT,N,W,LW,MODET, - * UPD1,YKSK,GSK,YRSR,LRESET) - DO 10 I = 1,N - R(I) = -G(I) - V(I) = 0.D0 - ZSOL(I) = 0.D0 -10 CONTINUE -C -C ************************************************************ -C MAIN ITERATION -C ************************************************************ -C - DO 30 K = 1,MAXIT - NLINCG = NLINCG + 1 - IF (MODET .GT. 1) WRITE(*,810) K -C -C CG ITERATION TO SOLVE SYSTEM OF EQUATIONS -C - IF (BOUNDS) CALL ZTIME(N,R,IPIVOT) - CALL MSOLVE(R,ZK,N,W,LW,UPD1,YKSK,GSK, - * YRSR,LRESET,FIRST) - IF (BOUNDS) CALL ZTIME(N,ZK,IPIVOT) - RZ = DDOT(N,R,1,ZK,1) - IF (RZ/RHSNRM .LT. TOL) GO TO 80 - IF (K .EQ. 1) BETA = 0.D0 - IF (K .GT. 1) BETA = RZ/RZOLD - DO 20 I = 1,N - V(I) = ZK(I) + BETA*V(I) -20 CONTINUE - IF (BOUNDS) CALL ZTIME(N,V,IPIVOT) - CALL GTIMS(V,GV,N,X,G,W,LW,SFUN,FIRST,DELTA,ACCRCY,XNORM) - IF (BOUNDS) CALL ZTIME(N,GV,IPIVOT) - NFEVAL = NFEVAL + 1 - VGV = DDOT(N,V,1,GV,1) - IF (VGV/RHSNRM .LT. TOL) GO TO 50 - CALL NDIA3(N,EMAT,V,GV,R,VGV,MODET) -C -C COMPUTE LINEAR STEP LENGTH -C - ALPHA = RZ / VGV - IF (MODET .GE. 1) WRITE(*,820) ALPHA -C -C COMPUTE CURRENT SOLUTION AND RELATED VECTORS -C - CALL DAXPY(N,ALPHA,V,1,ZSOL,1) - CALL DAXPY(N,-ALPHA,GV,1,R,1) -C -C TEST FOR CONVERGENCE -C - GTP = DDOT(N,ZSOL,1,G,1) - PR = DDOT(N,R,1,ZSOL,1) - QNEW = 5.D-1 * (GTP + PR) - QTEST = K * (1.D0 - QOLD/QNEW) - IF (QTEST .LT. 0.D0) GO TO 70 - QOLD = QNEW - IF (QTEST .LE. 5.D-1) GO TO 70 -C -C PERFORM CAUTIONARY TEST -C - IF (GTP .GT. 0) GO TO 40 - RZOLD = RZ -30 CONTINUE -C -C TERMINATE ALGORITHM -C - K = K-1 - GO TO 70 -C -C TRUNCATE ALGORITHM IN CASE OF AN EMERGENCY -C -40 IF (MODET .GE. -1) WRITE(*,830) K - CALL DAXPY(N,-ALPHA,V,1,ZSOL,1) - GTP = DDOT(N,ZSOL,1,G,1) - GO TO 90 -50 CONTINUE - IF (MODET .GT. -2) WRITE(*,840) -60 IF (K .GT. 1) GO TO 70 - CALL MSOLVE(G,ZSOL,N,W,LW,UPD1,YKSK,GSK,YRSR,LRESET,FIRST) - CALL NEGVEC(N,ZSOL) - IF (BOUNDS) CALL ZTIME(N,ZSOL,IPIVOT) - GTP = DDOT(N,ZSOL,1,G,1) -70 CONTINUE - IF (MODET .GE. -1) WRITE(*,850) K,RNORM - GO TO 90 -80 CONTINUE - IF (MODET .GE. -1) WRITE(*,860) - IF (K .GT. 1) GO TO 70 - CALL DCOPY(N,G,1,ZSOL,1) - CALL NEGVEC(N,ZSOL) - IF (BOUNDS) CALL ZTIME(N,ZSOL,IPIVOT) - GTP = DDOT(N,ZSOL,1,G,1) - GO TO 70 -C -C STORE (OR RESTORE) DIAGONAL PRECONDITIONING -C -90 CONTINUE - CALL DCOPY(N,EMAT,1,DIAGB,1) - RETURN -800 FORMAT(' ',//,' ENTERING MODLNP') -810 FORMAT(' ',//,' ### ITERATION ',I2,' ###') -820 FORMAT(' ALPHA',1PD16.8) -830 FORMAT(' G(T)Z POSITIVE AT ITERATION ',I2, - * ' - TRUNCATING METHOD',/) -840 FORMAT(' ',10X,'HESSIAN NOT POSITIVE-DEFINITE') -850 FORMAT(' ',/,8X,'MODLAN TRUNCATED AFTER ',I3,' ITERATIONS', - * ' RNORM = ',1PD14.6) -860 FORMAT(' PRECONDITIONING NOT POSITIVE-DEFINITE') - END -C -C - SUBROUTINE NDIA3(N,E,V,GV,R,VGV,MODET) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION E(N),V(N),GV(N),R(N),VGV,VR,DDOT -C -C UPDATE THE PRECONDITIOING MATRIX BASED ON A DIAGONAL VERSION -C OF THE BFGS QUASI-NEWTON UPDATE. -C - VR = DDOT(N,V,1,R,1) - DO 10 I = 1,N - E(I) = E(I) - R(I)*R(I)/VR + GV(I)*GV(I)/VGV - IF (E(I) .GT. 1.D-6) GO TO 10 - IF (MODET .GT. -2) WRITE(*,800) E(I) - E(I) = 1.D0 -10 CONTINUE - RETURN -800 FORMAT(' *** EMAT NEGATIVE: ',1PD16.8) - END -C -C SERVICE ROUTINES FOR OPTIMIZATION -C - SUBROUTINE NEGVEC(N,V) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER N - DOUBLE PRECISION V(N) -C -C NEGATIVE OF THE VECTOR V -C - INTEGER I - DO 10 I = 1,N - V(I) = -V(I) -10 CONTINUE - RETURN - END -C -C - SUBROUTINE LSOUT(ILOC,ITEST,XMIN,FMIN,GMIN,XW,FW,GW,U,A, - * B,TOL,EPS,SCXBD,XLAMDA) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION XMIN,FMIN,GMIN,XW,FW,GW,U,A,B, - * TOL,EPS,SCXBD,XLAMDA -C -C ERROR PRINTOUTS FOR GETPTC -C - DOUBLE PRECISION YA,YB,YBND,YW,YU - YU = XMIN + U - YA = A + XMIN - YB = B + XMIN - YW = XW + XMIN - YBND = SCXBD + XMIN - WRITE(*,800) - WRITE(*,810) TOL,EPS - WRITE(*,820) YA,YB - WRITE(*,830) YBND - WRITE(*,840) YW,FW,GW - WRITE(*,850) XMIN,FMIN,GMIN - WRITE(*,860) YU - WRITE(*,870) ILOC,ITEST - RETURN -800 FORMAT(///' OUTPUT FROM LINEAR SEARCH') -810 FORMAT(' TOL AND EPS'/2D25.14) -820 FORMAT(' CURRENT UPPER AND LOWER BOUNDS'/2D25.14) -830 FORMAT(' STRICT UPPER BOUND'/D25.14) -840 FORMAT(' XW, FW, GW'/3D25.14) -850 FORMAT(' XMIN, FMIN, GMIN'/3D25.14) -860 FORMAT(' NEW ESTIMATE'/2D25.14) -870 FORMAT(' ILOC AND ITEST'/2I3) - END -C -C - DOUBLE PRECISION FUNCTION STEP1(FNEW,FM,GTP,SMAX) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION FNEW,FM,GTP,SMAX -C -C ******************************************************** -C STEP1 RETURNS THE LENGTH OF THE INITIAL STEP TO BE TAKEN ALONG THE -C VECTOR P IN THE NEXT LINEAR SEARCH. -C ******************************************************** -C - DOUBLE PRECISION ALPHA,D,EPSMCH - DOUBLE PRECISION DABS,MCHPR1 - EPSMCH = MCHPR1() - D = DABS(FNEW-FM) - ALPHA = 1.D0 - IF (2.D0*D .LE. (-GTP) .AND. D .GE. EPSMCH) - * ALPHA = -2.D0*D/GTP - IF (ALPHA .GE. SMAX) ALPHA = SMAX - STEP1 = ALPHA - RETURN - END -C -C - DOUBLE PRECISION FUNCTION MCHPR1() - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION X -C -C RETURNS THE VALUE OF EPSMCH, WHERE EPSMCH IS THE SMALLEST POSSIBLE -C REAL NUMBER SUCH THAT 1.0 + EPSMCH .GT. 1.0 -C -C FOR VAX -C - MCHPR1 = 1.D-17 -C -C FOR SUN -C -C MCHPR1 = 1.0842021724855D-19 - RETURN - END -C -C - SUBROUTINE CHKUCP(LWTEST,MAXFUN,NWHY,N,ALPHA,EPSMCH, - * ETA,PEPS,RTEPS,RTOL,RTOLSQ,STEPMX,TEST, - * XTOL,XNORM,X,LW,SMALL,TINY,ACCRCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER LW,LWTEST,MAXFUN,NWHY,N - DOUBLE PRECISION ACCRCY,ALPHA,EPSMCH,ETA,PEPS,RTEPS,RTOL, - * RTOLSQ,STEPMX,TEST,XTOL,XNORM,SMALL,TINY - DOUBLE PRECISION X(N) -C -C CHECKS PARAMETERS AND SETS CONSTANTS WHICH ARE COMMON TO BOTH -C DERIVATIVE AND NON-DERIVATIVE ALGORITHMS -C - DOUBLE PRECISION DABS,DSQRT,MCHPR1 - EPSMCH = MCHPR1() - SMALL = EPSMCH*EPSMCH - TINY = SMALL - NWHY = -1 - RTEPS = DSQRT(EPSMCH) - RTOL = XTOL - IF (DABS(RTOL) .LT. ACCRCY) RTOL = 1.D1*RTEPS -C -C CHECK FOR ERRORS IN THE INPUT PARAMETERS -C - IF (LW .LT. LWTEST - * .OR. N .LT. 1 .OR. RTOL .LT. 0.D0 .OR. ETA .GE. 1.D0 .OR. - * ETA .LT. 0.D0 .OR. STEPMX .LT. RTOL .OR. - * MAXFUN .LT. 1) RETURN - NWHY = 0 -C -C SET CONSTANTS FOR LATER -C - RTOLSQ = RTOL*RTOL - PEPS = ACCRCY**0.6666D0 - XNORM = DNRM2(N,X,1) - ALPHA = 0.D0 - TEST = 0.D0 - RETURN - END -C -C - SUBROUTINE SETUCR(SMALL,NFTOTL,NITER,N,F,FNEW, - * FM,GTG,OLDF,SFUN,G,X) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER NFTOTL,NITER,N - DOUBLE PRECISION F,FNEW,FM,GTG,OLDF,SMALL - DOUBLE PRECISION G(N),X(N) - EXTERNAL SFUN -C -C CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET -C CONSTANTS FOR THE SUBSEQUENT MINIMIZATION -C - FM = F -C -C COMPUTE THE INITIAL FUNCTION VALUE -C - CALL SFUN(N,X,FNEW,G) - NFTOTL = 1 -C -C SET CONSTANTS FOR LATER -C - NITER = 0 - OLDF = FNEW - GTG = DDOT(N,G,1,G,1) - RETURN - END -C -C - SUBROUTINE GTIMS(V,GV,N,X,G,W,LW,SFUN,FIRST,DELTA,ACCRCY,XNORM) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION V(N),GV(N),DINV,DELTA,G(N) - DOUBLE PRECISION F,X(N),W(LW),ACCRCY,DSQRT,XNORM - LOGICAL FIRST - EXTERNAL SFUN - COMMON/SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, - * LHYR,LHG,LHYK,LPK,LEMAT,LWTEST -C -C THIS ROUTINE COMPUTES THE PRODUCT OF THE MATRIX G TIMES THE VECTOR -C V AND STORES THE RESULT IN THE VECTOR GV (FINITE-DIFFERENCE VERSION) -C - IF (.NOT. FIRST) GO TO 20 - DELTA = DSQRT(ACCRCY)*(1.D0+XNORM) - FIRST = .FALSE. -20 CONTINUE - DINV = 1.D0/DELTA - IHG = LHG - DO 30 I = 1,N - W(IHG) = X(I) + DELTA*V(I) - IHG = IHG + 1 -30 CONTINUE - CALL SFUN(N,W(LHG),F,GV) - DO 40 I = 1,N - GV(I) = (GV(I) - G(I))*DINV -40 CONTINUE - RETURN - END -C -C - SUBROUTINE MSOLVE(G,Y,N,W,LW,UPD1,YKSK,GSK, - * YRSR,LRESET,FIRST) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION G(N),Y(N),W(LW),YKSK,GSK,YRSR - LOGICAL UPD1,LRESET,FIRST -C -C THIS ROUTINE SETS UPT THE ARRAYS FOR MSLV -C - COMMON/SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, - * LHYR,LHG,LHYK,LPK,LEMAT,LWTEST - CALL MSLV(G,Y,N,W(LSK),W(LYK),W(LDIAGB),W(LSR),W(LYR),W(LHYR), - * W(LHG),W(LHYK),UPD1,YKSK,GSK,YRSR,LRESET,FIRST) - RETURN - END - SUBROUTINE MSLV(G,Y,N,SK,YK,DIAGB,SR,YR,HYR,HG,HYK, - * UPD1,YKSK,GSK,YRSR,LRESET,FIRST) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION G(N),Y(N) -C -C THIS ROUTINE ACTS AS A PRECONDITIONING STEP FOR THE -C LINEAR CONJUGATE-GRADIENT ROUTINE. IT IS ALSO THE -C METHOD OF COMPUTING THE SEARCH DIRECTION FROM THE -C GRADIENT FOR THE NON-LINEAR CONJUGATE-GRADIENT CODE. -C IT REPRESENTS A TWO-STEP SELF-SCALED BFGS FORMULA. -C - DOUBLE PRECISION DDOT,YKSK,GSK,YRSR,RDIAGB,YKHYK,GHYK, - * YKSR,YKHYR,YRHYR,GSR,GHYR - DOUBLE PRECISION SK(N),YK(N),DIAGB(N),SR(N),YR(N),HYR(N),HG(N), - * HYK(N),ONE - LOGICAL LRESET,UPD1,FIRST - IF (UPD1) GO TO 100 - ONE = 1.D0 - GSK = DDOT(N,G,1,SK,1) - IF (LRESET) GO TO 60 -C -C COMPUTE HG AND HY WHERE H IS THE INVERSE OF THE DIAGONALS -C - DO 57 I = 1,N - RDIAGB = 1.0D0/DIAGB(I) - HG(I) = G(I)*RDIAGB - IF (FIRST) HYK(I) = YK(I)*RDIAGB - IF (FIRST) HYR(I) = YR(I)*RDIAGB -57 CONTINUE - IF (FIRST) YKSR = DDOT(N,YK,1,SR,1) - IF (FIRST) YKHYR = DDOT(N,YK,1,HYR,1) - GSR = DDOT(N,G,1,SR,1) - GHYR = DDOT(N,G,1,HYR,1) - IF (FIRST) YRHYR = DDOT(N,YR,1,HYR,1) - CALL SSBFGS(N,ONE,SR,YR,HG,HYR,YRSR, - * YRHYR,GSR,GHYR,HG) - IF (FIRST) CALL SSBFGS(N,ONE,SR,YR,HYK,HYR,YRSR, - * YRHYR,YKSR,YKHYR,HYK) - YKHYK = DDOT(N,HYK,1,YK,1) - GHYK = DDOT(N,HYK,1,G,1) - CALL SSBFGS(N,ONE,SK,YK,HG,HYK,YKSK, - * YKHYK,GSK,GHYK,Y) - RETURN -60 CONTINUE -C -C COMPUTE GH AND HY WHERE H IS THE INVERSE OF THE DIAGONALS -C - DO 65 I = 1,N - RDIAGB = 1.D0/DIAGB(I) - HG(I) = G(I)*RDIAGB - IF (FIRST) HYK(I) = YK(I)*RDIAGB -65 CONTINUE - IF (FIRST) YKHYK = DDOT(N,YK,1,HYK,1) - GHYK = DDOT(N,G,1,HYK,1) - CALL SSBFGS(N,ONE,SK,YK,HG,HYK,YKSK, - * YKHYK,GSK,GHYK,Y) - RETURN -100 CONTINUE - DO 110 I = 1,N -110 Y(I) = G(I) / DIAGB(I) - RETURN - END -C -C - SUBROUTINE SSBFGS(N,GAMMA,SJ,YJ,HJV,HJYJ,YJSJ,YJHYJ, - * VSJ,VHYJ,HJP1V) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER N - DOUBLE PRECISION GAMMA,YJSJ,YJHYJ,VSJ,VHYJ - DOUBLE PRECISION SJ(N),YJ(N),HJV(N),HJYJ(N),HJP1V(N) -C -C SELF-SCALED BFGS -C - INTEGER I - DOUBLE PRECISION BETA,DELTA - DELTA = (1.D0 + GAMMA*YJHYJ/YJSJ)*VSJ/YJSJ - * - GAMMA*VHYJ/YJSJ - BETA = -GAMMA*VSJ/YJSJ - DO 10 I = 1,N - HJP1V(I) = GAMMA*HJV(I) + DELTA*SJ(I) + BETA*HJYJ(I) -10 CONTINUE - RETURN - END -C -C ROUTINES TO INITIALIZE PRECONDITIONER -C - SUBROUTINE INITPC(DIAGB,EMAT,N,W,LW,MODET, - * UPD1,YKSK,GSK,YRSR,LRESET) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION DIAGB(N),EMAT(N),W(LW) - DOUBLE PRECISION YKSK,GSK,YRSR - LOGICAL LRESET,UPD1 - COMMON/SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, - * LHYR,LHG,LHYK,LPK,LEMAT,LWTEST - CALL INITP3(DIAGB,EMAT,N,LRESET,YKSK,YRSR,W(LHYK), - * W(LSK),W(LYK),W(LSR),W(LYR),MODET,UPD1) - RETURN - END - SUBROUTINE INITP3(DIAGB,EMAT,N,LRESET,YKSK,YRSR,BSK, - * SK,YK,SR,YR,MODET,UPD1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DOUBLE PRECISION DIAGB(N),EMAT(N),YKSK,YRSR,BSK(N),SK(N), - * YK(N),COND,SR(N),YR(N),DDOT,SDS,SRDS,YRSK,TD,D1,DN - LOGICAL LRESET,UPD1 - IF (UPD1) GO TO 90 - IF (LRESET) GO TO 60 - DO 10 I = 1,N - BSK(I) = DIAGB(I)*SR(I) -10 CONTINUE - SDS = DDOT(N,SR,1,BSK,1) - SRDS = DDOT(N,SK,1,BSK,1) - YRSK = DDOT(N,YR,1,SK,1) - DO 20 I = 1,N - TD = DIAGB(I) - BSK(I) = TD*SK(I) - BSK(I)*SRDS/SDS+YR(I)*YRSK/YRSR - EMAT(I) = TD-TD*TD*SR(I)*SR(I)/SDS+YR(I)*YR(I)/YRSR -20 CONTINUE - SDS = DDOT(N,SK,1,BSK,1) - DO 30 I = 1,N - EMAT(I) = EMAT(I) - BSK(I)*BSK(I)/SDS+YK(I)*YK(I)/YKSK -30 CONTINUE - GO TO 110 -60 CONTINUE - DO 70 I = 1,N - BSK(I) = DIAGB(I)*SK(I) -70 CONTINUE - SDS = DDOT(N,SK,1,BSK,1) - DO 80 I = 1,N - TD = DIAGB(I) - EMAT(I) = TD - TD*TD*SK(I)*SK(I)/SDS + YK(I)*YK(I)/YKSK -80 CONTINUE - GO TO 110 -90 CONTINUE - CALL DCOPY(N,DIAGB,1,EMAT,1) -110 CONTINUE - IF (MODET .LT. 1) RETURN - D1 = EMAT(1) - DN = EMAT(1) - DO 120 I = 1,N - IF (EMAT(I) .LT. D1) D1 = EMAT(I) - IF (EMAT(I) .GT. DN) DN = EMAT(I) -120 CONTINUE - COND = DN/D1 - WRITE(*,800) D1,DN,COND -800 FORMAT(' ',//8X,'DMIN =',1PD12.4,' DMAX =',1PD12.4, - * ' COND =',1PD12.4,/) - RETURN - END -C -C - SUBROUTINE SETPAR(N) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER LSUB(14) - COMMON/SUBSCR/ LSUB,LWTEST -C -C SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE -C - DO 10 I = 1,14 - LSUB(I) = (I-1)*N + 1 -10 CONTINUE - LWTEST = LSUB(14) + N - 1 - RETURN - END -C -C LINE SEARCH ALGORITHMS OF GILL AND MURRAY -C - SUBROUTINE LINDER(N,SFUN,SMALL,EPSMCH,RELTOL,ABSTOL, - * TNYTOL,ETA,SFTBND,XBND,P,GTP,X,F,ALPHA,G,NFTOTL, - * IFLAG,W,LW) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER N,NFTOTL,IFLAG,LW - DOUBLE PRECISION SMALL,EPSMCH,RELTOL,ABSTOL,TNYTOL,ETA, - * SFTBND,XBND,GTP,F,ALPHA - DOUBLE PRECISION P(N),X(N),G(N),W(LW) -C -C - INTEGER I,IENTRY,ITEST,L,LG,LX,NUMF,ITCNT - DOUBLE PRECISION A,B,B1,BIG,E,FACTOR,FMIN,FPRESN,FU, - * FW,GMIN,GTEST1,GTEST2,GU,GW,OLDF,SCXBND,STEP, - * TOL,U,XMIN,XW,RMU,RTSMLL,UALPHA - LOGICAL BRAKTD -C -C THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE -C CALLED WITHIN LINDER -C - DOUBLE PRECISION DDOT,DSQRT - EXTERNAL SFUN -C -C ALLOCATE THE ADDRESSES FOR LOCAL WORKSPACE -C - LX = 1 - LG = LX + N - LSPRNT = 0 - NPRNT = 10000 - RTSMLL = DSQRT(SMALL) - BIG = 1.D0/SMALL - ITCNT = 0 -C -C SET THE ESTIMATED RELATIVE PRECISION IN F(X). -C - FPRESN = 10.D0*EPSMCH - NUMF = 0 - U = ALPHA - FU = F - FMIN = F - GU = GTP - RMU = 1.0D-4 -C -C FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY. -C - IENTRY = 1 -10 CONTINUE -C -C TEST FOR TOO MANY ITERATIONS -C - ITCNT = ITCNT + 1 - IFLAG = 1 - IF (ITCNT .GT. 20) GO TO 50 - IFLAG = 0 - CALL GETPTC(BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, - * FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, - * XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, - * BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST) -CLSOUT - IF (LSPRNT .GE. NPRNT) CALL LSOUT(IENTRY,ITEST,XMIN,FMIN,GMIN, - * XW,FW,GW,U,A,B,TOL,RELTOL,SCXBND,XBND) -C -C IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE -C CALCULATED. -C - IF (ITEST .NE. 1) GO TO 30 - UALPHA = XMIN + U - L = LX - DO 20 I = 1,N - W(L) = X(I) + UALPHA*P(I) - L = L + 1 -20 CONTINUE - CALL SFUN(N,W(LX),FU,W(LG)) - NUMF = NUMF + 1 - GU = DDOT(N,W(LG),1,P,1) -C -C THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS -C OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY -C LOWER THAN F AT THE ORIGIN. -C - IF (FU .LE. FMIN .AND. FU .LE. OLDF-UALPHA*GTEST1) - * CALL DCOPY(N,W(LG),1,G,1) - GOTO 10 -C -C IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND -C -30 CONTINUE - NFTOTL = NUMF - IFLAG = 1 - IF (ITEST .NE. 0) GO TO 50 -C -C IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE -C - IFLAG = 0 - F = FMIN - ALPHA = XMIN - DO 40 I = 1,N - X(I) = X(I) + ALPHA*P(I) -40 CONTINUE -50 RETURN - END -C -C - SUBROUTINE GETPTC(BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, - * FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, - * XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, - * BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - LOGICAL BRAKTD - INTEGER IENTRY,ITEST - DOUBLE PRECISION BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, - * FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, - * XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, - * GTEST1,GTEST2,TOL,DENOM -C -C ************************************************************ -C GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY -C ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC -C INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL -C IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A -C POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM. -C THE VALUE OF THE INTEGER PARAMETERS IENTRY DETERMINES THE PATH TAKEN -C THROUGH THE CODE. -C ************************************************************ -C - LOGICAL CONVRG - DOUBLE PRECISION ABGMIN,ABGW,ABSR,A1,CHORDM,CHORDU, - * D1,D2,P,Q,R,S,SCALE,SUMSQ,TWOTOL,XMIDPT - DOUBLE PRECISION ZERO, POINT1,HALF,ONE,THREE,FIVE,ELEVEN -C -C THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE CALLED -C WITHIN GETPTC -C - DOUBLE PRECISION DABS, DSQRT -C - ZERO = 0.D0 - POINT1 = 1.D-1 - HALF = 5.D-1 - ONE = 1.D0 - THREE = 3.D0 - FIVE = 5.D0 - ELEVEN = 11.D0 -C -C BRANCH TO APPROPRIATE SECTION OF CODE DEPENDING ON THE -C VALUE OF IENTRY. -C - GOTO (10,20), IENTRY -C -C IENTRY=1 -C CHECK INPUT PARAMETERS -C -10 ITEST = 2 - IF (U .LE. ZERO .OR. XBND .LE. TNYTOL .OR. GU .GT. ZERO) - * RETURN - ITEST = 1 - IF (XBND .LT. ABSTOL) ABSTOL = XBND - TOL = ABSTOL - TWOTOL = TOL + TOL -C -C A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS -C WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED. -C INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF -C FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH -C AT VALUES FOR LATEST ESTIMATE AT MINIMUM. -C - A = ZERO - XW = ZERO - XMIN = ZERO - OLDF = FU - FMIN = FU - FW = FU - GW = GU - GMIN = GU - STEP = U - FACTOR = FIVE -C -C THE MINIMUM HAS NOT YET BEEN BRACKETED. -C - BRAKTD = .FALSE. -C -C SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED -C EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.) SET THE UPPER BOUND -C ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND). -C - SCXBND = XBND - B = SCXBND + RELTOL*DABS(SCXBND) + ABSTOL - E = B + B - B1 = B -C -C COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA. -C - GTEST1 = -RMU*GU - GTEST2 = -ETA*GU -C -C SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION -C - IENTRY = 2 - GO TO 210 -C -C IENTRY = 2 -C -C UPDATE A,B,XW, AND XMIN -C -20 IF (FU .GT. FMIN) GO TO 60 -C -C IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT -C ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY. -C - CHORDU = OLDF - (XMIN + U)*GTEST1 - IF (FU .LE. CHORDU) GO TO 30 -C -C THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE -C CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND -C FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF -C UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES -C THE ROOT OF F(ALPHA)=CHORD(ALPHA). -C - CHORDM = OLDF - XMIN*GTEST1 - GU = -GMIN - DENOM = CHORDM-FMIN - IF (DABS(DENOM) .GE. 1.D-15) GO TO 25 - DENOM = 1.D-15 - IF (CHORDM-FMIN .LT. 0.D0) DENOM = -DENOM -25 CONTINUE - IF (XMIN .NE. ZERO) GU = GMIN*(CHORDU-FU)/DENOM - FU = HALF*U*(GMIN+GU) + FMIN - IF (FU .LT. FMIN) FU = FMIN - GO TO 60 -30 FW = FMIN - FMIN = FU - GW = GMIN - GMIN = GU - XMIN = XMIN + U - A = A-U - B = B-U - XW = -U - SCXBND = SCXBND - U - IF (GU .LE. ZERO) GO TO 40 - B = ZERO - BRAKTD = .TRUE. - GO TO 50 -40 A = ZERO -50 TOL = DABS(XMIN)*RELTOL + ABSTOL - GO TO 90 -C -C IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED -C BUT NEW POINT MAY NOW QUALIFY AS W. -C -60 IF (U .LT. ZERO) GO TO 70 - B = U - BRAKTD = .TRUE. - GO TO 80 -70 A = U -80 XW = U - FW = FU - GW = GU -90 TWOTOL = TOL + TOL - XMIDPT = HALF*(A + B) -C -C CHECK TERMINATION CRITERIA -C - CONVRG = DABS(XMIDPT) .LE. TWOTOL - HALF*(B-A) .OR. - * DABS(GMIN) .LE. GTEST2 .AND. FMIN .LT. OLDF .AND. - * (DABS(XMIN - XBND) .GT. TOL .OR. .NOT. BRAKTD) - IF (.NOT. CONVRG) GO TO 100 - ITEST = 0 - IF (XMIN .NE. ZERO) RETURN -C -C IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE -C CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA- -C UNIMODALITY CONSTANT, TOL. IF THE CHANGE IN F(X) IS LARGER THAN -C EXPECTED, REDUCE THE VALUE OF TOL. -C - ITEST = 3 - IF (DABS(OLDF-FW) .LE. FPRESN*(ONE + DABS(OLDF))) RETURN - TOL = POINT1*TOL - IF (TOL .LT. TNYTOL) RETURN - RELTOL = POINT1*RELTOL - ABSTOL = POINT1*ABSTOL - TWOTOL = POINT1*TWOTOL -C -C CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH -C -100 R = ZERO - Q = ZERO - S = ZERO - IF (DABS(E) .LE. TOL) GO TO 150 -C -C FIT CUBIC THROUGH XMIN AND XW -C - R = THREE*(FMIN-FW)/XW + GMIN + GW - ABSR = DABS(R) - Q = ABSR - IF (GW .EQ. ZERO .OR. GMIN .EQ. ZERO) GO TO 140 -C -C COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY -C WHICH AVOIDS UNDERFLOW AND OVERFLOW. -C - ABGW = DABS(GW) - ABGMIN = DABS(GMIN) - S = DSQRT(ABGMIN)*DSQRT(ABGW) - IF ((GW/ABGW)*GMIN .GT. ZERO) GO TO 130 -C -C COMPUTE THE SQUARE ROOT OF R*R + S*S. -C - SUMSQ = ONE - P = ZERO - IF (ABSR .GE. S) GO TO 110 -C -C THERE IS A POSSIBILITY OF OVERFLOW. -C - IF (S .GT. RTSMLL) P = S*RTSMLL - IF (ABSR .GE. P) SUMSQ = ONE +(ABSR/S)**2 - SCALE = S - GO TO 120 -C -C THERE IS A POSSIBILITY OF UNDERFLOW. -C -110 IF (ABSR .GT. RTSMLL) P = ABSR*RTSMLL - IF (S .GE. P) SUMSQ = ONE + (S/ABSR)**2 - SCALE = ABSR -120 SUMSQ = DSQRT(SUMSQ) - Q = BIG - IF (SCALE .LT. BIG/SUMSQ) Q = SCALE*SUMSQ - GO TO 140 -C -C COMPUTE THE SQUARE ROOT OF R*R - S*S -C -130 Q = DSQRT(DABS(R+S))*DSQRT(DABS(R-S)) - IF (R .GE. S .OR. R .LE. (-S)) GO TO 140 - R = ZERO - Q = ZERO - GO TO 150 -C -C COMPUTE THE MINIMUM OF FITTED CUBIC -C -140 IF (XW .LT. ZERO) Q = -Q - S = XW*(GMIN - R - Q) - Q = GW - GMIN + Q + Q - IF (Q .GT. ZERO) S = -S - IF (Q .LE. ZERO) Q = -Q - R = E - IF (B1 .NE. STEP .OR. BRAKTD) E = STEP -C -C CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH -C -150 A1 = A - B1 = B - STEP = XMIDPT - IF (BRAKTD) GO TO 160 - STEP = -FACTOR*XW - IF (STEP .GT. SCXBND) STEP = SCXBND - IF (STEP .NE. SCXBND) FACTOR = FIVE*FACTOR - GO TO 170 -C -C IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE -C WITHIN (A,B). -C -160 IF ((A .NE. ZERO .OR. XW .GE. ZERO) .AND. (B .NE. ZERO .OR. - * XW .LE. ZERO)) GO TO 180 -C -C IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE -C WITHIN (A1,B1). -C - D1 = XW - D2 = A - IF (A .EQ. ZERO) D2 = B -C THIS LINE MIGHT BE -C IF (A .EQ. ZERO) D2 = E - U = - D1/D2 - STEP = FIVE*D2*(POINT1 + ONE/U)/ELEVEN - IF (U .LT. ONE) STEP = HALF*D2*DSQRT(U) -170 IF (STEP .LE. ZERO) A1 = STEP - IF (STEP .GT. ZERO) B1 = STEP -C -C REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE -C REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED -C DURING THE LAST-BUT-ONE ITERATION. -C -180 IF (DABS(S) .LE. DABS(HALF*Q*R) .OR. - * S .LE. Q*A1 .OR. S .GE. Q*B1) GO TO 200 -C -C A CUBIC INTERPOLATION STEP -C - STEP = S/Q -C -C THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B. -C - IF (STEP - A .GE. TWOTOL .AND. B - STEP .GE. TWOTOL) GO TO 210 - IF (XMIDPT .GT. ZERO) GO TO 190 - STEP = -TOL - GO TO 210 -190 STEP = TOL - GO TO 210 -200 E = B-A -C -C IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO -C COMPUTE THE NEW POINT ON THE BOUNDARY). -C -210 IF (STEP .LT. SCXBND) GO TO 220 - STEP = SCXBND -C -C MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND. -C - SCXBND = SCXBND - (RELTOL*DABS(XBND)+ABSTOL)/(ONE + RELTOL) -220 U = STEP - IF (DABS(STEP) .LT. TOL .AND. STEP .LT. ZERO) U = -TOL - IF (DABS(STEP) .LT. TOL .AND. STEP .GE. ZERO) U = TOL - ITEST = 1 - RETURN - END //GO.SYSIN DD tn.f echo blas.f 1>&2 sed >blas.f <<'//GO.SYSIN DD blas.f' 's/^-//' -C%% TRUNCATED-NEWTON METHOD: BLAS -C NOTE: ALL ROUTINES HERE ARE FROM LINPACK WITH THE EXCEPTION -C OF DXPY (A VERSION OF DAXPY WITH A=1.0) -C WRITTEN BY: STEPHEN G. NASH -C OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. -C GEORGE MASON UNIVERSITY -C FAIRFAX, VA 22030 -C**************************************************************** - DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) -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(1),DY(1),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 - SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) -C -C CONSTANT TIMES A VECTOR PLUS A VECTOR. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1),DA - INTEGER I,INCX,INCY,IX,IY,M,MP1,N -C - IF(N.LE.0)RETURN - IF (DA .EQ. 0.0D0) 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) = DY(IY) + DA*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,4) - IF( M .EQ. 0 ) GO TO 40 - DO 30 I = 1,M - DY(I) = DY(I) + DA*DX(I) - 30 CONTINUE - IF( N .LT. 4 ) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,4 - DY(I) = DY(I) + DA*DX(I) - DY(I + 1) = DY(I + 1) + DA*DX(I + 1) - DY(I + 2) = DY(I + 2) + DA*DX(I + 2) - DY(I + 3) = DY(I + 3) + DA*DX(I + 3) - 50 CONTINUE - RETURN - END - DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - INTEGER NEXT - DOUBLE PRECISION DX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE - DATA ZERO, ONE /0.0D0, 1.0D0/ -C -C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE -C INCREMENT INCX . -C IF N .LE. 0 RETURN WITH RESULT = 0. -C IF N .GE. 1 THEN INCX MUST BE .GE. 1 -C -C C.L.LAWSON, 1978 JAN 08 -C -C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE -C HOPEFULLY APPLICABLE TO ALL MACHINES. -C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. -C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. -C WHERE -C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. -C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) -C V = LARGEST NO. (OVERFLOW LIMIT) -C -C BRIEF OUTLINE OF ALGORITHM.. -C -C PHASE 1 SCANS ZERO COMPONENTS. -C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO -C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO -C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M -C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. -C -C VALUES FOR CUTLO AND CUTHI.. -C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER -C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. -C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE -C UNIVAC AND DEC AT 2**(-103) -C THUS CUTLO = 2**(-51) = 4.44089E-16 -C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. -C THUS CUTHI = 2**(63.5) = 1.30438E19 -C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. -C THUS CUTLO = 2**(-33.5) = 8.23181D-11 -C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 -C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / -C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / - DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / -C - IF(N .GT. 0) GO TO 10 - DNRM2 = ZERO - GO TO 300 -C - 10 ASSIGN 30 TO NEXT - SUM = ZERO - NN = N * INCX -C BEGIN MAIN LOOP - I = 1 - 20 GO TO NEXT,(30, 50, 70, 110) - 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 - ASSIGN 50 TO NEXT - XMAX = ZERO -C -C PHASE 1. SUM IS ZERO -C - 50 IF( DX(I) .EQ. ZERO) GO TO 200 - IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 -C -C PREPARE FOR PHASE 2. - ASSIGN 70 TO NEXT - GO TO 105 -C -C PREPARE FOR PHASE 4. -C - 100 I = J - ASSIGN 110 TO NEXT - SUM = (SUM / DX(I)) / DX(I) - 105 XMAX = DABS(DX(I)) - GO TO 115 -C -C PHASE 2. SUM IS SMALL. -C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. -C - 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 -C -C COMMON CODE FOR PHASES 2 AND 4. -C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. -C - 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 - SUM = ONE + SUM * (XMAX / DX(I))**2 - XMAX = DABS(DX(I)) - GO TO 200 -C - 115 SUM = SUM + (DX(I)/XMAX)**2 - GO TO 200 -C -C -C PREPARE FOR PHASE 3. -C - 75 SUM = (SUM * XMAX) * XMAX -C -C -C FOR REAL OR D.P. SET HITEST = CUTHI/N -C FOR COMPLEX SET HITEST = CUTHI/(2*N) -C - 85 HITEST = CUTHI/FLOAT( N ) -C -C PHASE 3. SUM IS MID-RANGE. NO SCALING. -C - DO 95 J =I,NN,INCX - IF(DABS(DX(J)) .GE. HITEST) GO TO 100 - 95 SUM = SUM + DX(J)**2 - DNRM2 = DSQRT( SUM ) - GO TO 300 -C - 200 CONTINUE - I = I + INCX - IF ( I .LE. NN ) GO TO 20 -C -C END OF MAIN LOOP. -C -C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. -C - DNRM2 = XMAX * DSQRT(SUM) - 300 CONTINUE - RETURN - END - SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) -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(1),DY(1) - 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****************************************************************** -C SPECIAL BLAS FOR Y = X+Y -C****************************************************************** - SUBROUTINE DXPY(N,DX,INCX,DY,INCY) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) -C -C VECTOR PLUS A VECTOR. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C STEPHEN G. NASH 5/30/89. -C - DOUBLE PRECISION DX(1),DY(1) - 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) = 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,4) - IF( M .EQ. 0 ) GO TO 40 - DO 30 I = 1,M - DY(I) = DY(I) + DX(I) - 30 CONTINUE - IF( N .LT. 4 ) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,4 - DY(I) = DY(I) + DX(I) - DY(I + 1) = DY(I + 1) + DX(I + 1) - DY(I + 2) = DY(I + 2) + DX(I + 2) - DY(I + 3) = DY(I + 3) + DX(I + 3) - 50 CONTINUE - RETURN - END //GO.SYSIN DD blas.f echo mainc.f 1>&2 sed >mainc.f <<'//GO.SYSIN DD mainc.f' 's/^-//' -C*********************************************************************** -C CUSTOMIZED, NO BOUNDS -C*********************************************************************** -C MAIN PROGRAM TO MINIMIZE A FUNCTION (REPRESENTED BY THE ROUTINE SFUN) -C OF N VARIABLES X - CUSTOMIZED VERSION -C - DOUBLE PRECISION X(50), F, G(50), W(700) - DOUBLE PRECISION ETA, ACCRCY, XTOL, STEPMX, DSQRT - EXTERNAL SFUN -C -C SET UP FUNCTION AND VARIABLE INFORMATION -C N - NUMBER OF VARIABLES -C X - INITIAL ESTIMATE OF THE SOLUTION -C F - ROUGH ESTIMATE OF FUNCTION VALUE AT SOLUTION -C LW - DECLARED LENGTH OF THE ARRAY W -C - N = 10 - DO 10 I = 1,N - X(I) = I / FLOAT(N+1) -10 CONTINUE - F = 1.D0 - LW = 700 -C -C SET UP CUSTOMIZING PARAMETERS -C ETA - SEVERITY OF THE LINESEARCH -C MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS -C XTOL - DESIRED ACCURACY FOR THE SOLUTION X* -C STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH -C ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES -C MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT -C 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. -C MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP -C - MAXIT = N/2 - MAXFUN = 150*N - ETA = .25D0 - STEPMX = 1.D1 - ACCRCY = 1.D-15 - XTOL = DSQRT(ACCRCY) - MSGLVL = 1 -C -C MINIMIZE THE FUNCTION -C - CALL LMQN (IERROR, N, X, F, G, W, LW, SFUN, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) -C -C PRINT THE RESULTS -C - IF (IERROR .NE. 0) WRITE(*,800) IERROR - IF (MSGLVL .GE. 1) WRITE(*,810) - IF (MSGLVL .GE. 1) WRITE(*,820) (I,X(I),I=1,N) - STOP -800 FORMAT(//,' ERROR CODE =', I3,/) -810 FORMAT(10X, 'CURRENT SOLUTION IS ',/14X, 'I', 11X, 'X(I)') -820 FORMAT(10X, I5, 2X, 1PD22.15) - END -C -C - SUBROUTINE SFUN (N, X, F, G) - DOUBLE PRECISION X(N), G(N), F, T -C -C ROUTINE TO EVALUATE FUNCTION (F) AND GRADIENT (G) OF THE OBJECTIVE -C FUNCTION AT THE POINT X -C - F = 0.D0 - DO 10 I = 1,N - T = X(I) - I - F = F + T*T - G(I) = 2.D0 * T -10 CONTINUE - RETURN - END //GO.SYSIN DD mainc.f echo mainc.out 1>&2 sed >mainc.out <<'//GO.SYSIN DD mainc.out' 's/^-//' - - - NIT NF CG F GTG - - - 0 1 0 3.181818171522834D+02 1.27272727D+03 - 1 2 2 6.142878372311498D+01 2.45715135D+02 - 2 3 4 1.627901350786000D-02 6.51160540D-02 - 3 4 6 1.478295968387516D-19 5.91318387D-19 - 4 5 7 7.888609052210118D-31 3.15544362D-30 - CURRENT SOLUTION IS - I X(I) - 1 1.000000000000000D+00 - 2 2.000000000000000D+00 - 3 3.000000000000000D+00 - 4 4.000000000000000D+00 - 5 5.000000000000000D+00 - 6 6.000000000000000D+00 - 7 7.000000000000001D+00 - 8 8.000000000000000D+00 - 9 9.000000000000000D+00 - 10 1.000000000000000D+01 //GO.SYSIN DD mainc.out echo maincb.f 1>&2 sed >maincb.f <<'//GO.SYSIN DD maincb.f' 's/^-//' -C*********************************************************************** -C CUSTOMIZED, BOUNDS ON VARIABLES -C*********************************************************************** -C MAIN PROGRAM TO MINIMIZE A FUNCTION (REPRESENTED BY THE ROUTINE SFUN) -C SUBJECT TO BOUNDS ON THE VARIABLES X - CUSTOMIZED VERSION -C - DOUBLE PRECISION X(50), F, G(50), W(700), LOW(50), UP(50) - DOUBLE PRECISION ETA, ACCRCY, XTOL, STEPMX, DSQRT - INTEGER IPIVOT(50) - EXTERNAL SFUN -C -C SET UP FUNCTION AND VARIABLE INFORMATION -C N - NUMBER OF VARIABLES -C X - INITIAL ESTIMATE OF THE SOLUTION -C LOW - LOWER BOUNDS -C UP - UPPER BOUNDS -C F - ROUGH ESTIMATE OF FUNCTION VALUE AT SOLUTION -C LW - DECLARED LENGTH OF THE ARRAY W -C - N = 10 - DO 10 I = 1,N - X(I) = I / FLOAT(N+1) - LOW(I) = 0.D0 - UP(I) = 6.D0 -10 CONTINUE - F = 1.D0 - LW = 700 -C -C SET UP CUSTOMIZING PARAMETERS -C ETA - SEVERITY OF THE LINESEARCH -C MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS -C XTOL - DESIRED ACCURACY FOR THE SOLUTION X* -C STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH -C ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES -C MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT -C 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. -C MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP -C - MAXIT = N/2 - MAXFUN = 150*N - ETA = .25D0 - STEPMX = 1.D1 - ACCRCY = 1.D-15 - XTOL = DSQRT(ACCRCY) - MSGLVL = 1 -C -C MINIMIZE THE FUNCTION -C - CALL LMQNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) -C -C PRINT THE RESULTS -C - IF (IERROR .NE. 0) WRITE(*,800) IERROR - IF (MSGLVL .GE. 1) WRITE(*,810) - IF (MSGLVL .GE. 1) WRITE(*,820) (I,X(I),I=1,N) - STOP -800 FORMAT(//,' ERROR CODE =', I3,/) -810 FORMAT(10X, 'CURRENT SOLUTION IS ',/14X, 'I', 11X, 'X(I)') -820 FORMAT(10X, I5, 2X, 1PD22.15) - END -C -C - SUBROUTINE SFUN (N, X, F, G) - DOUBLE PRECISION X(N), G(N), F, T -C -C ROUTINE TO EVALUATE FUNCTION (F) AND GRADIENT (G) OF THE OBJECTIVE -C FUNCTION AT THE POINT X -C - F = 0.D0 - DO 10 I = 1,N - T = X(I) - I - F = F + T*T - G(I) = 2.D0 * T -10 CONTINUE - RETURN - END //GO.SYSIN DD maincb.f echo maincb.out 1>&2 sed >maincb.out <<'//GO.SYSIN DD maincb.out' 's/^-//' - - - NIT NF CG F GTG - - - 0 1 0 3.181818171522834D+02 1.27272727D+03 - 1 2 1 6.160000016784667D+01 1.82400001D+02 - 2 3 3 4.763276385151769D+01 9.05310554D+01 - 3 4 5 3.771160742015032D+01 3.48464297D+01 - 4 5 7 3.185626815692368D+01 7.42507263D+00 - 5 6 9 3.000000934245197D+01 3.73698079D-05 - 6 7 11 3.000000000007158D+01 2.86336350D-10 - 7 8 13 3.000000000000000D+01 3.32763631D-16 - 8 9 15 3.000000000000000D+01 3.32763631D-16 - 8 9 15 3.000000000000000D+01 3.32763631D-16 - - - ERROR CODE = 3 - - CURRENT SOLUTION IS - I X(I) - 1 9.999999933920257D-01 - 2 2.000000004649754D+00 - 3 2.999999996218138D+00 - 4 4.000000001695302D+00 - 5 5.000000000853727D+00 - 6 6.000000000000000D+00 - 7 6.000000000000000D+00 - 8 6.000000000000000D+00 - 9 6.000000000000000D+00 - 10 6.000000000000000D+00 //GO.SYSIN DD maincb.out echo mainfc.f 1>&2 sed >mainfc.f <<'//GO.SYSIN DD mainfc.f' 's/^-//' -C*********************************************************************** -C FINITE-DIFFERENCING, CUSTOMIZED, NO BOUNDS -C*********************************************************************** -C MAIN PROGRAM TO MINIMIZE A FUNCTION (REPRESENTED BY THE ROUTINE SFUN) -C OF N VARIABLES X - CUSTOMIZED VERSION, WITH FINITE DIFFERENCING -C - DOUBLE PRECISION X(50), F, G(50), W(700) - DOUBLE PRECISION ETA, ACCRCY, XTOL, STEPMX, FACCY, DSQRT - COMMON /FINDIF/ ACCRCY - EXTERNAL SFUN -C -C SET UP FUNCTION AND VARIABLE INFORMATION -C N - NUMBER OF VARIABLES -C X - INITIAL ESTIMATE OF THE SOLUTION -C F - ROUGH ESTIMATE OF FUNCTION VALUE AT SOLUTION -C LW - DECLARED LENGTH OF THE ARRAY W -C - N = 10 - DO 10 I = 1,N - X(I) = I / FLOAT(N+1) -10 CONTINUE - F = 1.D0 - LW = 700 -C -C SET UP CUSTOMIZING PARAMETERS -C ETA - SEVERITY OF THE LINESEARCH -C MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS -C XTOL - DESIRED ACCURACY FOR THE SOLUTION X* -C STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH -C FACCY - ACCURACY OF COMPUTED FUNCTION VALUE -C ACCRCY - ACCURACY OF COMPUTED FUNCTION AND GRADIENT VALUES -C (VALUES ADJUSTED BECAUSE OF FINITE DIFFERENCING OF GRADIENT) -C MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT -C 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. -C MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP -C - MAXIT = N/2 - MAXFUN = 150*N - ETA = .25D0 - STEPMX = 1.D1 - FACCY = 1.D-15 - ACCRCY = DSQRT(FACCY) - XTOL = DSQRT(ACCRCY) - MSGLVL = 1 -C -C MINIMIZE THE FUNCTION -C - CALL LMQN (IERROR, N, X, F, G, W, LW, SFUN, - * MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY, XTOL) -C -C PRINT THE RESULTS -C - IF (IERROR .NE. 0) WRITE(*,800) IERROR - IF (MSGLVL .GE. 1) WRITE(*,810) - IF (MSGLVL .GE. 1) WRITE(*,820) (I,X(I),I=1,N) - STOP -800 FORMAT(//,' ERROR CODE =', I3,/) -810 FORMAT(10X, 'CURRENT SOLUTION IS ',/14X, 'I', 11X, 'X(I)') -820 FORMAT(10X, I5, 2X, 1PD22.15) - END -C -C - SUBROUTINE SFUN (N, X, F, G) - DOUBLE PRECISION X(N), G(N), F, T, H, XH, FH, HINV - COMMON /FINDIF/ H -C -C ROUTINE TO COMPUTE FUNCTION (F) AND GRADIENT (G) OF THE OBJECTIVE -C FUNCTION AT THE POINT X -C GRADIENT OBTAINED VIA FINITE-DIFFERENCING -C FUNCTION OBTAINED FROM SUBROUTINE GETF -C - HINV = 1.D0 / H - CALL GETF (N, X, F) - DO 10 I = 1,N - XH = X(I) - X(I) = X(I) + H - CALL GETF (N, X, FH) - G(I) = (FH - F) * HINV - X(I) = XH -10 CONTINUE - RETURN - END -C -C - SUBROUTINE GETF (N, X, F) - DOUBLE PRECISION X(N), F, T -C -C ROUTINE TO EVALUATE FUNCTION (F) -C - F = 0.D0 - DO 10 I = 1,N - T = X(I) - I - F = F + T*T -10 CONTINUE - RETURN - END //GO.SYSIN DD mainfc.f echo mainfc.out 1>&2 sed >mainfc.out <<'//GO.SYSIN DD mainfc.out' 's/^-//' - - - NIT NF CG F GTG - - - 0 1 0 3.181818171522834D+02 1.27272745D+03 - 1 2 2 6.142887666611016D+01 2.45715506D+02 - 2 3 4 1.628235471185229D-02 6.51293747D-02 - 3 4 6 2.214977744825822D-08 8.86127752D-08 - 4 5 8 2.188893670127930D-12 8.77273472D-12 - CURRENT SOLUTION IS - I X(I) - 1 1.000000190659504D+00 - 2 1.999999534054855D+00 - 3 3.000000404161502D+00 - 4 3.999999966635100D+00 - 5 4.999999946566544D+00 - 6 6.000000036153342D+00 - 7 7.000000429448899D+00 - 8 8.000000284593261D+00 - 9 8.999998847693560D+00 - 10 1.000000041663653D+01 //GO.SYSIN DD mainfc.out echo mains.f 1>&2 sed >mains.f <<'//GO.SYSIN DD mains.f' 's/^-//' -C*********************************************************************** -C EASY TO USE, NO BOUNDS -C*********************************************************************** -C MAIN PROGRAM TO MINIMIZE A FUNCTION (REPRESENTED BY THE ROUTINE SFUN) -C OF N VARIABLES X -C - DOUBLE PRECISION X(50), F, G(50), W(700) - EXTERNAL SFUN -C -C DEFINE SUBROUTINE PARAMETERS -C N - NUMBER OF VARIABLES -C X - INITIAL ESTIMATE OF THE SOLUTION -C F - ROUGH ESTIMATE OF FUNCTION VALUE AT SOLUTION -C LW - DECLARED LENGTH OF THE ARRAY W -C - N = 10 - DO 10 I = 1,N - X(I) = I / FLOAT(N+1) -10 CONTINUE - F = 1.D0 - LW = 700 - CALL TN (IERROR, N, X, F, G, W, LW, SFUN) - STOP - END -C -C - SUBROUTINE SFUN (N, X, F, G) - DOUBLE PRECISION X(N), G(N), F, T -C -C ROUTINE TO EVALUATE FUNCTION (F) AND GRADIENT (G) OF THE OBJECTIVE -C FUNCTION AT THE POINT X -C - F = 0.D0 - DO 10 I = 1,N - T = X(I) - I - F = F + T*T - G(I) = 2.D0 * T -10 CONTINUE - RETURN - END //GO.SYSIN DD mains.f echo mains.out 1>&2 sed >mains.out <<'//GO.SYSIN DD mains.out' 's/^-//' - - - NIT NF CG F GTG - - - 0 1 0 3.181818171522834D+02 1.27272727D+03 - 1 2 2 6.142878372311498D+01 2.45715135D+02 - 2 3 4 1.627901350785958D-02 6.51160540D-02 - 3 4 6 1.465781507294254D-20 5.86312603D-20 - - - OPTIMAL FUNCTION VALUE = 1.465781507294254D-20 - CURRENT SOLUTION IS (AT MOST 10 COMPONENTS) - I X(I) - 1 1.000000000006170D+00 - 2 2.000000000012340D+00 - 3 3.000000000018510D+00 - 4 4.000000000024681D+00 - 5 5.000000000030851D+00 - 6 6.000000000037021D+00 - 7 7.000000000043193D+00 - 8 8.000000000049361D+00 - 9 9.000000000055532D+00 - 10 1.000000000006170D+01 //GO.SYSIN DD mains.out echo mainsb.f 1>&2 sed >mainsb.f <<'//GO.SYSIN DD mainsb.f' 's/^-//' -C*********************************************************************** -C EASY TO USE, BOUNDS ON VARIABLES -C*********************************************************************** -C MAIN PROGRAM TO MINIMIZE A FUNCTION (REPRESENTED BY THE ROUTINE SFUN) -C SUBJECT TO BOUNDS ON THE VARIABLES X -C - DOUBLE PRECISION X(50), F, G(50), W(700), LOW(50), UP(50) - INTEGER IPIVOT(50) - EXTERNAL SFUN -C -C DEFINE SUBROUTINE PARAMETERS -C N - NUMBER OF VARIABLES -C X - INITIAL ESTIMATE OF THE SOLUTION -C LOW - LOWER BOUNDS -C UP - UPPER BOUNDS -C F - ROUGH ESTIMATE OF FUNCTION VALUE AT SOLUTION -C LW - DECLARED LENGTH OF THE ARRAY W -C - N = 10 - DO 10 I = 1,N - LOW(I) = 0.D0 - UP(I) = 6.D0 - X(I) = I / FLOAT(N+1) -10 CONTINUE - F = 1.D0 - LW = 700 - CALL TNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT) - STOP - END -C -C - SUBROUTINE SFUN (N, X, F, G) - DOUBLE PRECISION X(N), G(N), F, T -C -C ROUTINE TO EVALUATE FUNCTION (F) AND GRADIENT (G) OF THE OBJECTIVE -C FUNCTION AT THE POINT X -C - F = 0.D0 - DO 10 I = 1,N - T = X(I) - I - F = F + T*T - G(I) = 2.D0 * T -10 CONTINUE - RETURN - END //GO.SYSIN DD mainsb.f echo mainsb.out 1>&2 sed >mainsb.out <<'//GO.SYSIN DD mainsb.out' 's/^-//' - - - NIT NF CG F GTG - - - 0 1 0 3.181818171522834D+02 1.27272727D+03 - 1 2 1 6.160000016784666D+01 1.82400001D+02 - 2 3 3 4.763276385235395D+01 9.05310554D+01 - 3 4 5 3.771160742225126D+01 3.48464297D+01 - 4 5 7 3.185626815775707D+01 7.42507263D+00 - 5 6 9 3.000000934245227D+01 3.73698091D-05 - 6 7 11 3.000000000007158D+01 2.86327855D-10 - 7 8 13 3.000000000000000D+01 3.52531318D-16 - - - OPTIMAL FUNCTION VALUE = 3.000000000000000D+01 - CURRENT SOLUTION IS (AT MOST 10 COMPONENTS) - I X(I) - 1 9.999999932280454D-01 - 2 2.000000005305688D+00 - 3 2.999999996494098D+00 - 4 4.000000001308861D+00 - 5 5.000000000344479D+00 - 6 6.000000000000000D+00 - 7 6.000000000000000D+00 - 8 6.000000000000000D+00 - 9 6.000000000000000D+00 - 10 6.000000000000000D+00 //GO.SYSIN DD mainsb.out