C ----------------------------------------------------------- C PROGRAM EQV C TO FIND AN OPTIMAL 4-PARAMETER SET (A1,A2,A3,A4) TO MAKE C F(X) = A1*X^(-A2)*EXP(-A3*X^A4) HAVE CLOSE APPROXIMATE C VALUES AND DERIVATIVES WITH (A1,A2,A3,A4)' AT VARIOUS X'S C (BOND LENGTHS), BUT PROJECTED TO A DIFFERENT BOUND. C C AUG. 14, 1998, JU LI (MIT) C ---------------------------- ------------------------------ PROGRAM EQV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NOP=4,MMS=20) COMMON /A/ A1,A2,A3,A4,NBOND,BOND(MMS),WB(MMS),NAME(MMS), A VALUE_OLD(MMS),DERIV_OLD(MMS),VALUE_NEW(MMS), A DERIV_NEW(MMS),DERIVATIVE_ERROR_RATIO,ITOT,MAXITER COMMON /B/ BOUND(NOP,2),GUESS(NOP) CHARACTER*16 NAME EXTERNAL FE ITOT = 1 READ (*,'("Old Parameters:")') READ *, A1,A2,A3,A4 READ (*,'("New Parameter bounds:")') READ *,BOUND(1,2),BOUND(2,2),BOUND(3,2),BOUND(4,2) READ *,BOUND(1,1),BOUND(2,1),BOUND(3,1),BOUND(4,1) READ (*,'(/,"Number of bond lengths: ",I)') NBOND READ (*,'("Name Bond Length Weight")') SUMWEIGHT = 0. DO N = 1,NBOND READ (*,'(A16,F16,F16)') NAME(N),BOND(N),WB(N) SUMWEIGHT = SUMWEIGHT + WB(N) ENDDO DO N = 1,NBOND WB(N) = WB(N) / SUMWEIGHT X = BOND(N) VALUE_OLD(N) = A1*X**(-A2)*DEXP(-A3*X**A4) DERIV_OLD(N) = VALUE_OLD(N)*(-A2/X-A3*A4*X**(A4-1.)) ENDDO READ (*,'(/,"Derivative error ratio: ",F)') A DERIVATIVE_ERROR_RATIO READ (*,'("DEL: ",F)') DEL READ (*,'("NDIV: ",I)') NDIV READ (*,'("Random number seed: ",I)') ISEED READ (*,'("Maximum number of iterations: ",I)') MAXITER IF (ISEED.EQ.0) THEN DO I=1,NOP GUESS(I) = (BOUND(I,1)+BOUND(I,2))/2. ENDDO ELSE C Initialize random number generator: KICK = MOD(ISEED,7985431) DO I = 1, KICK XX = RAN1(ISEED) ENDDO DO I=1,NOP CC = RAN1(ISEED) GUESS(I) = CC*BOUND(I,1)+(1-CC)*BOUND(I,2) ENDDO ENDIF IF (MAXITER.NE.0) CALL MINA A (FE, NOP, NDIV, DEL, BOUND, GUESS, GUESS, EMIN, IERR) CC = FE(GUESS) PRINT *,' ' PRINT *,'The average error for' WRITE (*,'(4F11.6,1x)') (GUESS(IP),IP=1,4) PRINT *,'is', cc PRINT *,' ' CALL WRAPUP STOP END SUBROUTINE WRAPUP IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NOP=4,MMS=20) COMMON /A/ A1,A2,A3,A4,NBOND,BOND(MMS),WB(MMS),NAME(MMS), A VALUE_OLD(MMS),DERIV_OLD(MMS),VALUE_NEW(MMS), A DERIV_NEW(MMS),DERIVATIVE_ERROR_RATIO,ITOT,MAXITER COMMON /B/ BOUND(NOP,2),GUESS(NOP) CHARACTER*16 NAME EXTERNAL FE PRINT *,' Length Value Deriv. F.Value F.Deriv.' DO N = 1,NBOND WRITE (*,'(F9.6,4(F11.6))') BOND(N), A VALUE_OLD(N),DERIV_OLD(N),VALUE_NEW(N),DERIV_NEW(N) ENDDO PRINT *, ' ' C Matlab script to plot out the new and old curves: print *, 'Matlab script to plot the two curves:' print *, ' ' print *, 'xmin=1.6; xmax=2.6; mesh = 1000;' print *, 'x = xmin:(xmax-xmin)/mesh:xmax;' write (*, '("alpha=[",4F11.6,"];")') A1,A2,A3,A4 write (*, '("beta=[",4F11.6,"];")') (GUESS(I),I=1,4) write (*, '("y = alpha(1) * x.^(-alpha(2)) .*", A "exp(-alpha(3)*x.^(alpha(4)));")') write (*, '("z = beta(1) * x.^(-beta(2)) .*", A "exp(-beta(3)*x.^(beta(4)));")') print *, 'clf; plot(x,y); hold on; plot(x,z,\'r\');' STOP END DOUBLE PRECISION FUNCTION FE (GUESS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NOP=4,MMS=20) COMMON /A/ A1,A2,A3,A4,NBOND,BOND(MMS),WB(MMS),NAME(MMS), A VALUE_OLD(MMS),DERIV_OLD(MMS),VALUE_NEW(MMS), A DERIV_NEW(MMS),DERIVATIVE_ERROR_RATIO,ITOT,MAXITER CHARACTER*16 NAME DOUBLE PRECISION GUESS(NOP) FE = 0. DO N = 1,NBOND X = BOND(N) VALUE_NEW(N) = GUESS(1)*X**(-GUESS(2))*DEXP(-GUESS(3)*X**GUESS(4)) DERIV_NEW(N) = VALUE_NEW(N)*( -GUESS(2)/X - A GUESS(3)*GUESS(4)*X**(GUESS(4)-1.) ) FE = FE + DERIVATIVE_ERROR_RATIO * WB(N) * A ABS((DERIV_NEW(N)-DERIV_OLD(N))/DERIV_OLD(N)) A + (1-DERIVATIVE_ERROR_RATIO) * WB(N) * A ABS((VALUE_NEW(N)-VALUE_OLD(N))/VALUE_OLD(N)) ENDDO WRITE (*,'(" step=",I6,", fe=", F6.2, 4F11.6)') A ITOT,FE,(GUESS(I),I=1,4) IF (ITOT.EQ.MAXITER) THEN PRINT *,'Maximum Iteration =',MAXITER,' reached;' PRINT *,'Optimization stop.' PRINT *,' ' CALL WRAPUP ENDIF ITOT = ITOT+1 RETURN END DOUBLE PRECISION FUNCTION RAN1(IDUM) INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV DOUBLE PRECISION AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1.D0/IM,IQ=127773,IR=2836, A NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2D-7,RNMX=1.D0-EPS) INTEGER J,K,IV(NTAB),IY SAVE IV,IY DATA IV /NTAB*0/, IY /0/ 54 IF (IDUM.LE.0.OR.IY.EQ.0) THEN IDUM=MAX(-IDUM,1) DO J=NTAB+8,1,-1 K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) IDUM=IDUM+IM IF (J.LE.NTAB) IV(J)=IDUM ENDDO IY=IV(1) ENDIF K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) IDUM=IDUM+IM J=1+IY/NDIV IY=IV(J) IV(J)=IDUM RAN1=MIN(AM*IY,RNMX) C MAKE IT (0,1) IF (RAN1.LE.0.D0) GOTO 54 IF (RAN1.GE.1.D0) GOTO 54 RETURN END SUBROUTINE MINA(FN,NV,NDIV,DEL,A,GUESS,X,FOFX,IERR) MINA.2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NOP=4) c ADDRESS.2 c SANDIA MATHEMATICAL PROGRAM LIBRARY ADDRESS.3 c APPLIED MATHEMATICS DIVISION 2613 ADDRESS. c SANDIA LABORATORIES ADDRESS.5 c ALBUQUERQUE, NEW MEXICO 87185 JUN0278.1 c CONTROL DATA 6600/7600 VERSION 7.2 MAY 1978 JUN0278.2 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.1 c ISSUED BY SANDIA LABORATORIES *MAR1378.2 c * A PRIME CONTRACTOR TO THE *MAR1378.3 c * UNITED STATES DEPARTMENT OF ENERGY *MAR1378.4 c * * * * * * * * * * * * * * * NOTICE * * * * * * * * * * * * * * * *MAR1378.5 c * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE *MAR1378.6 c * UNITED STATES GOVERNMENT. NEITHER THE UNITED STATES NOR THE *MAR1378.7 c * UNITED STATES DEPARTMENT OF ENERGY NOR ANY OF THEIR EMPLOYEES, *MAR1378.8 c * NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR EMPLOYEES *MAR1378.9 c * MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL *MAR1378.10 c * LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS OR *MAR1378.11 c * USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT OR PROCESS *MAR1378.12 c * DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE *MAR1378.13 c * OWNED RIGHTS. *MAR1378.14 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.15 c * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS *MAR1378.16 c * PART IS SAND77-1441. *MAR1378.17 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.18 c ADDRESS.9 c ORIGINAL ROUTINE WAS H2 SAND MIN, BY Z. BEISINGER AND S. BELL MINA.4 c PRESENT VERSION BY R E JONES MINA.5 c MINA.6 c ABSTRACT MINA.7 c MINA FINDS AN APPROXIMATE MINIMUM OF A REAL FUNCTION OF MINA.8 c NV VARIABLES, GIVEN AN INITIAL ESTIMATE OF THE POSITION OF MINA.9 c THE MINIMUM AND RANGES FOR EACH OF THE VARIABLES. MINA.10 c MINA USES A SELECTIVE DIRECTED SEARCH OF A SURROUNDING MINA.11 c NV-DIMENSIONAL GRID OF POINTS TO FIND A DIRECTION IN WHICH MINA.12 c THE FUNCTION DECREASES. IT THEN PROCEEDS IN THIS DIRECTION MINA.13 c AS FAR AS THE FUNCTION DECREASES, THEN DETERMINES A NEW MINA.14 c DIRECTION TO TRAVEL. WHEN NO SUCH DIRECTION IS FOUND THE MINA.15 c SEARCH INCREMENT FACTOR IS DECREASED AND THE PROCESS MINA.16 c IS REPEATED. MINA.17 c MINA.18 c DESCRIPTION OF ARGUMENTS MINA.19 c THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST MINA.20 c A(NV,2), GUESS(NV), X(NV) MINA.21 c MINA.22 c INPUT-- MINA.23 c FN - NAME OF FUNCTION OF NV VARIABLES TO BE MINIMIZED. MINA.24 c (THIS NAME MUST APPEAR IN AN EXTERNAL STATEMENT.) MINA.25 c FORM OF THE CALLING SEQUENCE MUST BE FUNCTION FN(X), MINA.26 c WHERE X IS AN ARRAY OF NV VARIABLE VALUES. THE MINA.27 c ORDERING OF THE VARIABLES IS ARBITRARY, EXCEPT MINA.28 c THAT IT MUST AGREE WITH THE ORDERING USED IN MINA.29 c ARRAYS A AND GUESS. MINA.30 c NV - NUMBER OF VARIABLES. (NV .GE. 1) MINA.31 c NDIV - NUMBER OF REFINEMENTS OF THE SEARCH INCREMENTS TO USE. MINA.32 c AT EACH REFINEMENT, THE INCREMENT IN EACH DIMENSION MINA.33 c IS DIVIDED BY 10. (USUALLY NDIV IS ABOUT 3 OR 4.) MINA.34 c DEL - FRACTION OF VARIABLE RANGE (IN EACH DIMENSION) TO USE MINA.35 c AS THE INITIAL INCREMENT (IN THAT DIMENSION) MINA.36 c A - ARRAY OF SEARCH BOUNDS, DIMENSIONED NV BY 2. MINA.37 c A(I,1) SHOULD BE THE LOWER BOUND OF THE I-TH VARIABLE. MINA.38 c A(I,2) SHOULD BE THE UPPER BOUND OF THE I-TH VARIABLE. MINA.39 c GUESS- ARRAY OF NV INITIAL VALUES. GUESS(I) SHOULD BE THE MINA.40 c INITIAL VALUE TO USE FOR THE I-TH VARIABLE. MINA.41 c MINA.42 c OUTPUT-- MINA.43 c X - ARRAY (DIMENSIONED NV) GIVING THE VALUES OF THE MINA.44 c VARIABLES AT THE MINIMUM. X(I) WILL BE THE VALUE MINA.45 c OF THE I-TH VARIABLE. MINA.46 c FOFX - FUNCTION VALUE AT THE MINIMUM MINA.47 c IERR - A STATUS CODE MINA.48 c -NORMAL CODE MINA.49 c =1 MEANS THE SEARCH FOR A MINIMUM PROCEEDED FOR THE MINA.50 c SPECIFIED NUMBER OF REFINEMENTS. MINA.51 c -ABNORMAL CODES MINA.52 c =2 MEANS NV IS GREATER THAN 50 MINA.53 c =3 MEANS A RANGE MINIMUM IS GREATER THAN THE MINA.54 c CORRESPONDING MAXIMUM MINA.55 c MINA.57 DIMENSION A(NOP,2),GUESS(NOP),X(NOP) MINA.58 DIMENSION XNOW(NOP),XNEW(NOP),R(NOP) MINA.59 c MESS.2 c CALL MLMESS(10HSANDIA7.2 ) MESS3.1 c MESS.4 c MINA.60 c INITIALIZE MINA.61 c MINA.62 IERR = 1 MINA.64 IF (NV.LE.NOP) GO TO 2 MINA.65 cc CALL ERRCHK(33,33HIN MINA , NV IS GREATER THAN 100.) MINA.66 write(6,*)'nv is greater than 100' IERR = 2 MINA.67 RETURN MINA.68 2 NX = NV MINA.69 IDIV = 0 MINA.70 DO 5 I=1,NX MINA.71 XNOW(I) = GUESS(I) MINA.72 IF (XNOW(I).LT.A(I,1)) XNOW(I) = A(I,1) MINA.73 IF (XNOW(I).GT.A(I,2)) XNOW(I) = A(I,2) MINA.74 IF (A(I,1)-A(I,2)) 5,5,4 MINA.75 4 write(6,*) 'range minimum greater than maximum' print *,A(i,1),A(i,2) c 4 CALL ERRCHK(46,46HIN MINA , RANGE MINIMUM GREATER THAN MAXIMUM.) MINA.76 IERR = 3 MINA.77 RETURN MINA.78 5 R(I) = A(I,2)-A(I,1) MINA.79 DELTA = DEL MINA.80 IF (DELTA.LE.0.0) DELTA = 0.1 MINA.81 FNOW = FN(XNOW) MINA.82 c MINA.83 c FIND NEW DIRECTION MINA.84 c MINA.85 7 DO 8 I=1,NX MINA.86 8 XNEW(I) = XNOW(I) MINA.87 FOLD = FNOW MINA.88 10 DO 40 I=1,NX MINA.89 IF (XNOW(I).GE.A(I,2)) GO TO 20 MINA.90 XNEW(I) = MIN(XNOW(I)+DELTA*R(I),A(I,2)) FNEW = FN(XNEW) MINA.92 IF (FNEW.LT.FNOW) GO TO 30 MINA.93 20 IF (XNOW(I) .LE. A(I,1)) GO TO 25 MINA.94 XNEW(I) = MAX(XNOW(I)-DELTA*R(I),A(I,1)) FNEW = FN(XNEW) MINA.96 IF (FNEW.LT.FNOW) GO TO 30 MINA.97 25 XNEW(I) = XNOW(I) MINA.98 GO TO 40 MINA.99 30 FNOW = FNEW MINA.100 40 CONTINUE MINA.101 ISTEP = 1 MINA.102 c MINA.103 c REFINE IF NEEDED MINA.104 c MINA.105 IF (FNOW.LT.FOLD) GO TO 50 MINA.106 IF (IDIV.GE.NDIV) GO TO 100 MINA.107 DELTA = DELTA*0.1 MINA.108 IDIV = IDIV+1 MINA.109 GO TO 10 MINA.110 c MINA.111 c TRY TO CONTINUE IN CHOSEN DIRECTION MINA.112 c MINA.113 50 ICHNG = 0 MINA.114 FAc = 1.0 MINA.115 IF ((ISTEP/10)*10.EQ.ISTEP) FAc = 2.0 MINA.116 DO 60 I=1,NX MINA.117 DX = (XNEW(I)-XNOW(I))*FAc MINA.118 XNOW(I) = XNEW(I) MINA.119 IF (DX) 52,54,56 MINA.120 52 XNEW(I) = MAX(XNOW(I)+DX,A(I,1)) IF (XNEW(I).LT.XNOW(I)) ICHNG = 1 MINA.122 GO TO 60 MINA.123 54 XNEW(I) = XNOW(I) MINA.124 GO TO 60 MINA.125 56 XNEW(I) = MIN(XNOW(I)+DX,A(I,2)) IF (XNEW(I).GT.XNOW(I)) ICHNG = 1 MINA.127 60 CONTINUE MINA.128 IF (ICHNG.EQ.0) GO TO 7 MINA.129 FNEW = FN(XNEW) MINA.130 IF (FNEW.GE.FNOW) GO TO 7 MINA.131 FNOW = FNEW MINA.132 ISTEP = ISTEP+1 MINA.133 GO TO 50 MINA.134 c MINA.135 c RETURN ANSWERS MINA.136 c MINA.137 100 FOFX = FOLD MINA.138 DO 110 I=1,NX MINA.139 110 X(I) = XNOW(I) MINA.140 RETURN MINA.141 END MINA.142 SUBROUTINE SIMIN (F,K,EPS,ANS,S,NEV,ICONT,Y) SIMIN.2 IMPLICIT double precision (A-H,O-Z) c ADDRESS.2 c SANDIA MATHEMATICAL PROGRAM LIBRARY ADDRESS.3 c APPLIED MATHEMATICS DIVISION 2613 ADDRESS.4 c SANDIA LABORATORIES ADDRESS.5 c ALBUQUERQUE, NEW MEXICO 87185 JUN0278.1 c CONTROL DATA 6600/7600 VERSION 7.2 MAY 1978 JUN0278.2 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.1 c ISSUED BY SANDIA LABORATORIES *MAR1378.2 c * A PRIME CONTRACTOR TO THE *MAR1378.3 c * UNITED STATES DEPARTMENT OF ENERGY *MAR1378.4 c * * * * * * * * * * * * * * * NOTICE * * * * * * * * * * * * * * * *MAR1378.5 c * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE *MAR1378.6 c * UNITED STATES GOVERNMENT. NEITHER THE UNITED STATES NOR THE *MAR1378.7 c * UNITED STATES DEPARTMENT OF ENERGY NOR ANY OF THEIR EMPLOYEES, *MAR1378.8 c * NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR EMPLOYEES *MAR1378.9 c * MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL *MAR1378.10 c * LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS OR *MAR1378.11 c * USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT OR PROCESS *MAR1378.12 c * DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE *MAR1378.13 c * OWNED RIGHTS. *MAR1378.14 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.15 c * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS *MAR1378.16 c * PART IS SAND77-1441. *MAR1378.17 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *MAR1378.18 c ADDRESS.9 c ORIGINAL ROUTINE BY L F SHAMPINE, AS DESCRIBED IN REF.1 BELOW. SIMIN.4 c PREPARATION FOR MATH LIBRARY BY R E JONES. SIMIN.5 c SIMIN.6 c ABSTRACT SIMIN.7 c SIMIN FINDS AN APPROXIMATE MINIMUM OF A REAL FUNCTION OF K SIMIN.8 c VARIABLES, GIVEN AN INITIAL ESTIMATE OF THE POSITION OF THE SIMIN.9 c MINIMUM. THE SIMPLEX METHOD IS USED. SEE REFERENCE 1 BELOW SIMIN.10 c FOR A FULL EXPLANATION OF THIS METHOD. BRIEFLY, A SET OF SIMIN.11 c K+1 POINTS IN K-DIMENSIONAL SPACE IS CALLED A SIMPLEX. SIMIN.12 c THE MINIMIZATION PROCESS ITERATES BY REPLACING THE POINT SIMIN.13 c WITH THE LARGEST FUNCTION VALUE BY A NEW POINT WITH A SIMIN.14 c SMALLER FUNCTION VALUE. ITERATION CONTINUES UNTIL ALL THE SIMIN.15 c POINTS CLUSTER SUFFICIENTLY CLOSE TO A MINIMUM. SIMIN.16 c SIMIN.17 c REFERENCES SIMIN.18 c 1. L F SHAMPINE, A ROUTINE FOR UNCONSTRAINED OPTIMIZATION, SIMIN.19 c SC-TM-72130 OR SC-RR-720657 SIMIN.20 c 2. J A NELDER AND R MEAD, A SIMPLEX METHOD FOR FUNCTION SIMIN.21 c MINIMIZATION, COMPUTER JOURNAL, 7(1965) 308-313 SIMIN.22 c SIMIN.23 c DESCRIPTION OF PARAMETERS SIMIN.24 c --INPUT-- SIMIN.25 c F - NAME OF FUNCTION OF K VARIABLES TO BE MINIMIZED. SIMIN.26 c (THIS NAME MUST APPEAR IN AN EXTERNAL STATEMENT.) SIMIN.27 c FORM OF THE CALLING SEQUENCE MUST BE FUNCTION F(X), SIMIN.28 c WHERE X IS AN ARRAY OF K VARIABLES. SIMIN.29 c K - THE NUMBER OF VARIABLES. K MUST BE AT LEAST 2. SIMIN.30 c NORMALLY K SHOULD BE LESS THAN ABOUT 10, AS SIMIN SIMIN.31 c BECOMES LESS EFFECTIVE FOR LARGER VALUES OF K. SIMIN.32 c EPS- THE CONVERGENCE CRITERION. LET YAVG BE THE AVERAGE SIMIN.33 c VALUE OF THE FUNCTION F AT THE K+1 POINTS OF THE SIMIN.34 c SIMPLEX, AND LET R BE THEIR STANDARD ERROR. (THAT IS, SIMIN.35 c THE ROOT-MEAN-SQUARE OF THE SET OF VALUES (Y(I)-YAVG), SIMIN.36 c WHERE Y(I) IS THE FUNCTION VALUE AT THE I-TH POINT OF SIMIN.37 c THE SIMPLEX.) THEN-- SIMIN.38 c IF EPS.GT.0, CONVERGENCE IS OBTAINED IF R.LE.EPS. SIMIN.39 c IF EPS.LT.0, CONVERGENCE IS IF R.LE.ABS(EPS*YAVG). SIMIN.40 c IF EPS=0, THE PROCESS WILL NOT CONVERGE BUT INSTEAD WILL SIMIN.41 c QUIT WHEN NEV FUNCTION EVALUATIONS HAVE BEEN USED. SIMIN.42 c ANS- AN ARRAY OF LENGTH K CONTAINING A GUESS FOR THE LOCATION SIMIN.43 c OF A MINIMUM OF F. SIMIN.44 c S - A SCALE PARAMETER, WHICH MAY BE A SIMPLE VARIABLE OR AN SIMIN.45 c ARRAY OF LENGTH K. USE OF AN ARRAY IS SIGNALLED BY SIMIN.46 c SETTING S(1) NEGATIVE. SIMIN.47 c -SIMPLE VARIABLE CASE. HERE S IS THE LENGTH OF EACH SIMIN.48 c SIDE OF THE INITIAL SIMPLEX. THUS, THE INITIAL SEARCH SIMIN.49 c RANGE IS THE SAME FOR ALL THE VARIABLES. SIMIN.50 c -ARRAY CASE. HERE THE LENGTH OF SIDE I OF THE INITIAL SIMIN.51 c SIMPLEX IS ABS(S(I)). THUS, THE INITIAL SEARCH RANGE SIMIN.52 c MAY BE DIFFERENT FOR DIFFERENT VARIABLES. SIMIN.53 c NOTE-- THE VALUE(S) USED FOR S ARE NOT VERY CRITICAL. SIMIN.54 c ANY REASONABLE GUESS SHOULD DO O.K. SIMIN.55 c NEV- THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS TO BE USED. SIMIN.56 c (THE ACTUAL NUMBER USED MAY EXCEED THIS SLIGHTLY SO THE SIMIN.57 c LAST SEARCH ITERATION MAY BE COMPLETED.) SIMIN.58 c ICONT - ICONT SHOULD BE ZERO ON ANY CALL TO SIMIN WHICH SIMIN.59 c IS NOT A CONTINUATION OF A PREVIOUS CALL. SIMIN.60 c IF ICONT=1 THE PROBLEM WILL BE CONTINUED. IN THIS SIMIN.61 c CASE THE WORK ARRAY Y MUST BE THE SAME ARRAY THAT WAS SIMIN.62 c USED IN THE CALL THAT IS BEING CONTINUED (AND THE VALUES SIMIN.63 c IN IT MUST BE UNCHANGED). THE REASON FOR THIS IS THAT SIMIN.64 c IF ICONT=1 THEN THE ARGUMENT S IS IGNORED AND THE SIMPLEX SIMIN.65 c AND RELATED FUNCTION VALUES THAT WERE STORED IN ARRAY Y SIMIN.66 c DURING A PREVIOUS EXECUTION ARE USED TO CONTINUE THAT SIMIN.67 c PREVIOUS PROBLEM. SIMIN.68 c Y - A WORK ARRAY CONTAINING AT LEAST K*K + 5*K + 1 WORDS. SIMIN.69 c IF ICONT=1 THIS MUST BE THE SAME ARRAY USED IN THE CALL SIMIN.70 c THAT IS BEING CONTINUED. SIMIN.71 c --OUTPUT-- SIMIN.72 c ANS- ANS WILL CONTAIN THE LOCATION OF THE POINT WITH THE SIMIN.73 c SMALLEST VALUE OF THE FUNCTION THAT WAS FOUND. SIMIN.74 c S - IN THE SIMPLE VARIABLE CASE S WILL BE RETURNED AS THE SIMIN.75 c AVERAGE DISTANCE FROM THE VERTICES TO THE CENTROID OF SIMIN.76 c THE SIMPLEX. SIMIN.77 c IN THE ARRAY CASE S(I) WILL BE RETURNED AS THE AVERAGE SIMIN.78 c DISTANCE IN THE I-TH DIMENSION OF VERTICES FROM SIMIN.79 c THE CENTROID. (S(1) WILL BE NEGATED.) SIMIN.80 c NOTE-- THE VALUE(S) RETURNED IN S ARE USEFUL FOR SIMIN.81 c ASSESSING THE FLATNESS OF THE FUNCTION NEAR THE SIMIN.82 c MINIMUM. THE LARGER THE VALUE OF S (FOR A GIVEN SIMIN.83 c VALUE OF EPS), THE FLATTER THE FUNCTION. SIMIN.84 c NEV- NEV WILL BE THE COUNT OF THE ACTUAL NUMBER OF FUNCTION SIMIN.85 c EVALUATIONS USED. SIMIN.86 c Y - WILL CONTAIN ALL DATA NEEDED TO CONTINUE THE MINIMIZATION SIMIN.87 c SEARCH EFFICIENTLY IN A SUBSEQUENT CALL. SIMIN.88 c NOTE -- THE FIRST K+1 ELEMENTS OF Y WILL CONTAIN THE SIMIN.89 c FUNCTION VALUES AT THE K+1 POINTS OF THE LATEST SIMPLEX. SIMIN.90 c THE NEXT K*(K+1) ELEMENTS OF Y WILL BE THE K+1 POINTS SIMIN.91 c OF THE SIMPLEX (IN EXACT CORRESPONDENSE TO THE ARRAY SIMIN.92 c P DISCUSSED IN REFERENCE 1 ABOVE). THE REMAINING 3*K SIMIN.93 c WORDS ARE TEMPORARY WORKING STORAGE ONLY. SIMIN.94 c SIMIN.96 DIMENSION ANS(K),S(K),Y(1) SIMIN.97 EXTERNAL F SIMIN.98 c MESS.2 c CALL MLMESS(10HSANDIA7.2 ) MESS3.1 c MESS.4 c SIMIN.99 IF(K.GE.2 .AND. S(1).NE.0.) GO TO 10 SIMIN.101 write(6,*) 's(1)=0 or k is less than 2' c CALL ERRCHK(39,39HIN SIMIN , S(1)=0. OR K IS LESS THAN 2.) SIMIN.102 RETURN SIMIN.103 c 10 IF (K.GT.100) CALL ERRCHK(31,31HIN SIMIN , K IS LARGER THAN 100.) SIMIN.104 10 write(6,*) 'k is large than 100' c SIMIN.105 IP = K+2 SIMIN.106 Ic = IP+K*(K+1) SIMIN.107 IR = IC+K SIMIN.108 IRR = IR+K SIMIN.109 CALL SIMINA(F,K,EPS,ANS,S,NEV,ICONT,Y,Y(IP),Y(IC),Y(IR),Y(IRR)) SIMIN.110 RETURN SIMIN.111 END SIMIN.112 SUBROUTINE SIMINA(F,K,EPS,ANS,S,NEV,ICONT,Y,P,PC,PR,PRR) SIMIN.113 IMPLICIT double precision (A-H,O-Z) DIMENSION ANS(K),S(K),Y(3),P(K,3),PC(K),PR(K),PRR(K) SIMIN.114 DATA ALPHA,BETA,GAMMA/1.0,0.5,2.0/ SIMIN.115 c SIMIN.116 c SIMINA IS A CORE MINIMIZATION ROUTINE CALLED ONLY BY SIMIN. SIMIN.117 c SIMIN.118 c INITIALIZATION SIMIN.119 c SIMIN.120 KOUNT = 0 SIMIN.122 KK = K SIMIN.123 IF (KK.LE.1) GO TO 99 SIMIN.124 ONEK = 1.0/FLOAT(KK) SIMIN.125 KP1 = KK+1 SIMIN.126 ONEKP1 = 1.0/FLOAT(KP1) SIMIN.127 TOL = FLOAT(KP1)*EPS**2 SIMIN.128 c SIMIN.129 c INITIAL SIMPLEX SIMIN.130 c SIMIN.131 IF (ICONT.GE.1) GO TO 10 SIMIN.132 IF (S(1)) 4,99,1 SIMIN.133 1 SKP1 = S(1)*ONEKP1 SIMIN.134 DO 2 I=1,KP1 SIMIN.135 DO 2 J=1,KK SIMIN.136 2 P(J,I) = ANS(J) - SKP1 SIMIN.137 DO 3 J=1,KK SIMIN.138 3 P(J,J+1) = P(J,J+1) + S(1) SIMIN.139 GO TO 7 SIMIN.140 c SIMIN.141 4 DO 5 I=1,KP1 SIMIN.142 DO 5 J=1,KK SIMIN.143 5 P(J,I) = ANS(J) - ABS(S(J))*ONEKP1 SIMIN.144 DO 6 J=1,KK SIMIN.145 6 P(J,J+1) = P(J,J+1) + ABS(S(J)) SIMIN.146 c SIMIN.147 c FUNCTION VALUES FOR INITIAL SIMPLEX SIMIN.148 c SIMIN.149 7 I1 = 1 SIMIN.150 DO 8 I=1,KP1 SIMIN.151 Y(I) = F(P(1,I)) SIMIN.152 IF (Y(I).GT.Y(I1)) I1 = I SIMIN.153 8 CONTINUE SIMIN.154 YANS = F(ANS) SIMIN.155 KOUNT = KP1+1 SIMIN.156 IF (YANS.GE.Y(I1)) GO TO 10 SIMIN.157 Y(I1) = YANS SIMIN.158 DO 9 J=1,KK SIMIN.159 9 P(J,I1) = ANS(J) SIMIN.160 c SIMIN.161 c RE-START / NEXT ITERATION SIMIN.162 c IF K.LT.0 VALUES IN THE P AND Y ARRAYS (AND ONLY THESE VALUES) SIMIN.163 c WILL NOT HAVE BEEN DEFINED IN THIS CALL. THIS IS NON-ANSI USAGE. SIMIN.164 c SIMIN.165 c FIRST FIND LARGEST, SECOND LARGEST, AND SMALLEST FUNCTION VALUES. SIMIN.166 c SIMIN.167 10 I1 = 1 SIMIN.168 IL = 1 SIMIN.169 DO 12 I=2,KP1 SIMIN.170 IF (Y(I).LT.Y(IL)) IL = I SIMIN.171 IF (Y(I).GT.Y(I1)) I1 = I SIMIN.172 12 CONTINUE SIMIN.173 I2 = IL SIMIN.174 DO 13 I=1,KP1 SIMIN.175 IF (I.EQ.I1) GO TO 13 SIMIN.176 IF (Y(I).GT.Y(I2)) I2 = I SIMIN.177 13 CONTINUE SIMIN.178 c SIMIN.179 c COMPUTE CENTROID, LEAVING OUT P(*,I1) SIMIN.180 c SIMIN.181 DO 15 J=1,KK SIMIN.182 SUM = 0.0 SIMIN.183 DO 14 I=1,KP1 SIMIN.184 IF (I.EQ.I1) GO TO 14 SIMIN.185 SUM = SUM + P(J,I) SIMIN.186 14 CONTINUE SIMIN.187 15 PC(J) = SUM*ONEK SIMIN.188 c SIMIN.189 c FORM REFLECTED POINT AND TEST SIMIN.190 c SIMIN.191 DO 20 J=1,KK SIMIN.192 20 PR(J) = PC(J) + ALPHA*(PC(J)-P(J,I1)) SIMIN.193 YR = F(PR) SIMIN.194 KOUNT = KOUNT+1 SIMIN.195 IF (YR.LT.Y(IL)) GO TO 30 SIMIN.196 IF (YR.GE.Y(I2)) GO TO 40 SIMIN.197 c SIMIN.198 c ACCEPT REFLECTED POINT SIMIN.199 c SIMIN.200 21 Y(I1) = YR SIMIN.201 DO 22 J=1,KK SIMIN.202 22 P(J,I1) = PR(J) SIMIN.203 GO TO 60 SIMIN.204 c SIMIN.205 c EXPAND IN FAVORABLE DIRECTION AND TEST SIMIN.206 c SIMIN.207 30 DO 31 J=1,KK SIMIN.208 31 PRR(J) = PR(J) + GAMMA*(PR(J)-PC(J)) SIMIN.209 YRR = F(PRR) SIMIN.210 KOUNT = KOUNT+1 SIMIN.211 IF (YRR.GE.YR) GO TO 21 SIMIN.212 c SIMIN.213 c ACCEPT EXPANDED POINT SIMIN.214 c SIMIN.215 Y(I1) = YRR SIMIN.216 DO 32 J=1,KK SIMIN.217 32 P(J,I1) = PRR(J) SIMIN.218 GO TO 60 SIMIN.219 c SIMIN.220 c DECIDE WHETHER TO ACCEPT REFLECTED POINT. SIMIN.221 c SIMIN.222 40 IF (YR.GE.Y(I1)) GO TO 42 SIMIN.223 Y(I1) = YR SIMIN.224 DO 41 J=1,KK SIMIN.225 41 P(J,I1) = PR(J) SIMIN.226 c SIMIN.227 c TRY CONTRACTION. SIMIN.228 c SIMIN.229 42 DO 43 J=1,KK SIMIN.230 43 PR(J) = PC(J) + BETA*(P(J,I1)-PC(J)) SIMIN.231 YCT = F(PR) SIMIN.232 KOUNT = KOUNT+1 SIMIN.233 IF (YCT.GT.Y(I1)) GO TO 50 SIMIN.234 Y(I1) = YCT SIMIN.235 DO 44 J=1,KK SIMIN.236 44 P(J,I1) = PR(J) SIMIN.237 GO TO 60 SIMIN.238 c SIMIN.239 c ALL EFFORTS FAILED. SHRINK THE SIMPLEX ABOUT BEST POINT. SIMIN.240 c SIMIN.241 50 DO 52 I=1,KP1 SIMIN.242 IF (I.EQ.IL) GO TO 52 SIMIN.243 DO 51 J=1,KK SIMIN.244 51 P(J,I) = 0.5*(P(J,I)+P(J,IL)) SIMIN.245 Y(I) = F(P(1,I)) SIMIN.246 52 CONTINUE SIMIN.247 KOUNT = KOUNT+KP1 SIMIN.248 c SIMIN.249 c CHECK FOR CONVERGENCE SIMIN.250 c SIMIN.251 60 IF (KOUNT.GE.NEV) GO TO 65 SIMIN.252 IF (EPS.EQ.0.0) GO TO 10 SIMIN.253 SUM = 0.0 SIMIN.254 DO 61 I=1,KP1 SIMIN.255 61 SUM = SUM + Y(I) SIMIN.256 YAVG = SUM*ONEKP1 SIMIN.257 SUM = 0.0 SIMIN.258 DO 62 I=1,KP1 SIMIN.259 62 SUM = SUM + (Y(I)-YAVG)**2 SIMIN.260 IF (EPS) 64,63,63 SIMIN.261 63 IF (SUM-TOL) 65,65,10 SIMIN.262 64 IF (SUM-TOL*ABS(YAVG)) 65,65,10 SIMIN.263 c SIMIN.264 c CONVERGENCE OBTAINED. SIMIN.265 c SIMIN.266 c COMPUTE CENTROID SIMIN.267 c SIMIN.268 65 DO 68 J=1,KK SIMIN.269 SUM = 0.0 SIMIN.270 DO 67 I=1,KP1 SIMIN.271 67 SUM = SUM+P(J,I) SIMIN.272 68 PC(J) = SUM*ONEKP1 SIMIN.273 IF (S(1)) 73,69,69 SIMIN.274 c SIMIN.275 c COMPUTE S(1) AS AVERAGE DISTANCE OF VERTICES FROM CENTROID. SIMIN.276 c SIMIN.277 69 DIST = 0.0 SIMIN.278 DO 71 I=1,KP1 SIMIN.279 SUM = 0.0 SIMIN.280 DO 70 J=1,KK SIMIN.281 70 SUM = SUM + (P(J,I)-PC(J))**2 SIMIN.282 71 DIST = DIST + SQRT(SUM) SIMIN.283 S(1) = DIST*ONEKP1 SIMIN.284 GO TO 80 SIMIN.285 c SIMIN.286 c COMPUTE S(J) AS AVERAGE DISTANCE IN J-TH DIMENSION OF SIMIN.287 c VERTICES FROM THE CENTROID. SIMIN.288 c SIMIN.289 73 DO 75 J=1,KK SIMIN.290 SUM = 0.0 SIMIN.291 DO 74 I=1,KP1 SIMIN.292 74 SUM = SUM + ABS(P(J,I)-PC(J)) SIMIN.293 75 S(J) = SUM*ONEKP1 SIMIN.294 S(1) = -S(1) SIMIN.295 c SIMIN.296 c RETURN P(*,IL) AS ANSWER SIMIN.297 c SIMIN.298 80 IL = 1 SIMIN.299 DO 82 I=2,KP1 SIMIN.300 IF (Y(I).LT.Y(IL)) IL = I SIMIN.301 82 CONTINUE SIMIN.302 DO 84 J=1,KK SIMIN.303 84 ANS(J) = P(J,IL) SIMIN.304 NEV = KOUNT SIMIN.305 RETURN SIMIN.306 c SIMIN.307 c ERROR MESSAGE SIMIN.308 c SIMIN.309 c 99 CALL ERRCHK(39,39HIN SIMINA, S(1)=0. OR K IS LESS THAN 2.) SIMIN.310 99 write(6,*) 's(1)=0. or k is less than 2.' RETURN SIMIN.311 END SIMIN.312