#! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/PackingList # Fortran90/ # Fortran90/Drivers/ # Fortran90/Drivers/Sp/ # Fortran90/Drivers/Sp/data1 # Fortran90/Drivers/Sp/data1.7 # Fortran90/Drivers/Sp/driver1.f90 # Fortran90/Drivers/Sp/driver2.f90 # Fortran90/Drivers/Sp/res.8 # Fortran90/Drivers/Sp/res1 # Fortran90/Drivers/Sp/res2 # Fortran90/Src/ # Fortran90/Src/Sp/ # Fortran90/Src/Sp/src.f90 # This archive created: Mon Oct 18 15:54:03 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'PackingList' then echo shar: will not over-write existing file "'PackingList'" else cat << SHAR_EOF > 'PackingList' Algorithm:734 Language:Fortran90 Precision:Sp Src:src.f90 Driver_0:@Src driver1.f90 DriverLib_0:port Data_0:stdin=data1 data1.7 Res_0:stdout=res1 res1.8 Driver_1:@Src driver2.f90 DriverLib_1:port Res_1:stdout=res2 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran90' then mkdir 'Fortran90' fi cd 'Fortran90' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << SHAR_EOF > 'data1' -3 SHAR_EOF fi # end of overwriting check if test -f 'data1.7' then echo shar: will not over-write existing file "'data1.7'" else cat << SHAR_EOF > 'data1.7' &Args state = 0, dotrace = false / &Args state = 2, dotrace = false / &Args state = 0, dotrace = false, derv = 4 / &Args state = 0, dotrace = false, derv = 2 / &Args state = 0, dotrace = false, meth = 3, mterms = 2 / &Args state = 2, dotrace = false, meth = 7 / &Args state = 0, dotrace = false, meth = 3, mterms = 2, update = 1 / &Args state = 0, dotrace = false, run8 = true / &Args state = 0, dotrace = false, meth = 6 / &Args state = 0, dotrace = false, expens = 300, meth = 6 / SHAR_EOF fi # end of overwriting check if test -f 'driver1.f90' then echo shar: will not over-write existing file "'driver1.f90'" else cat << SHAR_EOF > 'driver1.f90' C--**--CH2468--734--P:RW--22:9:1999 ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> test_functions.f90 MODULE test_functions USE precisions, ONLY : stnd, short IMPLICIT NONE PUBLIC !------ -------------! INTEGER, PARAMETER :: nargs = 10 REAL (stnd) :: arguments(nargs) INTEGER :: funcno !--------- DATA FOR TESTPACK FUNCTION ARGAUS REAL (stnd), DIMENSION (15), PARAMETER :: argasy = (/ 9.000D-4, & 4.400D-3, 1.750D-2, 5.400D-2, 1.295D-1, 2.420D-1, 3.521D-1, 3.989D-1, & 3.521D-1, 2.420D-1, 1.295D-1, 5.400D-2, 1.750D-2, 4.400D-3, 9.000D-4/) !--------- DATA FOR TESTPACK FUNCTION BARD70 REAL (stnd), DIMENSION (15), PARAMETER :: bard7y = (/ .14D0, .18D0, & .22D0, .25D0, .29D0, .32D0, .35D0, .39D0, .37D0, .58D0, .73D0, .96D0, & 1.34D0, 2.10D0, 4.39D0/) !--------- DATA FOR TESTPACK FUNCTION CHNRSN REAL (stnd), DIMENSION (50), PARAMETER :: chnrsn = (/ 1.25D0, 1.40D0, & 2.40D0, 1.40D0, 1.75D0, 1.20D0, 2.25D0, 1.20D0, 1.00D0, 1.10D0, & 1.50D0, 1.60D0, 1.25D0, 1.25D0, 1.20D0, 1.20D0, 1.40D0, 0.50D0, & 0.50D0, 1.25D0, 1.80D0, 0.75D0, 1.25D0, 1.40D0, 1.60D0, 2.00D0, & 1.00D0, 1.60D0, 1.25D0, 2.75D0, 1.25D0, 1.25D0, 1.25D0, 3.00D0, & 1.50D0, 2.00D0, 1.25D0, 1.40D0, 1.80D0, 1.50D0, 2.20D0, 1.40D0, & 1.50D0, 1.25D0, 2.00D0, 1.50D0, 1.25D0, 1.40D0, 0.60D0, 1.50D0/) !--------- DATA FOR TESTPACK FUNCTION HIMM32 REAL (stnd), DIMENSION (7), PARAMETER :: him32a = (/ 0.0D0, 4.28D-4, & 1.0D-3, 1.61D-3, 2.09D-3, 3.48D-3, 5.25D-3/) REAL (stnd), DIMENSION (7), PARAMETER :: him32b = (/ 7.391D0, 1.118D1, & 1.644D1, 1.62D1, 2.22D1, 2.402D1, 3.132D1/) !--------- DATA FOR TESTPACK FUNCTION KOWOSB REAL (stnd), DIMENSION (11), PARAMETER :: kowosu = (/ 4.0D0, 2.0D0, & 1.0D0, 0.5D0, 0.25D0, 0.167D0, 0.125D0, 0.1D0, 0.0833D0, 0.0714D0, & 0.0625D0/) REAL (stnd), DIMENSION (11), PARAMETER :: kowosy = (/ 0.1957D0, & 0.1947D0, 0.1735D0, 0.1600D0, 0.0844D0, 0.0627D0, 0.0456D0, 0.0342D0, & 0.0323D0, 0.0235D0, 0.0246D0/) !--------- DATA FOR TESTPACK FUNCTION MEYER REAL (stnd), DIMENSION (16), PARAMETER :: mey = (/ 3.478D4, 2.861D4, & 2.365D4, 1.963D4, 1.637D4, 1.372D4, 1.154D4, 9.744D3, 8.261D3, & 7.030D3, 6.005D3, 5.147D3, 4.427D3, 3.820D3, 3.307D3, 2.872D3/) !--------- DATA FOR TESTPACK FUNCTION ORTOIT REAL (stnd), DIMENSION (33), PARAMETER :: orbeta = (/ 1.0D0, 1.5D0, & 1.0D0, 0.1D0, 1.5D0, 2.0D0, 1.0D0, 1.5D0, 3.0D0, 2.0D0, 1.0D0, 3.0D0, & 0.1D0, 1.5D0, 0.15D0, 2.0D0, 1.0D0, 0.1D0, 3.0D0, 0.1D0, 1.2D0, 1.0D0, & 0.1D0, 2.0D0, 1.2D0, 3.0D0, 1.5D0, 3.0D0, 2.0D0, 1.0D0, 1.2D0, 2.0D0, & 1.0D0/) REAL (stnd), DIMENSION (33), PARAMETER :: od = (/ 5.0D0, 5.0D0, 5.0D0, & 2.5D0, 6.0D0, 6.0D0, 5.0D0, 6.0D0, 10.0D0, 6.0D0, 5.0D0, 9.0D0, 2.0D0, & 7.0D0, 2.5D0, 6.0D0, 5.0D0, 2.0D0, 9.0D0, 2.0D0, 5.0D0, 5.0D0, 2.5D0, & 5.0D0, 6.0D0, 10.0D0, 7.0D0, 10.0D0, 6.0D0, 5.0D0, 4.0D0, 4.0D0, & 4.0D0/) INTEGER, DIMENSION (50), PARAMETER :: a = (/ -31, -1, -2, -4, -6, -8, & -10, -12, + 11, + 13, -14, -16, + 9, -18, + 5, + 20, -21, -19, -23, & + 7, -25, -28, -29, -32, + 3, -33, -35, -36, + 30, -37, + 38, -39, & -40, -41, -44, -46, + 42, + 45, + 48, -50, + 26, + 34, -43, + 15, & + 17, + 24, -47, -49, -22, -27 /) INTEGER, DIMENSION (56), PARAMETER :: b = (/ -1, + 2, -3, + 4, -5, + 6, & -7, + 8, -9, + 10, -11, + 12, -13, + 14, -15, + 16, -17, + 18, -19, & -20, 0, + 22, + 23, -24, + 25, -26, + 27, -28, + 29, -30, + 31, -32, & + 33, -34, -35, + 21, -36, + 37, -38, -39, -40, + 41, -42, + 43, + 44, & -50, + 45, + 46, -47, -48, -49, 0, 0, 0, 0, 0 /) !--------- DATA FOR TESTPACK FUNCTION OSBRN1 REAL (stnd), DIMENSION (33), PARAMETER :: osb1y = (/ .844D0, .908D0, & .932D0, .936D0, .925D0, .908D0, .881D0, .850D0, .818D0, .784D0, & .751D0, .718D0, .685D0, .658D0, .628D0, .603D0, .580D0, .558D0, & .538D0, .522D0, .506D0, .490D0, .478D0, .467D0, .457D0, .448D0, & .438D0, .431D0, .424D0, .420D0, .414D0, .411D0, .406D0/) !--------- DATA FOR TESTPACK FUNCTION OSBRN2 REAL (stnd), DIMENSION (65), PARAMETER :: osb2y = (/ 1.366D0, 1.191D0, & 1.112D0, 1.013D0, .991D0, .885D0, .831D0, .847D0, .786D0, .725D0, & .746D0, .679D0, .608D0, .655D0, .616D0, .606D0, .602D0, .626D0, & .651D0, .724D0, .649D0, .649D0, .694D0, .644D0, .624D0, .661D0, & .612D0, .558D0, .533D0, .495D0, .50D0, .423D0, .395D0, .375D0, .372D0, & .391D0, .396D0, .405D0, .428D0, .429D0, .523D0, .562D0, .607D0, & .653D0, .672D0, .708D0, .633D0, .668D0, .645D0, .632D0, .591D0, & .559D0, .597D0, .625D0, .739D0, .710D0, .729D0, .720D0, .636D0, & .581D0, .428D0, .292D0, .162D0, .098D0, .054D0/) CONTAINS SUBROUTINE functions(x,f,g,ifg) USE reals, ONLY : c100, c20, c200, c36, c40, c400, c600, eight, fifth, & five, four, one, six, ten, tenth, three, two, zero USE supp_codes, ONLY : justf, justg, noforg, ok USE num_constants, ONLY : machhuge IMPLICIT NONE !-------------! ! ARGUMENTS: INTEGER (short), INTENT (INOUT) :: ifg REAL (stnd), INTENT (IN) :: x(:) REAL (stnd) :: work(2*SIZE(x)) REAL (stnd), INTENT (OUT) :: f, g(SIZE(x)) ! STATUS: ! SYSTEM DEPENDENCE: NONE ! DESCRIPTION: ! THIS TEST FUNCTION EVALUATES ONE OF THE STANDARD TEST ! FUNCTIONS PROVIDED WITH TESTPACK. THE ARGUMENTS IN THE CALLING ! SEQUENCE HAVE PRECISELY THE SAME MEANING AS IN THE ROUTINE EVALUATE_F. ! THE TEST FUNCTION TO USE IS SELECTED THRU THE MODULE TEST_FUNCTIONS. ! THE VALUE OF THE INTEGER, FUNCNO, SPECIFIES WHICH OF THE TEST FUNCTIONS ! IS TO BE USED; THE FUNCTION IS CHOSEN USING A CASE CONSTRUCT. ! SOME OF THE FUNCTIONS NEED SPECIAL ARGUMENTS (OTHER THAN THE ! VALUE OF X); THESE ARE PROVIDED THROUGH THE MODULE TEST_FUNCTIONS. A ! MAXIMUM OF TEN ARGUMENTS ARE PROVIDED. IF THE MAXIMUM NUMBER OF ! ARGUMENTS IS TO BE INCREASED, THE PARAMETER NARGS SHOULD BE INCREASED. ! ALL FUNCTION ARGUMENTS ARE REAL. INTEGER VALUES MAY BE PASSED ! BY ASSIGNING THE INTEGER VALUE TO A REAL ARGUMENT AND THEN USING ! NINT TO RECOVER THE INTEGER VALUE. ! THE AMOUNT OF SPACE AVAILABLE IN THE ARRAY WORK IS ALLOCATED ! DYNAMICALLY. THIS MEANS THAT IT DOES NOT HAVE TO ! BE PROVIDED IN THE CALL TO FUNCTIONS OR IN THE CALL TO EVALUATE_F. ! SUBROUTINES ! PREDEFINED FUNCTIONS : SIN, COS, TAN, ACOS, ATAN, ABS, MAX, NINT ! EXP, LOG, MIN, MOD, SIGN, SQRT, REAL(DBLE) ! PARAMETERS INTEGER, PARAMETER :: alpha = 5, beta = 14, gamma = 3 ! LOCAL DECLARATIONS INTEGER :: i, j, k, n, i1, i2 INTEGER (short) :: ret LOGICAL :: fonly, gonly REAL (stnd) :: pi, biggst !--------- VARIABLES FOR THE TEST FUNCTIONS. REAL (stnd) :: x1, x2, x3, x4, x5, x6, g1, g2, g3, g4, g5, g6, w1, w2, & w3, w4, w5, w6, w7, w8, w9, w10, w11, w12, t, ri, ti, yi, rf1, rf2 ! EQUIVALENCES: None. ! COMMON: None. ! EXECUTION pi = ACOS(-one) biggst = LOG(machhuge) !--------- SET LOGICAL FLAGS AND SELECT FUNCTION. fonly = ifg == justf gonly = ifg == justg n = SIZE(x) ret = ok ! THE TEST FUNCTIONS APPEAR HERE IN ALPHABETICAL ORDER. ! Case( 26 ); 3600 BARD70 ! Case( 11 ); 2100 BIGGS6 ! Case( 16 ); 2600 BOX663 ! Case( 27 ); 3700 CRGLVY ! Case( 25 ); 3500 ENGVL2 ! Case( 40 ); 5000 PENAL1 ! Case( 41 ); 5100 PENAL2 ! Case( 3 ); 1300 PWSING ! Case( 1 ); 1100 ROSENB ! Case( 24 ); 3400 SCHMVT ! Case( 81 ); 9100 SARSEB ! new quadratic - not yet in Mintest ! package !>>>>> NOTE : IF WE SUPPOSE THAT EACH OF THESE TEST FUNCTIONS HAD !>>>>> BEEN CODED AS A SEPARATE ROUTINE, THEN, UNLESS !>>>>> OTHERWISE SPECIFIED, ALL TEST FUNCTIONS WOULD HAVE !>>>>> HAD AN ARGUMENT LIST AS FOLLOWS: !>>>>> !>>>>> ( X, F, G, CASE ) !>>>>> !>>>>> THOSE WHICH WOULD REQUIRE ADDITIONAL ARGUMENTS ARE !>>>>> NOTED BY GIVING A SUITABLE CALLING SEQUENCE. THIS !>>>>> SERVES TO DEFINE THE SPECIAL ARGUMENTS FOR THOSE TEST !>>>>> FUNCTIONS. SEE FOR EXAMPLE PENAL2 AT 5100. SELECT CASE (funcno) !--------- TESTPACK FUNCTION ROSENB CASE (1) ! 1100 x1 = x(1) w1 = one - x1 w2 = x(2) - x1*x1 IF ( .NOT. gonly) f = c100*w2*w2 + w1*w1 IF ( .NOT. fonly) THEN g(1) = -c400*w2*x1 - two*w1 g(2) = c200*w2 END IF !--------- TESTPACK FUNCTION PWSING CASE (3) ! 1300 IF ( .NOT. gonly) f = zero IF (4*(n/4)/=n) THEN IF ( .NOT. fonly) g = zero ELSE DO i = 1, n/4 j = 4*i w1 = x(j-3) w2 = x(j-2) w3 = x(j-1) w4 = x(j) w5 = w1 + ten*w2 w6 = w3 - w4 w2 = w2 - two*w3 w3 = w2*w2*w2 w1 = w1 - w4 w4 = w1*w1*w1 IF ( .NOT. gonly) f = f + w5*w5 + five*w6*w6 + w2*w3 + ten*w1*w4 IF ( .NOT. fonly) THEN g(j-3) = two*w5 + c40*w4 g(j-2) = c20*w5 + four*w3 g(j-1) = ten*w6 - eight*w3 g(j) = -ten*w6 - c40*w4 END IF END DO END IF !--------- TESTPACK FUNCTION BIGGS6 ( X, F, G, IFG, NINT(ARGUMENTS(1))) !--------- NINT(ARGUMENTS(1)) IS M CASE (11) ! 2100 x1 = x(1) x2 = x(2) x3 = x(3) x4 = x(4) x5 = x(5) x6 = x(6) IF ( .NOT. gonly) f = zero IF ( .NOT. fonly) THEN g1 = zero g2 = zero g3 = zero g4 = zero g5 = zero g6 = zero END IF DO i = 1, NINT(arguments(1)) t = REAL(i) ti = t/ten IF (MAX(-t,-ti*four,-ti*x1,-ti*x2,-ti*x5)<=biggst) THEN yi = EXP(-ti) - five*EXP(-t) + three*EXP(-four*ti) w3 = EXP(-ti*x1) w4 = EXP(-ti*x2) w5 = EXP(-ti*x5) ELSE ret = noforg GO TO 10 END IF ri = x3*w3 - x4*w4 + x6*w5 - yi IF ( .NOT. gonly) f = f + ri*ri IF ( .NOT. fonly) THEN w1 = ti*ri g1 = g1 - w3*w1 g2 = g2 + w4*w1 g3 = g3 + w3*ri g4 = g4 - w4*ri g5 = g5 - w5*w1 g6 = g6 + w5*ri END IF END DO IF ( .NOT. fonly) THEN g(1) = two*x3*g1 g(2) = two*x4*g2 g(3) = two*g3 g(4) = two*g4 g(5) = two*x6*g5 g(6) = two*g6 END IF !--------- TESTPACK FUNCTION BOX663 ( X, F, G, IFG, NINT(ARGUMENTS(1))) !--------- NINT(ARGUMENTS(1)) IS M CASE (16) ! 2600 IF ( .NOT. gonly) f = zero IF ( .NOT. fonly) THEN g1 = zero g2 = zero g3 = zero END IF DO i = 1, NINT(arguments(1)) w2 = REAL(i) ti = w2/ten IF (MAX(-w2,-ti,-ti*x(1),-ti*x(2))<=biggst) THEN w3 = EXP(-ti*x(1)) w4 = EXP(-ti*x(2)) w5 = EXP(-ti) - EXP(-w2) ELSE ret = noforg GO TO 10 END IF ri = w3 - w4 - w5*x(3) IF ( .NOT. gonly) THEN IF (ABS(ri)<=SQRT(machhuge-MAX(f,zero))) THEN f = f + ri*ri ELSE ret = noforg GO TO 10 END IF END IF IF ( .NOT. fonly) THEN w2 = ti*ri g1 = g1 - w3*w2 g2 = g2 + w4*w2 g3 = g3 - w5*ri END IF END DO IF ( .NOT. fonly) THEN g(1) = two*g1 g(2) = two*g2 g(3) = two*g3 END IF !--------- TESTPACK FUNCTION SCHMVT CASE (24) ! 3400 x1 = x(1) x2 = x(2) x3 = x(3) w1 = x1 - x2 w2 = x1 + x3 w3 = one + w1*w1 w4 = (pi*x2+x3)/two w5 = (w2/x2) - two IF (-w5**2<=biggst) THEN w6 = EXP(-w5*w5) ELSE ret = noforg GO TO 10 END IF IF ( .NOT. gonly) f = -((one/w3)+SIN(w4)+w6) IF ( .NOT. fonly) THEN w3 = two*w1/(w3*w3) w4 = COS(w4)/two w6 = two*w5*w6/x2 g(1) = w3 + w6 g(2) = -w3 - pi*w4 - w6*w2/x2 g(3) = -w4 + w6 END IF !--------- TESTPACK FUNCTION ENGVL2 CASE (25) ! 3500 x1 = x(1) x2 = x(2) x3 = x(3) w1 = x1*x1 w2 = x1*w1 w3 = x2*x2 w4 = x3*x3 w5 = x3 - two w6 = five*x3 - x1 + one w7 = w1 + w3 - one w8 = w7 + w4 w9 = w7 + w5*w5 w10 = x1 + x2 + x3 - one w11 = x1 + x2 - x3 + one w12 = w2 + three*w3 + w6*w6 - c36 IF ( .NOT. gonly) f = w8*w8 + w9*w9 + w10*w10 + w11*w11 + w12*w12 IF ( .NOT. fonly) THEN w10 = w8 + w9 g(1) = two*(two*x1*w10+two*(x1+x2)+w12*(three*w1-two*w6)) g(2) = two*(two*x2*w10+two*(x1+x2)+six*w12*x2) g(3) = two*(two*(w8*x3+w5*w9)+two*x3-two+ten*w12*w6) END IF !--------- TESTPACK FUNCTION BARD70 CASE (26) ! 3600 x1 = x(1) x2 = x(2) x3 = x(3) IF ( .NOT. gonly) f = zero IF ( .NOT. fonly) THEN g1 = zero g2 = zero g3 = zero END IF DO i = 1, 15 w1 = REAL(i) w2 = REAL(16-i) w3 = MIN(w1,w2) w4 = x2*w2 + x3*w3 ri = bard7y(i) - (x1+w1/w4) w4 = w4*w4 IF ( .NOT. gonly) f = f + ri*ri IF ( .NOT. fonly) THEN w4 = ri*w1/w4 g1 = g1 - ri g2 = g2 + w2*w4 g3 = g3 + w3*w4 END IF END DO IF ( .NOT. fonly) THEN g(1) = g1*two g(2) = g2*two g(3) = g3*two END IF !--------- TESTPACK FUNCTION CRGLVY CASE (27) ! 3700 x1 = x(1) x2 = x(2) x3 = x(3) x4 = x(4) w1 = x2 - x3 w2 = x3 - x4 w3 = x4 - one IF (x1<=biggst) THEN w4 = EXP(x1) ELSE ret = noforg GO TO 10 END IF w5 = w4 - x2 w6 = TAN(w2) IF ( .NOT. gonly) f = w5**4 + c100*w1**6 + w6**4 + x1**8 + w3*w3 IF ( .NOT. fonly) THEN w2 = COS(w2) w5 = four*w5**3 w1 = c600*w1**5 w6 = four*w6**3/(w2*w2) g(1) = w4*w5 + eight*x1**7 g(2) = -w5 + w1 g(3) = -w1 + w6 g(4) = -w6 + two*w3 END IF !--------- TESTPACK FUNCTION PENAL1 ( X, F, G, IFG, ! ARGUMENTS(1), ARGUMENTS(2) ) !--------- ARGUMENTS(1) IS A !--------- ARGUMENTS(2) IS B CASE (40) ! 5000 rf1 = arguments(1) rf2 = arguments(2) w1 = -one/four w2 = zero DO j = 1, n w3 = x(j) w1 = w1 + w3*w3 w3 = w3 - one w2 = w2 + w3*w3 END DO IF ( .NOT. gonly) f = rf1*w2 + rf2*w1*w1 IF ( .NOT. fonly) THEN w1 = four*rf2*w1 w2 = two*rf1 DO j = 1, n w3 = x(j) g(j) = w2*(w3-one) + w3*w1 END DO END IF !--------- TESTPACK FUNCTION PENAL2 ( X, F, G, IFG ) !--------- ARGUMENTS(1) IS A !--------- ARGUMENTS(2) IS B CASE (41) ! 5100 rf1 = arguments(1) rf2 = arguments(2) IF (SIZE(work)<2*n) THEN f = zero g = zero GO TO 10 END IF w1 = EXP(tenth) w2 = EXP(-tenth) w3 = zero i1 = 0 i2 = n DO k = 1, n w4 = x(k) w3 = w3 + REAL(n-k+1)*w4*w4 IF (tenth*w4<=biggst) THEN w5 = EXP(tenth*w4) ELSE ret = noforg GO TO 10 END IF IF (k==1) THEN w6 = zero w7 = one ELSE w7 = w9*w1 w10 = w5 + w8 - (w7+w9) w11 = w5 - w2 IF ( .NOT. fonly) THEN work(i1+k) = w10 work(i2+k) = w11 END IF IF ( .NOT. gonly) w6 = w6 + w10*w10 + w11*w11 END IF w8 = w5 w9 = w7 END DO w1 = x(1) - fifth w2 = w3 - one IF ( .NOT. gonly) f = rf1*w6 + rf2*(w1*w1+w2*w2) IF ( .NOT. fonly) THEN w3 = fifth*rf1 w2 = four*rf2*w2 DO k = 1, n ! --NOTE THAT W8 DOES NOT NEED TO BE PRE-DEFINED WHEN K = 1. w4 = x(k) IF (tenth*w4<=biggst) THEN w5 = EXP(tenth*w4) ELSE ret = noforg GO TO 10 END IF IF (k/=1) THEN w6 = w8 w7 = work(i2+k) END IF IF (k> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> toms.f90 PROGRAM test_minimize_f USE minimize, ONLY : done, long, minimizestate, minimize_f USE min_codes, ONLY : available, ident, normal, revcommfandg, & revcommrestart, sumform USE support, ONLY : evaluate_f USE supp_codes, ONLY : analytic, both USE supp_defs, ONLY : defaultevalstate USE test_functions, ONLY : arguments, funcno, functions, short, stnd USE reals, ONLY : c13, c1_m5, c20, half, one, ten, three, two, zero USE general, ONLY : cpusecs, my_date_time USE true_false, ONLY : false, true USE min_states, ONLY : tracelist ! DESCRIPTION: ! This is a routine provided to test the minimization routine minimize_f ! after installation on a particular system. It calls minimize_f ! to minimize a collection of 10 test functions which are provided ! in the module Test_Functions. ! It also serves as a model to illustrate the use of some of the ! features of the minimization algorithm implementation. ! For an example of the coding of a test function, see the routine ! Functions in the file test_functions.f90. ! Each function is minimized several times. Tests involve analytic ! and differenced derivatives and use both forward and reverse ! communication. Both the conjugate gradient and quasi-Newton codes ! are tried as well as the Nocedal updates. These are the standard ! tests performed with the namelist file provided. A total of 73 ! test runs will be executed. ! Note that test #8 runs a function with n = 500 using the quasi-Newton ! method, if enough memory is available. On a Sun Sparc10 (which runs at ! about 5 megaflops), this function takes about 130 seconds to minimize at ! double precision. On a Sparc1 (about 0.7 megaflops), it takes over ! 660 secs. ! If it is preferred that this long test not be run, simply set the ! single occurrence of the argument run8 in the file namelist.inp ! to false. In this case, only 72 test runs will be made instead of 73. ! The time required on a SPARCstation 1 is then reduced to about 110 secs. !----- Summary of Tests. ! Standard tests: ! 1. Forward call, analytic derivatives, method = Available. ! 2. Reverse call, analytic derivatives, method = Available. ! 3. Forward call, finite differences, method = Available. ! 4. Forward call, derivative testing, method = Available. ! 5. Forward call, analytic derivatives, method = FixTerms. ! 6. Reverse call, analytic derivatives, method = QN. ! 7. Forward call, analytic derivatives, method = Available,Nocedal updates. ! 8. Forward call, analytic derivatives, method = Available, function #11 ! only (dimension large). ! 9. Forward call, analytic derivatives, method = Dynamic, expense = 1, ! function #11 only (dimension large). ! 10. Forward call, analytic derivatives, method = Dynamic, expense = 300, ! function #11 only (dimension large). ! The auxiliary tests listed below provide extra testing, but require a more ! extensive namelist file. They also use a quadratic function that is not ! used in the standard set of tests. The auxiliary tests are *not* used ! in testing the minimization routine after installation on a new system. ! Auxiliary tests: ! 11. Forward call, analytic derivatives, method = ConMin. ! 12. Forward call, analytic derivatives, method = Variable, Nocedal updates. ! 13. Forward call, analytic derivatives, method = SD, ! exact line search with quadratic function sarseb only. ! 14. Forward call, analytic derivatives, method = CG, ! exact line search with quadratic function sarseb only. ! SUBROUTINES: ! minimize_f ! evaluate_f ! My_Date_Time ! CpuSecs ! functions ! bard70, biggs6, box663, crglvy, engvl2 |a subset of a full coll- ! penal1, penal2, pwsing, rosenb, schmvt |ection of test functions ! sarseb ! PARAMETERS: IMPLICIT NONE !-------------! ! maximal dimension ! number of tests ! number of problems INTEGER, PARAMETER :: mxn = 500, tests = 14, nprobs = 12, & ttests = tests*nprobs ! interactive input ! interactive output ! input from data file INTEGER (short), PARAMETER :: inpt = 7, outpt = 8 ! output from tests INTEGER (short):: trminp, trmout INTEGER :: i1mach ! print first/last points INTEGER (long), PARAMETER :: freq = 1000, memory = 130000 ! available memory ! LOCAL DECLARATIONS: !---- Places to hold some statistics and stuff from the tests. INTEGER :: icnts(ttests), fcnts(ttests), fnct, grct, indx(nprobs), & compnt(nprobs) REAL (stnd) :: fvals(ttests), acctim, acc !---- Various declarations needed to run the tests. TYPE (tracelist) :: dotraces INTEGER :: n, error, i, m, contrl, from, to, pdone, tfncs, titers, & mprobs INTEGER, PARAMETER :: DIM(nprobs) = (/ 3, 6, 3, 4, 3, 10, 10, 40, 2, 3, & 500, 4 /), ifnc(nprobs) = (/ 26, 11, 16, 27, 25, 40, 41, 3, 1, 24, 3, & 81 /) REAL (stnd), PARAMETER :: rpar1(nprobs) = (/ zero, c13, ten, zero, zero, & c1_m5, c1_m5, zero, zero, zero, zero, zero/), rpar2(nprobs) = (/ zero, & zero, zero, zero, zero, one, one, zero, zero, zero, zero, zero/) REAL (stnd) :: x(mxn), g(mxn), dererr(nprobs), fx, time, averrs(nprobs) !---- Declarations for remaining Minimize_f arguments. INTEGER (short) :: status, state, case, derv, meth, update, set_h0, & expens INTEGER (long) :: mterms, chkpoint LOGICAL :: dotrace, exact CHARACTER (30) :: chkfile TYPE (minimizestate) :: c LOGICAL :: firsttime, run8 CHARACTER (41) :: date ! Output identification. CHARACTER (60) :: title(tests) = (/ & ' analytic mode, forward calls; meth= Available. ', & ' analytic mode, reverse calls; meth= Available. ', & ' differencing, forward calls; meth= Available. ', & ' testing mode, forward calls; meth= Available. ', & ' analytic mode, forward calls; meth= FixTerms. ', & ' analytic mode, reverse calls; meth= QN. ', & ' analytic mode, forward calls, meth= Available, Nocedal ups.', & ' analytic mode, forward calls; meth= Available, big n. ', & ' analytic mode, forward calls; meth= Dynamic, big n. ', & ' analytic mode, forward calls; meth= Dynamic, big expense,n.', & ' analytic mode, forward calls, meth= ConMin. ', & ' analytic mode, forward calls, meth= Variable, Nocedal ups. ', & ' analytic mode, forward calls, meth= SD, exact line search. ', & ' analytic mode, forward calls, meth= CG, exact line search. '/) NAMELIST /args/state, dotrace, chkpoint, chkfile, exact, expens, derv, & meth, mterms, update, set_h0, run8 ! EXECUTION: acctim = zero !---- Initialize timing. trminp = i1mach(1) trmout = i1mach(2) OPEN (inpt,file='data1.7',action='READ') OPEN (outpt,file='res1.8') CALL my_date_time(date) WRITE (outpt,90020) ' starting test at ', date error = 0 ! Initialize counts. pdone = 0 tfncs = 0 titers = 0 dererr = zero compnt = 0 indx = 0 averrs = zero fvals = zero fcnts = 0 icnts = 0 mprobs = 10 ! Do minimizations. Run each of the test types. to = mprobs contrl = 0 ! Initialize accuracy. IF (stnd==SELECTED_REAL_KIND(12)) THEN acc = 5.0D-04 ! accuracy for normal precision ELSE IF (stnd==SELECTED_REAL_KIND(6)) THEN acc = 4.0D-03 ! accuracy for low precision END IF TESTSET: DO m = 1, tests ! Default settings for some of the arguments of Minimize_f. ! These arguments can be overridden by namelist input: state = normal dotrace = false chkpoint = 0 chkfile = 'ChkPt' exact = false expens = 1 derv = analytic meth = available mterms = 0 update = sumform set_h0 = ident run8 = true READ (inpt,nml=args,end=30) ! read namelist data ! Write title, unless no output requested. WRITE (trmout,'('' Ready for run #'',i2,'':'', a)') m, title(m) IF (freq/=0) WRITE (outpt,90080) m, title(m) ! Here we have a chance to choose a subset of the problems to test. IF ( .NOT. run8) CYCLE ! If run8 = false, do not run test #8. 10 IF (contrl/=-3) THEN contrl = -3 WRITE (trmout,*) ' control: 0 quit' WRITE (trmout,*) ' -1 skip to next set,' WRITE (trmout,*) ' -2 finish this set' WRITE (trmout,*) ' -3 (or eof) finish full run' WRITE (trmout,*) ' n > 0 do problem #n' READ (trminp,'(bn,i2)',end=20) contrl END IF 20 SELECT CASE (contrl) CASE (1:) from = contrl to = contrl CASE (0) GO TO 30 CASE (-1) CYCLE CASE (-2) from = MOD(to,mprobs) + 1 to = mprobs CASE (-3) from = MOD(to,mprobs) + 1 to = mprobs END SELECT ! For tests 8 through 10 just use function 11 (pwsing with n=500). ! For tests 13 and 14 just use function 12 (sarseb). SELECT CASE (m) CASE (8,9,10) from = 11 to = 11 mprobs = 11 CASE (11,12) IF (contrl<1) THEN from = 1 to = 10 END IF mprobs = 10 CASE (13,14) from = nprobs to = nprobs mprobs = nprobs END SELECT time = cpusecs() ! Start timing acctim = acctim - time EACHFUNCTION: DO i = from, to !Repeat for each test function selected. pdone = pdone + 1 IF (pdone>ttests) THEN WRITE (outpt,*) ' too many tests: stopping.' GO TO 40 END IF IF (freq/=0) WRITE (outpt,90030) i n = DIM(i) ! Set dimension, etc. funcno = ifnc(i) arguments(1) = rpar1(i) arguments(2) = rpar2(i) SELECT CASE (i) ! Set starting point. CASE (1) x(1:n) = (/ one, one, one/) CASE (2) x(1:n) = (/ one, two, one, one, one, one/) CASE (3) x(1:n) = (/ zero, ten, c20/) CASE (4) x(1:n) = (/ one, two, two, two/) CASE (5) x(1:n) = (/ one, two, zero/) CASE (6) x(1:n) = (/ (i,i=1,n) /) CASE (7) x(1:n) = half CASE (8) x(1:n) = (/ (three,-one,zero,one,i=1,n/4) /) CASE (9) x(1:n) = (/ -1.2_stnd, one/) CASE (10) x(1:n) = half CASE (11) x(1:n) = (/ (three,-one,zero,one,i=1,n/4) /) CASE (12) x(1:n) = (/ c20, c20, c20, c20/) END SELECT ! ---Set up calls to minimize. IF (dotrace) THEN dotraces = tracelist(6,0,true,true,true,true,true,true,true,true, & true,true,true) ELSE dotraces = tracelist(6,0,false,false,false,false,false,false, & false,false,false,false,false) END IF SELECT CASE (m) CASE (1,3,4,5,7:) ! Forward calls status = state CALL minimize_f(functions,x(1:n),fx,g,acc,status,memory,c, & traces=dotraces,checkpoint=chkpoint,checkfile=chkfile, & exactls=exact,expense=expens,relativetof0=true, & relativetog0=true,usegrad=false,usestep=true,printunit=8_short, & frequency=freq,derivatives=derv,method=meth,terms=mterms, & updateform=update,seth0=set_h0) IF (m==4) THEN dererr(i) = c%evalers%worst averrs(i) = c%evalers%average compnt(i) = c%evalers%index indx(i) = c%evalers%gradcnt END IF CASE (2,6) ! Reverse communication status = state c%eval = defaultevalstate firsttime = true DO case = both CALL evaluate_f(functions,x(1:n),fx,g(1:n),case,c%eval, & c%evalcts,first=firsttime) CALL minimize_f(functions,x(1:n),fx,g,acc,status,memory,c, & traces=dotraces,checkpoint=chkpoint,checkfile=chkfile, & exactls=exact,expense=expens,relativetof0=true, & relativetog0=true,usegrad=false,usestep=true, & printunit=8_short,frequency=freq,derivatives=derv,method=meth, & terms=mterms,updateform=update,seth0=set_h0) IF (status/=revcommfandg) EXIT status = revcommrestart firsttime = false END DO END SELECT IF (status/=done) THEN ! Add number of errors. error = error + 1 END IF fnct = c%evalcts%fevals grct = c%evalcts%gevals fvals(pdone) = fx fcnts(pdone) = c%evalcts%fevals icnts(pdone) = c%shared%it tfncs = tfncs + c%evalcts%fevals titers = titers + c%shared%it IF (freq/=0) WRITE (outpt,90060) status END DO EACHFUNCTION time = cpusecs() acctim = acctim + time IF (m==4) THEN IF (freq/=0) WRITE (outpt,90070) (dererr(i),compnt(i),indx(i), & averrs(i),i=1,mprobs) END IF IF (to/=mprobs) GO TO 10 END DO TESTSET 30 WRITE (outpt,90000) WRITE (outpt,90010) (i,icnts(i),fvals(i),fcnts(i),i=1,pdone) WRITE (outpt,90050) pdone, error, tfncs, titers WRITE (outpt,90040) acctim WRITE (trmout,*) ' ' WRITE (trmout,*) ' ' WRITE (trmout,*) ' ' WRITE (trmout,90040) acctim CALL my_date_time(date) WRITE (outpt,90020) ' test ended at ', date WRITE (trmout,90020) ' test ended at ', date ! EXIT: 40 STOP ! FORMATS: 90000 FORMAT (//' rn its funct. value fns |',' rn its funct. value fns |', & ' rn its funct. value fns'/' ------------------------|', & '-------------------------|','------------------------') 90010 FORMAT ((2(I3,I4,1X,E12.6,I4,' |'),(I3,I4,1X,E12.6,I4))) 90020 FORMAT (4A//) 90030 FORMAT (/' function #',I2/) 90040 FORMAT (' time taken was ',F12.3,' seconds.') 90050 FORMAT (//' *****Test Finished**** problems done ',I3, & '; number of errors is ',I2,'.'/ & ' total function calls = ',I4,' total iterations = ', & I4//) 90060 FORMAT (/' ************* Run Complete, status = ',I3,'.'/) 90070 FORMAT (//' Testing mode derivative estimation errors'// & ' max error component iterate av. decimals '// & 10(E10.2,I7,7X,I5,F9.2/)) 90080 FORMAT ('1 beginning run #',I2,':',A) !END: END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f90' then echo shar: will not over-write existing file "'driver2.f90'" else cat << SHAR_EOF > 'driver2.f90' C--**--CH2470--734--P:RW--22:9:1999 ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> rosenbrockfunc.f90 MODULE RosenbrockFunc USE Precisions, only : stnd, short USE Supp_Codes, only : justf, justg, both, OK USE Reals, only : one, two, c100, c200, c400 CONTAINS SUBROUTINE Rosenbrock ( x, f, g, job ) Implicit None Real(stnd), intent(IN) :: x(:) Real(stnd), intent(OUT) :: f, g(size(x)) Integer(short), intent(INOUT) :: job Real(stnd) :: w1, w2 ! temporaries Logical :: dof, dog dof = job == justf .or. job == both dog = job == justg .or. job == both w1 = x(2) - x(1)**2; w2 = one - x(1) IF ( dof ) then f = c100*w1**2 + w2**2 end if IF ( dog ) then g(1) = -c400 * w1 * x(1) - two*w2 g(2) = c200 * w1 end if job = OK return end Subroutine Rosenbrock end Module RosenbrockFunc ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> example.f90 Program Min_Rosenbrock USE Precisions, only : stnd, short USE Minimize USE RosenbrockFunc USE Reals, only : one Implicit None Integer(short),PARAMETER :: n = 2 Integer(short) :: state Integer :: nout, i1mach Real(stnd), PARAMETER :: acc = 0.001 Real(stnd) :: f, x(n), g(n) Type(MinimizeState) :: C x = (/ -1.2_stnd, one /) ! set the initial guess. state = Normal nout = i1mach(2) Call Minimize_f (Rosenbrock, x, f, g, acc, state, 0, C, & Frequency = 10, EvalLimit = 200 ) IF ( state == Done ) then write(nout,'(1x,a,e16.8,a)') & 'Least function value: ',f, '.' write(nout,'(1x,a,i5,a)') & 'Function evaluation count: ',C%EvalCts%FEvals,'.' write(nout,'(1x,a,i5,a)') & 'Iteration count: ',C%Shared%it,'.' else write(nout,'(1x,a,i5)') & 'Failure: state = ', state end if stop end SHAR_EOF fi # end of overwriting check if test -f 'res.8' then echo shar: will not over-write existing file "'res.8'" else cat << SHAR_EOF > 'res.8' starting test at 2:27 p.m., Wednesday, September 22, 1999 1 beginning run # 1: analytic mode, forward calls; meth= Available. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 9; f= 0.9979563513253989E-02(# 10) ||g||=.46E-02(# 10); 0.000 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.000 s [Prnt] pt 43; f= 0.5655650634086426E-02(# 46) ||g||=.10E-03(# 46); 0.003 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1155940389680377E-03(# 23) ||g||=.46E-03(# 23); 0.001 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 37; f= 0.2598111368329512E-06(# 39) ||g||=.66E-04(# 39); 0.000 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 26; f= 0.1268884276519305E-11(# 30) ||g||=.10E-03(# 30); 0.000 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446886540286596E-04(# 20) ||g||=.29E-04(# 20); 0.001 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947559660561767E-03(# 24) ||g||=.13E-04(# 24); 0.000 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 35; f= 0.1830050294747119E-05(# 36) ||g||=.79E-03(# 36); 0.002 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 35; f= 0.3316606905652118E-10(# 43) ||g||=.26E-03(# 43); 0.003 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 13; f= -2.999999999962793 (# 18) ||g||=.35E-04(# 18); 0.000 s ************* Run Complete, status = 0. 1 beginning run # 2: analytic mode, reverse calls; meth= Available. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 9; f= 0.9979563513253989E-02(# 10) ||g||=.46E-02(# 10); 0.001 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.000 s [Prnt] pt 43; f= 0.5655650634086426E-02(# 46) ||g||=.10E-03(# 46); 0.004 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1155940389680377E-03(# 23) ||g||=.46E-03(# 23); 0.001 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 37; f= 0.2598111368329512E-06(# 39) ||g||=.66E-04(# 39); 0.002 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 26; f= 0.1268884276519305E-11(# 30) ||g||=.10E-03(# 30); 0.001 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446886540286596E-04(# 20) ||g||=.29E-04(# 20); 0.000 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947559660561767E-03(# 24) ||g||=.13E-04(# 24); 0.000 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 35; f= 0.1830050294747119E-05(# 36) ||g||=.79E-03(# 36); 0.000 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 35; f= 0.3316606905652118E-10(# 43) ||g||=.26E-03(# 43); 0.000 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 13; f= -2.999999999962793 (# 18) ||g||=.35E-04(# 18); 0.000 s ************* Run Complete, status = 0. 1 beginning run # 3: differencing, forward calls; meth= Available. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 9; f= 0.9979563487512386E-02(# 10) ||g||=.46E-02(# 10); 0.001 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.001 s [Prnt] pt 43; f= 0.5655650641108877E-02(# 46) ||g||=.10E-03(# 46); 0.018 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1155940142572029E-03(# 23) ||g||=.46E-03(# 23); 0.007 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 37; f= 0.2598207363761166E-06(# 39) ||g||=.66E-04(# 39); 0.002 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 26; f= 0.7350504582619462E-10(# 30) ||g||=.10E-03(# 30); 0.002 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446897745426434E-04(# 20) ||g||=.29E-04(# 20); 0.002 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947559661113872E-03(# 24) ||g||=.13E-04(# 24); 0.006 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 35; f= 0.1845932705826346E-05(# 36) ||g||=.11E-02(# 36); 0.025 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 35; f= 0.5161578294273414E-10(# 43) ||g||=.24E-03(# 43); 0.002 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 13; f= -2.999999999963106 (# 18) ||g||=.35E-04(# 18); 0.001 s ************* Run Complete, status = 0. 1 beginning run # 4: testing mode, forward calls; meth= Available. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.001 s [Prnt] pt 9; f= 0.9979563513253989E-02(# 10) ||g||=.46E-02(# 10); 0.002 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.001 s [Prnt] pt 43; f= 0.5655650634086426E-02(# 46) ||g||=.10E-03(# 46); 0.027 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1155940389680377E-03(# 23) ||g||=.46E-03(# 23); 0.007 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 37; f= 0.2598111368329512E-06(# 39) ||g||=.66E-04(# 39); 0.000 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 26; f= 0.1268884276519305E-11(# 30) ||g||=.10E-03(# 30); 0.001 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446886540286596E-04(# 20) ||g||=.29E-04(# 20); 0.005 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947559660561767E-03(# 24) ||g||=.13E-04(# 24); 0.010 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.001 s [Prnt] pt 35; f= 0.1830050294747119E-05(# 36) ||g||=.79E-03(# 36); 0.048 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 35; f= 0.3316606905652118E-10(# 43) ||g||=.26E-03(# 43); 0.006 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 13; f= -2.999999999962793 (# 18) ||g||=.35E-04(# 18); 0.001 s ************* Run Complete, status = 0. Testing mode derivative estimation errors max error component iterate av. decimals 0.19E-07 2 3 8.13 0.66E-07 4 2 7.84 0.37E-04 2 3 7.54 -0.11E-06 1 2 7.97 0.28E-05 2 24 8.08 0.29E-05 1 2 7.93 0.49E-05 1 4 8.18 0.36E-03 7 2 8.74 0.77E-07 1 18 8.11 0.21E-06 3 2 8.05 1 beginning run # 5: analytic mode, forward calls; meth= FixTerms. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 8; f= 0.9979046504339782E-02(# 9) ||g||=.33E-02(# 9); 0.000 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.000 s [Prnt] pt 72; f= 0.6137461027584269E-02(# 94) ||g||=.33E-02(# 94); 0.007 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1224664189148026E-03(# 25) ||g||=.57E-03(# 25); 0.002 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 24; f= 0.1885171251667609E-05(# 31) ||g||=.21E-03(# 31); 0.001 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 27; f= 0.7521861258578740E-04(# 36) ||g||=.22E-01(# 36); 0.002 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446692822484672E-04(# 20) ||g||=.29E-04(# 20); 0.000 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 16; f= 0.1948184130374005E-03(# 21) ||g||=.39E-05(# 21); 0.000 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 32; f= 0.1888895835532543E-03(# 45) ||g||=.80E-02(# 45); 0.002 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 44; f= 0.2146041211599826E-07(# 63) ||g||=.15E-02(# 63); 0.001 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 12; f= -2.999999998873182 (# 20) ||g||=.50E-04(# 20); 0.001 s ************* Run Complete, status = 0. 1 beginning run # 6: analytic mode, reverse calls; meth= QN. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 9; f= 0.9979563513253989E-02(# 10) ||g||=.46E-02(# 10); 0.000 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.000 s [Prnt] pt 43; f= 0.5655650634086426E-02(# 46) ||g||=.10E-03(# 46); 0.004 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1155940389680377E-03(# 23) ||g||=.46E-03(# 23); 0.001 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 37; f= 0.2598111368329512E-06(# 39) ||g||=.66E-04(# 39); 0.000 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 26; f= 0.1268884276519305E-11(# 30) ||g||=.10E-03(# 30); 0.001 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446886540286596E-04(# 20) ||g||=.29E-04(# 20); 0.001 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947559660561767E-03(# 24) ||g||=.13E-04(# 24); 0.000 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 35; f= 0.1830050294747119E-05(# 36) ||g||=.79E-03(# 36); 0.000 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 35; f= 0.3316606905652118E-10(# 43) ||g||=.26E-03(# 43); 0.000 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 13; f= -2.999999999962793 (# 18) ||g||=.35E-04(# 18); 0.001 s ************* Run Complete, status = 0. 1 beginning run # 7: analytic mode, forward calls, meth= Available, Nocedal ups. function # 1 [Prnt] pt 0; f= 41.68169586167802 (# 1) ||g||=.85E+02(# 1); 0.000 s [Prnt] pt 10; f= 0.9977313147861876E-02(# 11) ||g||=.33E-02(# 11); 0.000 s ************* Run Complete, status = 0. function # 2 [Prnt] pt 0; f= 0.7790700756559704 (# 1) ||g||=.26E+01(# 1); 0.000 s [Prnt] pt 45; f= 0.5689490105560070E-02(# 59) ||g||=.36E-02(# 59); 0.003 s ************* Run Complete, status = 0. function # 3 [Prnt] pt 0; f= 1031.153810609399 (# 1) ||g||=.15E+03(# 1); 0.001 s [Prnt] pt 18; f= 0.1026354148996716E-03(# 25) ||g||=.43E-03(# 25); 0.002 s ************* Run Complete, status = 0. function # 4 [Prnt] pt 0; f= 2.266182511289055 (# 1) ||g||=.12E+02(# 1); 0.000 s [Prnt] pt 40; f= 0.1673668525013668E-06(# 46) ||g||=.44E-04(# 46); 0.002 s ************* Run Complete, status = 0. function # 5 [Prnt] pt 0; f= 629.0000000000000 (# 1) ||g||=.46E+03(# 1); 0.000 s [Prnt] pt 38; f= 0.4834610081576529E-03(# 51) ||g||=.43E-01(# 51); 0.000 s ************* Run Complete, status = 0. function # 6 [Prnt] pt 0; f= 148032.5653500000 (# 1) ||g||=.30E+05(# 1); 0.000 s [Prnt] pt 19; f= 0.7446694567243586E-04(# 20) ||g||=.29E-04(# 20); 0.001 s ************* Run Complete, status = 0. function # 7 [Prnt] pt 0; f= 162.6526819887061 (# 1) ||g||=.50E+03(# 1); 0.000 s [Prnt] pt 18; f= 0.1947769195679623E-03(# 24) ||g||=.60E-05(# 24); 0.000 s ************* Run Complete, status = 0. function # 8 [Prnt] pt 0; f= 2150.000000000000 (# 1) ||g||=.15E+04(# 1); 0.000 s [Prnt] pt 73; f= 0.1977431543973824E-04(# 95) ||g||=.11E-02(# 95); 0.001 s ************* Run Complete, status = 0. function # 9 [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 39; f= 0.1112474255870166E-09(# 57) ||g||=.77E-04(# 57); 0.000 s ************* Run Complete, status = 0. function #10 [Prnt] pt 0; f= -2.860065561048750 (# 1) ||g||=.84E+00(# 1); 0.000 s [Prnt] pt 10; f= -2.999999933497466 (# 18) ||g||=.50E-03(# 18); 0.000 s ************* Run Complete, status = 0. 1 beginning run # 8: analytic mode, forward calls; meth= Available, big n. function #11 [Prnt] pt 0; f= 26875.00000000000 (# 1) ||g||=.51E+04(# 1); 0.000 s [Prnt] pt 39; f= 0.1826181873111509E-04(# 40) ||g||=.10E-02(# 40); 0.006 s ************* Run Complete, status = 0. 1 beginning run # 9: analytic mode, forward calls; meth= Dynamic, big n. function #11 [Prnt] pt 0; f= 26875.00000000000 (# 1) ||g||=.51E+04(# 1); 0.000 s [Prnt] pt 20; f= 0.2851056671899325E-05(# 37) ||g||=.19E-03(# 37); 0.008 s ************* Run Complete, status = 0. 1 beginning run #10: analytic mode, forward calls; meth= Dynamic, big expense,n. function #11 [Prnt] pt 0; f= 26875.00000000000 (# 1) ||g||=.51E+04(# 1); 0.055 s [Prnt] pt 20; f= 0.2851056671899325E-05(# 37) ||g||=.19E-03(# 37); 1.917 s ************* Run Complete, status = 0. rn its funct. value fns | rn its funct. value fns | rn its funct. value fns ------------------------|-------------------------|------------------------ 1 9 0.997956E-02 10 | 2 43 0.565565E-02 46 | 3 18 0.115594E-03 23 4 37 0.259811E-06 39 | 5 26 0.126888E-11 30 | 6 19 0.744689E-04 20 7 18 0.194756E-03 24 | 8 35 0.183005E-05 36 | 9 35 0.331661E-10 43 10 13 -.300000E+01 18 | 11 9 0.997956E-02 10 | 12 43 0.565565E-02 46 13 18 0.115594E-03 23 | 14 37 0.259811E-06 39 | 15 26 0.126888E-11 30 16 19 0.744689E-04 20 | 17 18 0.194756E-03 24 | 18 35 0.183005E-05 36 19 35 0.331661E-10 43 | 20 13 -.300000E+01 18 | 21 9 0.997956E-02 10 22 43 0.565565E-02 46 | 23 18 0.115594E-03 23 | 24 37 0.259821E-06 39 25 26 0.735050E-10 30 | 26 19 0.744690E-04 20 | 27 18 0.194756E-03 24 28 35 0.184593E-05 36 | 29 35 0.516158E-10 43 | 30 13 -.300000E+01 18 31 9 0.997956E-02 10 | 32 43 0.565565E-02 46 | 33 18 0.115594E-03 23 34 37 0.259811E-06 39 | 35 26 0.126888E-11 30 | 36 19 0.744689E-04 20 37 18 0.194756E-03 24 | 38 35 0.183005E-05 36 | 39 35 0.331661E-10 43 40 13 -.300000E+01 18 | 41 8 0.997905E-02 9 | 42 72 0.613746E-02 94 43 18 0.122466E-03 25 | 44 24 0.188517E-05 31 | 45 27 0.752186E-04 36 46 19 0.744669E-04 20 | 47 16 0.194818E-03 21 | 48 32 0.188890E-03 45 49 44 0.214604E-07 63 | 50 12 -.300000E+01 20 | 51 9 0.997956E-02 10 52 43 0.565565E-02 46 | 53 18 0.115594E-03 23 | 54 37 0.259811E-06 39 55 26 0.126888E-11 30 | 56 19 0.744689E-04 20 | 57 18 0.194756E-03 24 58 35 0.183005E-05 36 | 59 35 0.331661E-10 43 | 60 13 -.300000E+01 18 61 10 0.997731E-02 11 | 62 45 0.568949E-02 59 | 63 18 0.102635E-03 25 64 40 0.167367E-06 46 | 65 38 0.483461E-03 51 | 66 19 0.744669E-04 20 67 18 0.194777E-03 24 | 68 73 0.197743E-04 95 | 69 39 0.111247E-09 57 70 10 -.300000E+01 18 | 71 39 0.182618E-04 40 | 72 20 0.285106E-05 37 73 20 0.285106E-05 37 | *****Test Finished**** problems done 73; number of errors is 0. total function calls = 2329 total iterations = 1926 time taken was 20.146 seconds. test ended at 2:28 p.m., Wednesday, September 22, 1999 SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << SHAR_EOF > 'res1' Ready for run # 1: analytic mode, forward calls; meth= Available. control: 0 quit -1 skip to next set, -2 finish this set -3 (or eof) finish full run n > 0 do problem #n Ready for run # 2: analytic mode, reverse calls; meth= Available. Ready for run # 3: differencing, forward calls; meth= Available. Ready for run # 4: testing mode, forward calls; meth= Available. Ready for run # 5: analytic mode, forward calls; meth= FixTerms. Ready for run # 6: analytic mode, reverse calls; meth= QN. Ready for run # 7: analytic mode, forward calls, meth= Available, Nocedal ups. Ready for run # 8: analytic mode, forward calls; meth= Available, big n. Ready for run # 9: analytic mode, forward calls; meth= Dynamic, big n. Ready for run #10: analytic mode, forward calls; meth= Dynamic, big expense,n. time taken was 20.146 seconds. test ended at 2:28 p.m., Wednesday, September 22, 1999 SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << SHAR_EOF > 'res2' [Prnt] pt 0; f= 24.19999999999999 (# 1) ||g||=.23E+03(# 1); 0.000 s [Prnt] pt 10; f= 1.866316831898981 (# 15) ||g||=.27E+01(# 15); 0.000 s [Prnt] pt 20; f= 0.4875838224001645 (# 32) ||g||=.12E+01(# 32); 0.000 s [Prnt] pt 30; f= 0.1136216817825096E-03(# 51) ||g||=.47E+00(# 51); 0.001 s [Prnt] pt 36; f= 0.2005359678709316E-08(# 59) ||g||=.69E-03(# 59); 0.001 s Least function value: 0.20053597E-08. Function evaluation count: 59. Iteration count: 36. SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f90' then echo shar: will not over-write existing file "'src.f90'" else cat << SHAR_EOF > 'src.f90' ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> true_false.f90 MODULE true_false ! some convenient true false abbreviations IMPLICIT NONE PUBLIC !----- -------------! LOGICAL, PARAMETER :: true = .TRUE., false = .FALSE. END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> low.f90 MODULE low ! precision specification for real computations PUBLIC :: stnd, extd ! This requests the processor to use a real implementation 'stnd' ! which provides at least 6 decimal digits of precision and an ! exponent range of at least 10 ^ +- 35. This would be suitable for ! low accuracy computations. It is expected that this precision will ! be available on all machines. INTEGER, PARAMETER :: stnd = SELECTED_REAL_KIND(6,35) !------------------------- ! A few computations are preferably done in higher precision 'extd'. The ! numbers chosen here should be such that the underlying hardware will ! select a higher precision for kind 'extd' than for kind 'stnd', if ! this is feasible. If a higher precision is not readily available, ! the same values may be used as are given above for 'stnd'. It is ! anticipated that on most machines this higher precision will also ! be available. INTEGER, PARAMETER :: extd = SELECTED_REAL_KIND(12,35) !------------------------- END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> normal.f90 MODULE normal ! precision specification for real computations PUBLIC :: stnd, extd ! This requests the processor to use a real implementation 'stnd' ! which provides at least 12 decimal digits of precision and an ! exponent range of at least 10 ^ +- 50. It is expected that this ! precision will be available on all machines. INTEGER, PARAMETER :: stnd = SELECTED_REAL_KIND(12,50) !------------------------- ! A few computations are preferably done in higher precision 'extd'. The ! numbers chosen here should be such that the underlying hardware will ! select a higher precision for kind 'extd' than for kind 'stnd', if ! this is feasible. If a higher precision is not readily available, ! the same values may be used as are given above for 'stnd'. It is ! anticipated that on many machines this higher precision may ! not be available. !Integer, PARAMETER :: extd = Selected_Real_Kind ( 20, 50 ) ! preferred INTEGER, PARAMETER :: extd = SELECTED_REAL_KIND(12,50) ! NAG f90: Sun 4 !------------------------- END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> extended.f90 MODULE extended ! precision specification for real computations PUBLIC :: stnd, extd ! This requests the processor to use a real implementation 'stnd' ! which provides at least 20 decimal digits of precision and an ! exponent range of at least 10 ^ +- 80. It is expected that this ! precision may not be available on all machines. INTEGER, PARAMETER :: stnd = SELECTED_REAL_KIND(20,80) !------------------------- ! A few computations are preferably done in higher precision 'extd'. The ! numbers chosen here should be such that the underlying hardware will ! select a higher precision for kind 'extd' than for kind 'stnd', if ! this is feasible. If a higher precision is not readily available, ! the same values may be used as are given above for 'stnd'. It is ! anticipated that on many machines, such an even higher precision may ! not be available. INTEGER, PARAMETER :: extd = SELECTED_REAL_KIND(30,80) !------------------------- END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> long.f90 MODULE long ! precision specification for real computations TYPE :: verylong ! not yet implemented PRIVATE INTEGER :: value TYPE (verylong), POINTER :: next END TYPE END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> integers.f90 MODULE integers ! precision specification for integer computations PUBLIC :: short, long ! This is provided for those machines where using short integers ! has some advantage. The range here is from -999 to +999. Note that ! 999 = 10^3 -1. No harm will be done if short integers are made ! the same as long integers. INTEGER, PARAMETER :: short = SELECTED_INT_KIND(3) !------------------------- ! The range here is at least from -9 999 999 to +9 999 999. Note ! that 10^7 - 1 = 9,999,999. This may limit the largest possible value ! of the dimension n of problems which can be solved. If n is to ! be larger, the 7 should be replaced with a number k so that ! n is considerably less than 10^k. INTEGER, PARAMETER :: long = SELECTED_INT_KIND(7) !------------------------- END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> precision.f90 MODULE precisions USE integers, ONLY : short, long USE normal, ONLY : stnd, extd ! This provides a convenient way of selecting the precision ! required for a computation. By simply ensuring that a leading '!' ! appears on all but exactly one of the following USE statements, ! and then recompiling all routines, the precision of an entire ! computation can be altered. ! USE Low ! USE Extended ! USE Long PRIVATE PUBLIC :: stnd, extd, short, long END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> num_constants.f90 MODULE num_constants ! Values of machine numeric characteristics ! This provides simple names for the various machine dependent ! constants with intrinsic values in Fortran 90. All are for ! precision 'stnd'. USE precisions, ONLY : stnd IMPLICIT NONE PUBLIC !----- -------------! PRIVATE :: stnd, x REAL (stnd) :: x !dummy value REAL (stnd), PARAMETER :: machtol = EPSILON(x), machhuge = HUGE(x), & machmaxexp = MAXEXPONENT(x), machminexp = MINEXPONENT(x), & machdecprec = PRECISION(x), machbase = RADIX(x), & machdecexpr = RANGE(x), machtiny = TINY(x) END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> reals.f90 MODULE reals ! Real constants at a given precision ! This provides names for all required real values. By only ! using real values as defined within this module, all ! problems associated with the precision of real literal ! values can be totally avoided. ! All values are at 'stnd' precision. There are values for the ! integral values to 100, the fractions in tenths from 0 to 1, ! the reciprocals of the integers to 10, some common ! fractions, and whatever miscellaneous real values are ! necessary. A standard naming convention is used. USE precisions, ONLY : stnd IMPLICIT NONE PUBLIC !----- -------------! PRIVATE :: stnd ! digits REAL (stnd), PARAMETER :: zero = 0, one = 1, two = 2, three = 3, & four = 4, five = 5, six = 6, seven = 7, eight = 8, nine = 9 ! tenths REAL (stnd), PARAMETER :: c0_1 = 0.1_stnd, c0_2 = 0.2_stnd, & c0_3 = 0.3_stnd, c0_4 = 0.4_stnd, c0_5 = 0.5_stnd, c0_6 = 0.6_stnd, & c0_7 = 0.7_stnd, c0_8 = 0.8_stnd, c0_9 = 0.9_stnd ! reciprocals REAL (stnd), PARAMETER :: tenth = 0.1_stnd, ninth = one/nine, & eighth = 0.125_stnd, seventh = one/seven, sixth = one/six, & fifth = 0.2_stnd, quarter = 0.25_stnd, third = one/three, & half = 0.5_stnd ! fractions a/b named as fa_b REAL (stnd), PARAMETER :: f1_29 = one/29, f2_3 = two/three, & f4_3 = four/three, f7_3 = seven/three ! integral values to 99 REAL (stnd), PARAMETER :: ten = 10, c10 = 10, c40 = 40, c70 = 70, & c11 = 11, c41 = 41, c71 = 71, c12 = 12, c42 = 42, c72 = 72, c13 = 13, & c43 = 43, c73 = 73, c14 = 14, c44 = 44, c74 = 74, c15 = 15, c45 = 45, & c75 = 75, c16 = 16, c46 = 46, c76 = 76, c17 = 17, c47 = 47, c77 = 77, & c18 = 18, c48 = 48, c78 = 78, c19 = 19, c49 = 49, c79 = 79, c20 = 20, & c50 = 50, c80 = 80, c21 = 21, c51 = 51, c81 = 81, c22 = 22, c52 = 52, & c82 = 82, c23 = 23, c53 = 53, c83 = 83, c24 = 24, c54 = 54, c84 = 84, & c25 = 25, c55 = 55, c85 = 85, c26 = 26, c56 = 56, c86 = 86, c27 = 27, & c57 = 57, c87 = 87, c28 = 28, c58 = 58, c88 = 88, c29 = 29, c59 = 59, & c89 = 89, c30 = 30, c60 = 60, c90 = 90, c31 = 31, c61 = 61, c91 = 91, & c32 = 32, c62 = 62, c92 = 92, c33 = 33, c63 = 63, c93 = 93, c34 = 34, & c64 = 64, c94 = 94, c35 = 35, c65 = 65, c95 = 95, c36 = 36, c66 = 66, & c96 = 96, c37 = 37, c67 = 67, c97 = 97, c38 = 38, c68 = 68, c98 = 98, & c39 = 39, c69 = 69, c99 = 99 ! miscellaneous integral values REAL (stnd), PARAMETER :: c100 = 100, c180 = 180, c200 = 200, & c256 = 256, c360 = 360, c400 = 400, c600 = 600, c681 = 681, & c991 = 991, c1162 = 1162, c2324 = 2324, c10000 = 10000, c40000 = 40000 ! miscellaneous real values ! form: d.dd named as cd_dd ! nn ! form: d.dd x 10 named as cd_ddEnn ! -nn ! form: d.dd x 10 named as cd_ddMnn REAL (stnd), PARAMETER :: c1_e6 = 1.0E6_stnd, c2_m6 = 2.0E-6_stnd, & c1_m13 = 1.0E-13_stnd, c1_m5 = 1.0E-5_stnd, c1_m4 = 1.0E-4_stnd, & c4_m2 = 4.0E-2_stnd, c1_m2 = 1.0E-2_stnd, c0_1136 = 0.1136_stnd, & c1_0001 = 1.0001_stnd, c1_2 = 1.2_stnd, c1_5 = 1.5_stnd, & c2_25 = 2.25_stnd, c2_5 = 2.5_stnd, c2_625 = 2.625_stnd, & c7_5 = 7.5_stnd, c10_1 = 10.1_stnd, c19_8 = 19.8_stnd, & c20_2 = 20.2_stnd END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> strings.f90 MODULE strings ! Routines for simple string stuff USE precisions, ONLY : short IMPLICIT NONE PUBLIC !------ -------------! ! Names for common characters CHARACTER (1), PARAMETER :: ampersand = '&', apostrophe = '''', & atsign = '@', backslash = '\\', backquote = '`', bang = '!', & blank = ' ', caret = '^', cbrace = '}', cbracket = ']', cparen = ')', & colon = ':', comma = ',', dash = '-', dollar = '$', equals = '=', & exclamation = '!', greaterthan = '>', hash = '#', lessthan = '<', & minus = '-', obrace = '{', obracket = '[', oparen = '(', & percent = '%', period = '.', plus = '+', quesmark = '?', quote = '"', & semicolon = ';', slash = '/', star = '*', tilde = '~', vertbar = '|', & underscore = '_' CHARACTER (0), PARAMETER :: null = '' ! Codes for case conversions INTEGER (short), PARAMETER :: toupper = 1, tolower = 2, capitalize = 3 ! conversion INTEGER (short), PARAMETER, PRIVATE :: shift = ICHAR('a') - ICHAR('A') ! This assumes that the relation between a lower case letter and ! its corresponding upper case letter is the same for every letter. CONTAINS SUBROUTINE ascii_case_change(string,type) IMPLICIT NONE !-------------! ! ARGUMENTS: CHARACTER (*), INTENT (INOUT) :: string INTEGER (short), INTENT (IN) :: type ! DESCRIPTION: ! This converts each lower case alphabetic letter in string to upper ! case, or vice versa. Specifically, ! If type = ToUpper, conversion is lower to upper ! If type = ToLower, conversion is upper to lower ! If type = Capitalize, use upper for first letter; lower for rest ! Definitions of ToUpper, ToLower and Capitalize may be obtained from ! the host module Strings. ! All non-alphabetic characters are left unchanged. ! It uses the ASCII character set. ! LOCAL DECLARATIONS: INTEGER (short) :: i ! EXECUTION: DO i = 1, LEN(string) SELECT CASE (type) CASE (toupper) SELECT CASE (string(i:i)) CASE ('a':'z') string(i:i) = ACHAR(IACHAR(string(i:i))+shift) END SELECT CASE (tolower,capitalize) SELECT CASE (string(i:i)) CASE ('A':'Z') string(i:i) = ACHAR(IACHAR(string(i:i))-shift) END SELECT END SELECT END DO IF (type==capitalize) THEN SELECT CASE (string(1:1)) CASE ('a':'z') string(1:1) = ACHAR(IACHAR(string(1:1))+shift) END SELECT END IF ! EXIT: RETURN ! FORMATS: none. END SUBROUTINE SUBROUTINE case_change(string,type) IMPLICIT NONE !-------------! ! ARGUMENTS: CHARACTER (*), INTENT (INOUT) :: string INTEGER (short), INTENT (IN) :: type ! DESCRIPTION: ! This converts each lower case alphabetic letter in string to upper ! case, or vice versa. Specifically, ! If type = ToUpper, conversion is lower to upper ! If type = ToLower, conversion is upper to lower ! If type = Capitalize, use upper for first letter; lower for rest ! Definitions of ToUpper, ToLower and Capitalize may be obtained from ! the host module Strings. ! All non-alphabetic characters are left unchanged. ! It uses the underlying machine character set. ! LOCAL DECLARATIONS: INTEGER (short) :: i ! EXECUTION: DO i = 1, LEN(string) SELECT CASE (type) CASE (toupper) SELECT CASE (string(i:i)) CASE ('a':'i') string(i:i) = CHAR(ICHAR(string(i:i))+shift) CASE ('j':'r') string(i:i) = CHAR(ICHAR(string(i:i))+shift) CASE ('s':'z') string(i:i) = CHAR(ICHAR(string(i:i))+shift) END SELECT CASE (tolower,capitalize) SELECT CASE (string(i:i)) CASE ('A':'I') string(i:i) = CHAR(ICHAR(string(i:i))-shift) CASE ('J':'R') string(i:i) = CHAR(ICHAR(string(i:i))-shift) CASE ('S':'Z') string(i:i) = CHAR(ICHAR(string(i:i))-shift) END SELECT END SELECT END DO IF (type==capitalize) THEN SELECT CASE (string(1:1)) CASE ('a':'i') string(1:1) = CHAR(ICHAR(string(1:1))+shift) CASE ('j':'r') string(1:1) = CHAR(ICHAR(string(1:1))+shift) CASE ('s':'z') string(1:1) = CHAR(ICHAR(string(1:1))+shift) END SELECT END IF ! EXIT: RETURN ! FORMATS: none. END SUBROUTINE SUBROUTINE mid_shift(string,from,to,number) IMPLICIT NONE !-------------! ! ARGUMENTS: CHARACTER (*), INTENT (INOUT) :: string INTEGER (short), INTENT (IN) :: from, to, number ! DESCRIPTION: ! This routine performs a shift of characters within string. The ! number of characters shifted is number and they are shifted so ! that the character in position 'from' is moved to position 'to'. ! Characters in the to position are overwritten. Blanks replace ! characters in the from position. Shifting may be left or right, ! and the from and to positions may overlap. Care is taken not to ! alter or use any characters beyond the defined limits of the string. ! LOCAL DECLARATIONS: INTEGER (short) :: end, end1, n, shorten, slen CHARACTER (len=number) :: substring ! EXECUTION: slen = LEN(string) IF (from/=to) THEN end1 = from + number - 1 ! end1 is the initial position of the last ! character in the substring specified to ! be shifted. ! n is the number of characters that will be shifted and that will ! remain within the confines of string. n may have to be reduced ! for a shift to the right to ensure that characters beyond the ! end of string are not dealt with (this is not necessary for left ! shifts since the substring moved to the left can never extend ! beyond the end of string). IF (end1>slen) THEN shorten = end1 - slen n = number - shorten ELSE n = number END IF IF (fromleft) THEN string(left:) = ADJUSTL(string(left:)) ! move left ELSE string(1:left+1) = ADJUSTR(string(1:left+1)) ! move right END IF END IF ! EXIT: RETURN ! FORMATS: None. END SUBROUTINE END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> general.f90 MODULE general ! General purpose routines USE strings, ONLY : blank, mid_shift USE precisions, ONLY : stnd, long, short USE reals, ONLY : one !!!!USE F90_UNIX ! for use with the second version of CpuSecs below IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: entering, leaving, indent, cpusecs, my_date_time, writemat, & writevec INTEGER (short), PARAMETER :: spaces = 2 INTEGER (short) :: lastcount = 0 CHARACTER (9), PARAMETER :: months(12) = (/ 'January ', 'February ', & 'March ', 'April ', 'May ', 'June ', 'July ', & 'August ', 'September', 'October ', 'November ', 'December '/), & days(0:6) = (/ 'Sunday ', 'Monday ', 'Tuesday ', 'Wednesday', & 'Thursday ', 'Friday ', 'Saturday '/) INTEGER (short), PARAMETER :: defunit = 6 ! For WriteMat and WriteVec INTEGER (short), PARAMETER :: defline = 80 ! Default values INTEGER (short), PARAMETER :: defindent = 0 CHARACTER (1), PARAMETER :: deff = 'f' INTEGER (short), PARAMETER :: defw = 12 INTEGER (short), PARAMETER :: defd = 6 INTEGER (short), PARAMETER :: defs = 3 CONTAINS FUNCTION entering(string,trace,unit,level) ! Upon entering a procedure, this function will be called. ! It will return a prefix string suitable for indenting ! output lines from the procedure. It takes the given string ! and prepends 'level' blanks, followed by a '[', and appends ! the character ']'. For example, if string were 'hi' and level ! were 7, it would return ' [hi]'. Level is then also ! incremented by the value of Spaces. ! If trace is true, it also outputs a message that the ! routine identified by string was entered. ! ARGUMENTS: CHARACTER (*), INTENT (IN) :: string LOGICAL, INTENT (IN) :: trace INTEGER (short), INTENT (IN) :: unit INTEGER (short), INTENT (INOUT) :: level ! LOCAL DECLARATIONS: CHARACTER (len=LEN(string)+2+level) :: entering ! EXECUTION: entering(1:level) = blank IF (trace) WRITE (unit,'(A)') entering(1:level) // '[Into ' // & string // ']' entering(level+1:) = '[' // string // ']' level = level + spaces END FUNCTION SUBROUTINE leaving(string,trace,unit,level) ! This is the 'opposite' to Entering. It should be called ! just before leaving a routine. Level is reduced by ! Spaces and if trace is true, an exit message is output. ! ARGUMENTS: CHARACTER (*), INTENT (IN) :: string LOGICAL, INTENT (IN) :: trace INTEGER (short), INTENT (IN) :: unit INTEGER (short), INTENT (INOUT) :: level ! EXECUTION: level = level - spaces IF (trace) WRITE (unit,'(A)') REPEAT(blank,INT(level)) // '[Done ' // & string // ']' END SUBROUTINE SUBROUTINE indent(id,unit,level) ! This is also used to indent output, albeit in a manner ! different from Entering and Leaving. It simply writes ! out 'level' blanks followed by the string id in [], and ! leaves the output file marker where it is. It uses ! nonadvancing output. If level is not present, just ! the id part is output; i.e. level is treated as zero. ! ARGUMENTS: CHARACTER (*), INTENT (IN) :: id INTEGER (short), INTENT (IN) :: unit INTEGER (short), INTENT (IN), OPTIONAL :: level ! LOCAL DECLARATIONS: INTEGER (short) :: lev ! EXECUTION: IF (PRESENT(level)) THEN lev = level ELSE lev = 0 END IF WRITE (unit,'(A)',advance='no') REPEAT(blank,INT(lev)) // '[' // id // & ']' END SUBROUTINE !!!!FUNCTION CpuSecs() ! This function obtains, from a C system routine cputime, ! the current value of the system CPU usage clock. This value (in seconds) ! is then returned as a standard precision real value. !!!! Real(stnd) :: CpuSecs, cputime !!!! CpuSecs = cputime() !!!! return !!!!end Function CpuSecs !!!!FUNCTION CpuSecs() ! This function obtains, from the Unix module F90_UNIX, ! the current value of the system CPU usage clock. This value (in seconds) ! is then returned as a standard precision real value. ! The module F90_UNIX, provided with NAG compiler v2.0, contains the ! functions times and clock_ticks_per_second, as well as the defined type ! TMS. The USE statement for the F90_UNIX module, located at line 6 of ! general.f90, must be uncommented when using this version of CpuSecs. !!!! Real(stnd) :: CpuSecs !!!! Type(TMS) :: buffer !!!! CpuSecs = times(buffer) !!!! CpuSecs = real((buffer%utime+buffer%stime),stnd) / & !!!! clock_ticks_per_second() !!!! return !!!!end Function CpuSecs FUNCTION cpusecs() ! This function obtains, from the intrinsic routine system_clock, ! the current value of the system CPU usage clock. This value ! is then converted to seconds and returned as a standard precision ! real value. ! LOCAL DECLARATIONS: REAL (stnd) :: cpusecs INTEGER :: count, count_rate, count_max REAL (stnd) :: secs ! EXECUTION: CALL SYSTEM_CLOCK(count,count_rate,count_max) secs = REAL(count,stnd)/REAL(count_rate,stnd) ! wraparound of clock ticker IF (count=13) THEN k = k - 12 tdate(ptampm:ptampm) = 'p' ELSE IF (k==12) THEN tdate(ptampm:ptampm) = 'p' ELSE IF (k==0) THEN k = k + 12 tdate(ptampm:ptampm) = 'a' ELSE tdate(ptampm:ptampm) = 'a' END IF WRITE (tdate(pthour:pthour+1),'(i2)') k tdate(ptdayn:ptdayn+8) = days(dayno) k = LEN_TRIM(days(dayno)) IF (k/=9) THEN ! ==> shift over blanks. CALL mid_shift(tdate,ptdayn+9_short,ptdayn+k,modlen) END IF modlen = MIN(modlen,INT(LEN_TRIM(chdate),short)) chdate(1:modlen) = tdate RETURN END SUBROUTINE ! Optional SUBROUTINE writemat(x,a,b,y,f,w,d,s,unit,name,indent,line) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:,:) REAL (stnd), INTENT (IN), OPTIONAL :: y(SIZE(x,1),SIZE(x,2)), a, b INTEGER (short), INTENT (IN), OPTIONAL :: unit, w, d, s, indent, line CHARACTER (1), INTENT (IN), OPTIONAL :: f CHARACTER (*), INTENT (IN), OPTIONAL :: name ! DESCRIPTION: ! Print out a submatrix of a matrix with given format, as below. ! Print a title for the matrix: name ! If the value of the argument indent is positive, then each output ! line is preceded by indent blank characters. If indent is negative, ! then only the matrix output will be indented (not the title). ! Print the values from the matrix a*X + b*Y. ! a and b are scalars; X and Y are matrices. ! Print each entry in format fw.d ! ..f is a character 'f', 'g', 'e' or 'd' ! ..w and d are integers ! ..s is the number of spaces between each entry. ! ..line is the number of characters per line. ! (if line=0, then line is replaced with 80) ! All output is on the unit 'unit'. ! Defaults are defined for all optional arguments. See start of module. ! LOCAL DECLARATIONS: REAL (stnd) :: fa, fb INTEGER (short) :: actunit, actw, actd, acts, actindent, actline, i, j CHARACTER (1) :: actf CHARACTER (34), PARAMETER :: defform = & '((??x,??(? ??.??, ?? x), ? ??.??))' ! e.g. ((02x,05(f 09.04, 03 x), f 09.04)) ! 1234567890123456789012345678901234 CHARACTER (34) :: form ! EXECUTION: ! First set appropriate default values. form = defform IF (PRESENT(unit)) THEN actunit = unit ELSE actunit = defunit END IF IF (PRESENT(line)) THEN actline = line ELSE actline = defline END IF IF (PRESENT(indent)) THEN actindent = indent ELSE actindent = defindent END IF IF (PRESENT(f)) THEN actf = f ELSE actf = deff END IF IF (PRESENT(w)) THEN actw = w ELSE actw = defw END IF IF (PRESENT(d)) THEN actd = d ELSE actd = defd END IF IF (PRESENT(s)) THEN acts = s ELSE acts = defs END IF IF (PRESENT(a)) THEN fa = a ELSE fa = one END IF IF (PRESENT(b)) THEN fb = b ELSE fb = one END IF IF (actindent==0) THEN WRITE (form(03:06),90010) blank ELSE WRITE (form(03:04),90020) ABS(actindent) END IF WRITE (form(07:08),90020) (actline-ABS(actindent)-actw)/(actw+acts) WRITE (form(10:10),90000) actf WRITE (form(12:13),90020) actw WRITE (form(15:16),90020) actd IF (acts/=0) THEN WRITE (form(19:20),90020) acts ELSE WRITE (form(17:22),90010) blank END IF WRITE (form(26:26),90000) actf WRITE (form(28:29),90020) actw WRITE (form(31:32),90020) actd DO i = 1, actindent WRITE (unit,90010,advance='no') blank END DO IF (PRESENT(name)) WRITE (unit,90010) name DO i = 1, SIZE(x,1) IF (PRESENT(y)) THEN WRITE (unit,form) (fa*x(i,j)+fb*y(i,j),j=1,SIZE(x,2)) ELSE WRITE (unit,form) (fa*x(i,j),j=1,SIZE(x,2)) END IF END DO RETURN ! FORMATS: Also see character string form. 90000 FORMAT (A1) 90010 FORMAT (A) 90020 FORMAT (I2) END SUBROUTINE ! Optional SUBROUTINE writevec(x,a,b,y,f,w,d,s,unit,name,indent,line) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:) REAL (stnd), INTENT (IN), OPTIONAL :: y(SIZE(x)), a, b INTEGER (short), INTENT (IN), OPTIONAL :: unit, w, d, s, indent, line CHARACTER (1), INTENT (IN), OPTIONAL :: f CHARACTER (*), INTENT (IN), OPTIONAL :: name ! DESCRIPTION: This functions just as WriteMat, except that it ! prints a vector instead of a matrix. ! LOCAL DECLARATIONS: REAL (stnd) :: fa, fb INTEGER (short) :: actunit, actw, actd, acts, actindent, actline, i CHARACTER (1) :: actf CHARACTER (34), PARAMETER :: defform = & '((??x,??(? ??.??, ?? x), ? ??.??))' ! e.g. ((02x,05(f 09.04, 03 x), f 09.04)) ! 1234567890123456789012345678901234 CHARACTER (34) :: form ! EXECUTION: ! First set appropriate default values. form = defform IF (PRESENT(unit)) THEN actunit = unit ELSE actunit = defunit END IF IF (PRESENT(line)) THEN actline = line ELSE actline = defline END IF IF (PRESENT(indent)) THEN actindent = indent ELSE actindent = defindent END IF IF (PRESENT(f)) THEN actf = f ELSE actf = deff END IF IF (PRESENT(w)) THEN actw = w ELSE actw = defw END IF IF (PRESENT(d)) THEN actd = d ELSE actd = defd END IF IF (PRESENT(s)) THEN acts = s ELSE acts = defs END IF IF (PRESENT(a)) THEN fa = a ELSE fa = one END IF IF (PRESENT(b)) THEN fb = b ELSE fb = one END IF IF (actindent==0) THEN WRITE (form(03:06),90010) blank ELSE WRITE (form(03:04),90020) ABS(actindent) END IF WRITE (form(07:08),90020) (actline-ABS(actindent)-actw)/(actw+acts) WRITE (form(10:10),90000) actf WRITE (form(12:13),90020) actw WRITE (form(15:16),90020) actd IF (acts/=0) THEN WRITE (form(19:20),90020) acts ELSE WRITE (form(17:22),90010) blank END IF WRITE (form(26:26),90000) actf WRITE (form(28:29),90020) actw WRITE (form(31:32),90020) actd DO i = 1, actindent WRITE (unit,90010,advance='no') blank END DO IF (PRESENT(name)) WRITE (unit,90010) name IF (PRESENT(y)) THEN WRITE (unit,form) (fa*x(i)+fb*y(i),i=1,SIZE(x)) ELSE WRITE (unit,form) (fa*x(i),i=1,SIZE(x)) END IF RETURN ! FORMATS: Also see character string form. 90000 FORMAT (A1) 90010 FORMAT (A) 90020 FORMAT (I2) END SUBROUTINE END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> inner_product.f90 MODULE inner_product USE precisions, ONLY : stnd, long, short ! This contains routines for computing the normal inner product of ! two vectors, as given by their dot product, as well as a routine ! for doing a generalized inner product. Without modification, the ! generalized inner product is a duplicate of the normal one. ! An operator is provided for applying each inner product, namely ! x .GIP. y for the generalized inner product of x and y ! x .IP. y for the normal inner product of x and y ! The operator .GIP. is in fact overloaded; there must be a second, ! specialized version which computes the generalized inner product ! of the vector x with the vector e[j] = (0,0,...,0,1,0,...,0), i.e. ! with the canonical jth unit vector. It is accessed as ! x .GIP. j ! For each inner product, the norm ||.|| it induces is also available, ! as defined by ! 1/2 ! .GNORM. x = ||x|| = (x .GIP. x) for the generalized norm ! 1/2 ! .NORM. x = ||x|| = (x .IP. x) for the usual norm ! In the case of the usual norm, a special routine is used for ! computing the value with due regard to avoiding unnecessary overflow. !---------------------------------------------------------------------- ! Users who wish may replace the general inner product routine with ! one of their own. In this case, the interface should NOT be changed. ! Only the code WITHIN the functions GenInner, GenInnerJ and GenNorm ! should be modified; i.e. all changes should be restricted to these ! functions. ! It is up to the user to verify that a legal inner product is used. ! Users who do not understand this discussion may safely ignore this; ! anyone who may need this capability will undoubtedly know what this ! is all about. ! If the general inner product routine is altered, the general norm ! function must also be changed, in accordance with the definition above. ! This may be done by activating the line so marked in the routine GenNorm. ! If overflow is likely to be a problem when the generalized norm is computed ! by taking the square root of the inner product, a special general norm ! routine may be written. It is the responsibility of the user. !---------------------------------------------------------------------- IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: OPERATOR (.ip.), OPERATOR (.gip.), OPERATOR (.norm.), & OPERATOR (.gnorm.) INTERFACE OPERATOR (.gip.) MODULE PROCEDURE geninner, geninnerj END INTERFACE INTERFACE OPERATOR (.gnorm.) MODULE PROCEDURE gennorm END INTERFACE INTERFACE OPERATOR (.ip.) MODULE PROCEDURE inner END INTERFACE INTERFACE OPERATOR (.norm.) MODULE PROCEDURE norm2 END INTERFACE CONTAINS REAL (stnd) FUNCTION geninner(x,y) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:) REAL (stnd), INTENT (IN) :: y(SIZE(x)) ! DESCRIPTION: ! This is made available for users wishing to redefine the ! underlying metric of the space involved by redefining the ! inner product used for the updating of the quasi-Newton ! matrices. ! EXECUTION: ! Replace code below this line if required--------------- geninner = DOT_PRODUCT(x,y) ! Replace code above this line if required--------------- END FUNCTION REAL (stnd) FUNCTION geninnerj(x,j) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:) INTEGER (long), INTENT (IN) :: j ! DESCRIPTION: ! This is made available for users wishing to redefine the ! underlying metric of the space involved by redefining the ! inner product used for the updating of the quasi-Newton ! matrices. Specifically, this routine computes the general ! inner product of the vector x and the jth unit coordinate ! vector [0,...,0,1,0,...,0], where the 1 is in the jth ! position. In general, this may *not* just equal x[j]. ! EXECUTION: ! Replace code below this line if required--------------- geninnerj = x(j) ! Replace code above this line if required--------------- END FUNCTION REAL (stnd) FUNCTION norm2(v) USE reals, ONLY : zero, one USE num_constants, ONLY : machtiny, machhuge, machtol IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: v(:) ! DESCRIPTION: ! This computes the 2-norm (i.e. the Euclidean norm) of the vector v ! of length n, with due regard to avoiding overflow and underflow. ! The routine is based on snrm2 from the blas (in linpack), but this ! version is written in Fortran 90. It is machine independent. ! The machine constants MachTiny (the smallest magnitude), MachHuge(the ! largest magnitude), and MachTol (epsilon) are used to calculate the ! constants cutlo and cuthi. Three different cases must be considered ! when calculating the norm: ! (1) All components of v are below cutlo. ! To avoid underflow, each component is divided by sqrt(min)/n ! and then the regular Euclidean norm of this modified vector ! is calculated. This result is then multiplied by ! sqrt(min)/n in order to get the correct value for the norm. ! (2) One or more components are greater than cuthi. ! To avoid overflow, the same method as in case (1) is used ! with a scaling factor of sqrt(max)*n . ! (3) All components are less than cuthi, with at least one ! component greater than cutlo. ! The regular formula for the Euclidean norm is used. ! PARAMETERS: INTEGER (short), PARAMETER :: null = 0, small = 1, normal = 2, & large = 3 ! LOCAL DECLARATIONS: INTEGER (long) :: i, n INTEGER (short) :: range REAL (stnd) :: cutlo, cuthi, summ, xmax ! EXECUTION: n = SIZE(v) IF (n<=0) THEN norm2 = zero RETURN END IF cutlo = SQRT(machtiny/machtol) cuthi = SQRT(machhuge)/n summ = zero range = null ! Evaluate the norm by accumulating a scaled sum of squares and ! adjusting the scaling as numbers of increasingly large magnitude ! are found. DO i = 1, n SELECT CASE (range) CASE (normal) IF (ABS(v(i))=cuthi) THEN range = large xmax = ABS(v(i)) summ = one + (summ/v(i))/v(i) ELSE range = normal summ = (summ*xmax)*xmax + v(i)**2 END IF CASE (large) IF (ABS(v(i))<=xmax) THEN summ = summ + (v(i)/xmax)**2 ELSE summ = one + summ*(xmax/v(i))**2 xmax = ABS(v(i)) END IF CASE (null) IF (ABS(v(i))==zero) THEN ! just fall through... ELSE IF (ABS(v(i))<=cutlo) THEN range = small xmax = ABS(v(i)) summ = one ELSE IF (ABS(v(i))>=cuthi) THEN range = large xmax = ABS(v(i)) summ = one ELSE range = normal summ = v(i)**2 END IF END SELECT END DO SELECT CASE (range) CASE (normal,null) norm2 = SQRT(summ) CASE DEFAULT norm2 = xmax*SQRT(summ) END SELECT ! EXIT: RETURN ! FORMATS: none are defined. END FUNCTION REAL (stnd) FUNCTION gennorm(x) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:) ! DESCRIPTION: ! This is made available for users wishing to redefine the ! underlying metric of the space involved by redefining the ! inner product used for the updating of the quasi-Newton ! matrices. It is used in conjunction with GenInner and indeed ! this routine computes the norm of x associated with the ! redefined inner product. ! EXECUTION: ! GenNorm = sqrt ( x .GIP. x ) ! Activate this line if desired. ! The above line *defines* the generalized norm; the user may wish ! however to evaluate it in some other fashion, possibly to prevent ! unwanted overflows, as in the routine Norm2. gennorm = norm2(x) ! Remove this line if inner product changed. END FUNCTION REAL (stnd) FUNCTION inner(x,y) IMPLICIT NONE !-------------! ! ARGUMENTS: REAL (stnd), INTENT (IN) :: x(:) REAL (stnd), INTENT (IN) :: y(SIZE(x)) ! DESCRIPTION: ! This is provided so that an operator may be defined for ! computation of the standard inner product. ! EXECUTION: inner = DOT_PRODUCT(x,y) END FUNCTION END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> supp_codes.f90 MODULE supp_codes USE precisions, ONLY : short ! See Module Min_Codes for an explanation of how these codes work. IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: str_evaluation, str_derivative, str_returnstate !================= ! For Evaluate_f | !================= !Str_Evaluation ! Function evaluation ! just evaluate f at x ! evaluate both f and g at x ! just evaluate g at x INTEGER (short), PARAMETER, PUBLIC :: justf = 0, both = 1, justg = 2, & noop = 3 ! Evaluate_f just calls the ! user routine and then returns. !Str_Derivative ! Derivative calculation modes ! Compute from analytic formulae ! Compute both and compare ! Like CompareTest, just iteration 1 INTEGER (short), PARAMETER, PUBLIC :: analytic = 1, comparetest = 2, & firsttest = 3, differences = 4 ! Compute using finite differences !Str_ReturnState ! Return from Evaluation ! f and/or g succesfully computed. ! User requested abort. ! Function count limit exceeded. ! f could not be computed. ! g could not be computed. INTEGER (short), PARAMETER, PUBLIC :: ok = 0, abort = -1, limit = -2, & nof = -3, nog = -4, noforg = -5 ! Neither f nor g could be computed. !================ ! For Test_Done | !================ ! type of norm ! use absolute sum norm ! use sqrt(sum of squares) norm ! use max absolute term norm INTEGER (short), PARAMETER, PUBLIC :: l1 = 1, l2 = 2, linf = 3, g2 = 4 ! use sqrt(x .GIP. x) norm CONTAINS FUNCTION str_evaluation(i) ! ARGUMENTS: INTEGER (short) :: i ! LOCAL DECLARATIONS: CHARACTER (5), PARAMETER :: str(justf:noop) = (/ 'justF', 'both ', & 'justG', 'noOp '/) CHARACTER (LEN_TRIM(str(i))) :: str_evaluation ! DESCRIPTION: This function returns the specified value for the evaluation ! to be done. ! EXECUTION: str_evaluation = str(i) RETURN END FUNCTION FUNCTION str_derivative(i) ! ARGUMENTS: INTEGER (short) :: i ! LOCAL DECLARATIONS: CHARACTER (11), PARAMETER :: str(analytic:differences) = (/ & 'Analytic ', 'CompareTest', 'FirstTest ', 'Differences'/) CHARACTER (LEN_TRIM(str(i))) :: str_derivative ! DESCRIPTION: This function returns the specified value for the derivative ! method. ! EXECUTION: str_derivative = str(i) RETURN END FUNCTION FUNCTION str_returnstate(i) ! ARGUMENTS: INTEGER (short) :: i ! LOCAL DECLARATIONS: CHARACTER (6), PARAMETER :: str(noforg:ok) = (/ 'NoForG', 'NoG ', & 'NoF ', 'Limit ', 'Abort ', 'OK '/) CHARACTER (LEN_TRIM(str(i))) :: str_returnstate ! DESCRIPTION: This function returns the specified value for the return ! condition from the function evaluation. ! EXECUTION: str_returnstate = str(i) RETURN END FUNCTION END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> supp_states.f90 MODULE supp_states USE precisions, ONLY : stnd, long, short IMPLICIT NONE PUBLIC !------ -------------! !================= ! For Evaluate_f | !================= TYPE :: evalstate INTEGER (short) :: scalef, derivatives, expense, evallimit, & evaltraceunit LOGICAL :: tracef, traceg, tracedervtest END TYPE TYPE :: evalcounts INTEGER (long) :: fevals, gevals REAL (stnd) :: time END TYPE TYPE :: evalerrors REAL (stnd) :: worst, average, accumsum INTEGER (long) :: gradcnt, index, accumnumber END TYPE !==================== ! For Print_Iterate | !==================== TYPE :: printstate INTEGER (long) :: frequency INTEGER (short) :: printunit LOGICAL :: printx, printgrad END TYPE TYPE :: lastprint INTEGER (short) :: unit INTEGER (long) :: iter LOGICAL :: prx LOGICAL :: prg END TYPE TYPE :: printvalues REAL (stnd) :: time INTEGER (long) :: nextpoint TYPE (lastprint) :: last END TYPE !================ ! For Test_Done | !================ TYPE :: termstate INTEGER (short) :: thenorm LOGICAL :: traceterm INTEGER (short) :: termtraceunit LOGICAL :: usegrad, usestep, useshanno, usefunc REAL (stnd) :: fatx0, normgatx0 END TYPE TYPE :: termvalues REAL (stnd) :: normgsq, normxsq, diffsq END TYPE END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> supp_defs.f90 MODULE supp_defs ! Default values for support routines USE reals, ONLY : zero, one USE true_false, ONLY : false, true USE supp_codes, ONLY : analytic, l2 USE supp_states, ONLY : evalstate, lastprint, printvalues, evalerrors, & evalcounts, termstate, printstate USE precisions, ONLY : short, stnd, long IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: defaultevalstate, evalstate, defaultprintstate, printstate, & defaulttermstate, termstate, initevalcts, evalcounts, initevalers, & evalerrors, initprvals, printvalues !================= ! For Evaluate_f | !================= ! Default values INTEGER (short), PARAMETER, PRIVATE :: scalef = 0, & derivatives = analytic, expense = 1, evaltraceunit = 6 INTEGER (long), PARAMETER, PRIVATE :: evallimit = 300 LOGICAL, PARAMETER, PRIVATE :: tracef = false, traceg = false, & tracedervtest = false TYPE (evalstate) :: defaultevalstate = evalstate(scalef,derivatives, & expense,evallimit,evaltraceunit,tracef,traceg,tracedervtest) ! Initial values INTEGER (long), PARAMETER, PRIVATE :: fevals = 0, gevals = 0 REAL (stnd), PARAMETER, PRIVATE :: etime = zero TYPE (evalcounts) :: initevalcts = evalcounts(fevals,gevals,etime) ! Initial values for derivative est. REAL (stnd), PARAMETER, PRIVATE :: worst = zero, average = zero, & accumsum = zero INTEGER (long), PARAMETER, PRIVATE :: iteration = 0, index = 0, & accumnumber = 0 TYPE (evalerrors) :: initevalers = evalerrors(worst,average,accumsum, & iteration,index,accumnumber) !==================== ! For Print_Iterate | !==================== ! Default Values INTEGER (long), PARAMETER, PRIVATE :: frequency = 10000 ! Default Values INTEGER (short), PARAMETER, PRIVATE :: printunit = 6 LOGICAL, PARAMETER, PRIVATE :: printx = false, printgrad = false TYPE (printstate) :: defaultprintstate = printstate(frequency,printunit, & printx,printgrad) ! Initial Values REAL (stnd), PARAMETER, PRIVATE :: time = zero INTEGER (long), PARAMETER, PRIVATE :: nextpoint = 0 TYPE (printvalues) :: initprvals = printvalues(time,nextpoint, & lastprint(printunit,-1,false,false)) !================ ! For Test_Done | !================ ! Default Values INTEGER (short), PARAMETER, PRIVATE :: thenorm = l2, termtraceunit = 6 LOGICAL, PARAMETER, PRIVATE :: traceterm = false, usegrad = true, & usestep = false, useshanno = true, usefunc = false REAL (stnd), PARAMETER, PRIVATE :: fatx0 = one, normgatx0 = one TYPE (termstate) :: defaulttermstate = termstate(thenorm,traceterm, & termtraceunit,usegrad,usestep,useshanno,usefunc,fatx0,normgatx0) END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> fscale.f90 MODULE fscale ! 'USE'ed in modules Support and EvaluateF USE reals, ONLY : one, three, two USE supp_codes, ONLY : justf, justg, both USE true_false, ONLY : true USE precisions, ONLY : stnd, short IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: f_scale CONTAINS SUBROUTINE f_scale(ft,scale,fscaled,gscale,which) ! ARGUMENTS: REAL (stnd), INTENT (IN) :: ft INTEGER (short), INTENT (IN) :: scale REAL (stnd), INTENT (OUT) :: fscaled, gscale INTEGER (short), INTENT (IN), OPTIONAL :: which ! DESCRIPTION: ! This subroutine applies one of several scalings (linear or nonlinear) ! to a function, thereby modifying the function and/or gradient value. ! On Entry: ! Ft - the present function value ! Scale - the code for the type of scale desired, where the scaling ! function is one of the following: ! 1: ff(z) = 1 + z ! 2: ff(z) = z*z ! 3: ff(z) = -1 / (1 + z*z) ! 4: ff(z) = sqrt(1 + z*z) ! 5: ff(z) = z*z*z ! which - See Evaluate_F. If it is omitted, the value both is assumed. ! On Exit: ! FScaled - the scaled function value. ! GScale - gradient scaling factor. ! LOCAL DECLARATIONS: LOGICAL :: dof, dog ! EXECUTION: IF ( .NOT. PRESENT(which)) THEN dof = true dog = true ELSE dof = which == justf .OR. which == both dog = which == justg .OR. which == both END IF SELECT CASE (INT(scale)) CASE (1) ! ff(z) = 1 + f(z) IF (dof) fscaled = ft + one IF (dog) gscale = one CASE (2) ! ff(z) = z*z IF (dof) fscaled = ft*ft IF (dog) gscale = two*ft CASE (3) ! ff(z) = -1/(1+z**2) IF (dof) fscaled = -one/(one+ft**2) IF (dog) gscale = two*ft*(one/(one+ft**2))**2 CASE (4) ! ff(z) = sqrt(1+z**2) IF (dof) fscaled = SQRT(one+ft**2) IF (dog) gscale = ft/SQRT(one+ft**2) CASE (5) ! ff(z) = z*z*z IF (dof) fscaled = ft*ft*ft IF (dog) gscale = three*ft*ft CASE DEFAULT ! ff(z) = z IF (dof) fscaled = ft IF (dog) gscale = one END SELECT RETURN ! FORMATS: none. END SUBROUTINE END ! <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> <<=>> evaluatef.f90 MODULE evaluatef ! 'USE'ed in module Support USE general, ONLY : indent, writevec, cpusecs USE fscale, ONLY : f_scale USE reals, ONLY : zero, one, c100 USE num_constants, ONLY : machtol USE supp_codes, ONLY : justf, noforg, nog, nof, abort, limit, ok, & analytic, comparetest, firsttest, both, justg USE supp_defs, ONLY : initevalcts, initevalers USE true_false, ONLY : false USE precisions, ONLY : stnd, long, short USE supp_states, ONLY : evalstate, evalerrors, evalcounts IMPLICIT NONE PRIVATE !------ -------------! PUBLIC :: evaluate_f, evalstate, evalcounts, evalerrors CONTAINS ! optional RECURSIVE SUBROUTINE evaluate_f(fname,x,fx,g,job,ev,cts,first,test, & level) INTERFACE SUBROUTINE fname(x,fx,g,job) USE precisions, ONLY : stnd, short REAL (stnd), INTENT (IN) :: x(:) REAL (stnd), INTENT (OUT) :: fx, g(SIZE(x)) INTEGER (short), INTENT (INOUT) :: job END SUBROUTINE END INTERFACE ! ARGUMENTS: REAL (stnd), INTENT (INOUT) :: x(:), fx REAL (stnd), INTENT (OUT) :: g(SIZE(x)) INTEGER (short), INTENT (INOUT) :: job TYPE (evalstate), INTENT (IN) :: ev TYPE (evalcounts), INTENT (INOUT) :: cts LOGICAL, INTENT (IN), OPTIONAL :: first TYPE (evalerrors), INTENT (INOUT), OPTIONAL :: test INTEGER (short), INTENT (IN), OPTIONAL :: level ! DESCRIPTION: ! This subroutine evaluates a test function at the given point x. It returns ! the value of the function and/or the value of the gradient at x. It allows ! the application of a nonlinear scaling to the function if desired (see ! ScaleF below). It also allows the use of finite differences (see Derivatives ! below). It can also act as a noop, i.e. as a do nothing routine (see ! job below). ! On Entry: ! -------- ! First This must be present and true each time a new problem is begun. ! It causes initialization of counts described below to be done. ! FName is the name of the function to evaluate. There must be a subroutine ! provided of the form ! Subroutine FName(x,fx,g,job) ! (x, fx, g and job mean the same as in this subroutine.) ! The precise interface definition for FName appears above. ! x contains the value of the n-coordinates x(1),...,x(n) at which to ! evaluate the function. ! job = justf only evaluate the function. ! = both evaluate both. ! = justg only evaluate the gradient. ! = noop actually, if job has any value other than one of the first ! three, then just call FName with this same code for job; ! i.e. Evaluate_f should do nothing. This is intended for ! the possible convenience of the writer of FName. ! Ev This record contains a number of control fields. These are explained ! below. ! level If present, indent trace output by this number of spaces. ! Cts See below. ! Test See Finite Difference Error Estimates below. ! On Exit: ! -------- ! fx contains the function value (with the scaling applied if required). ! g contains the gradient value (with the scaling applied if required). ! job = OK The request made on the call was completed satisfactorily. ! fx and/or g are available as requested. ! Abort The minimization routine which called Evaluate_f is hereby ! requested to exit immediately to the routine which called ! it. This can be used by the routine FName to trigger ! premature termination due to circumstances of which the ! minimization routine may not be aware. ! Limit Terminate the minimization; the preset limit on the number ! of function evaluations allowed has been exceeded. ! NoF The function value could not be determined. ! NoG The gradient value could not be determined. ! NoForG Neither fx nor g could be evaluated. ! Cts See below. ! Control values in Ev: Ev contains the following fields (amongst others). ! ------------------- ! ScaleF controls the nonlinear scaling of FName. ! = 0 no effect. ! = k > 0 This routine computes and returns ff( FName(x) ), where ff ! is the k-th of the nonlinear functions of one variable ! defined in the routine F_Scale. ! Note that for certain scalings, the function value is needed to ! scale the gradient vector. So, in this case, if you call ! Evaluate_f just for a gradient value, it is necessary to supply ! the function value in fx as well, in order to do the scaling. ! Derivatives This specifies the method by which derivatives are to be ! computed, when requested. The choice is between: ! Analytic Use analytic formulae which must be coded and available in the ! user routine FName. ! Differences Use finite difference approximations. In this case, the user ! routine FName may ignore calls with job /= justf, and need only ! be able to compute function values. Further comments appear in ! the discussion of finite difference computations (below). ! CompareTest In this case both analytic and finite differences are ! computed. They are then compared and a record is kept to see to ! what extent they disagree. A record of the level of agreement is ! available through the optional argument Test. A more complete ! description is also given below. ! FirstTest This case is precisely the same as for CompareTest, with ! the sole exception that the testing only takes place on the ! first call to Evaluate_f. ! EvalLimit The maximum value allowed for the count Cts%FEvals. ! = 0 On entry specifies no maximum, i.e. any number of function ! evaluations may be done. This is inadvisable. ! = k > 0 Specifies the maximum number of times that FName may be ! called. If the function evaluation count in FEvals is greater ! than or equal to EvalLimit on entry to Evaluate_f, then the ! function is not evaluated and the return code job is set as ! above. Note that the count in FEvals does NOT include function ! evaluations used for computing finite difference gradients. ! Expense This may be used to artificially increase the expense of a ! function evaluation. If Expense is >1, then the function is ! simply evaluated Expense times. Of course, Expense-1 of these ! are totally redundant. This can be used for testing purposes. ! It is possible to trace parts of the execution of Evaluate_f. ! TraceF if true, the function value is to be printed ! TraceG if true, the gradient value is to be printed ! TraceDervTest if true, the function values are to be printed during testing ! mode ! EvalTraceUnit the unit number for output of computed values ! Note that an error message is printed when the maximum number of function ! evaluations is exceeded, provided either TraceF or TraceG is true. ! Counts available on exit. ! ------------------------ ! The following counts are available in the record Cts on exit. As noted ! above, Evaluate_f must be called once with First = true to ensure that these ! are correctly initialized. ! FEvals the number of calls to evaluate the function, i.e. calls with job ! = justf or both. Note that calls to evaluate the function for the ! purpose of doing finite differences are not included in this count. ! GEvals the number of calls to evaluate the gradient, i.e. calls with ! job = justg or both. ! time the amount of cpu time spent in Evaluate_f. The time used in ! the final scaling is included in the timing. Timing commences on ! entry to Evaluate_f, and ends just before return from Evaluate_f. ! Finite Difference Error Estimates ! --------------------------------- ! If the derivative mode is CompareTest, then the accuracy of finite ! difference/analytic derivative calculations is monitored. On return, if ! Test is present, one may find in the fields of Test: ! Worst The worst error which occurred. ! Average Estimate of the average number of figures of agreement. ! GradCnt Count of gradient evaluations where the worst error occurred. ! Index Component of the gradient where the worst error occurred. ! To be specific, when in test mode, each component of the analytic derivative ! is computed, and these are returned in g as the gradient. As well, for each ! component, a finite difference approximation is computed (as described ! below) and the relative difference between that and the analytic component ! is determined. This quantity is monitored, and the largest such value is ! recorded. In addition, Index records in which component of that gradient the ! error occurred, and GradCnt tells which gradient evaluation was in ! progress when the error occurred; i.e. GradCnt just records the current ! value of Cts%GEvals. Note that Index and GradCnt only refer to the point ! at which the largest error occurred. ! If the function and gradient evaluations are correct, one would normally ! expect the relative error to be of the order of 10**-(t/2), where t is the ! number of figures of relative accuracy of the machine in use. However, as ! the minimum is approached and the gradient components generally become very ! small, this relative accuracy may be much worse than expected. Therefore we ! also maintain an estimate of the average agreement. Here, for each component ! of each gradient computation, we compute the base 10 log of the relative ! accuracy; this is roughly the number of significant figures of agreement ! between the two values. This quantity is monitored and Average is returned ! as the average value of the number of significant figures of agreement. ! When function and gradient computations are correct, Worst will generally be ! at least as small as 10**(-t/2), although it can be more like 10**(-t/4). ! Gross blunders will usually give the worst error a value very near to 1, but ! not always. If all is well, Average will usually be about t/2; blunders ! will often result in Average being near 0 or 1. ! Finite difference computations ! For first derivatives, simple forward differences are used. ! To estimate the i-th component of the gradient of f, we compute ! ( f(x + h*e[i]) - f(x) ) / h, ! where h = eps * abs(x[i]). When x[i] = 0, we just choose h = eps. ! Here eps is sqrt(MachTol), where MachTol is the machine accuracy. This is ! used when Derivatives = Differences, CompareTest or FirstTest. ! When Derivatives = CompareTest, more information is required; thus we also ! compute f(x + sqrt(h)*e[i]). This means that when in test mode, twice as ! many function evaluations are needed. This is required to eliminate ! scaling effects in the estimate of figures of agreement. ! SUBROUTINES: ! FName ...the user routine. ! F_Scale ...performs scaling. ! PARAMETERS: CHARACTER (*), PARAMETER :: id = 'Eval' ! LOCAL DECLARATIONS: LOGICAL :: firsttime, docompare, dof, dog INTEGER (short) :: rememberbad, dervs, job_return, calls INTEGER (long) :: count, kk, zdcnt, zindex, zgcnt REAL (stnd) :: ft, fv, tt, scale, rh, zerr, zserr, fval, fval2, h, & terr ! EXECUTION: IF (ev%tracef .OR. ev%traceg) THEN CALL indent(id,ev%evaltraceunit,level) WRITE (ev%evaltraceunit,90000) job, justf, justg, both END IF IF (job/=justf .AND. job/=justg .AND. job/=both) THEN ! noop call. CALL fname(x,fx,g,job) GO TO 10 END IF dervs = ev%derivatives IF (PRESENT(first)) THEN firsttime = first ELSE firsttime = false END IF IF (firsttime) THEN cts = initevalcts IF (PRESENT(test)) test = initevalers IF (ev%derivatives==firsttest) dervs = comparetest ELSE IF (ev%derivatives==firsttest) dervs = analytic IF (ev%evallimit>0 .AND. cts%fevals>=ev%evallimit) THEN GO TO 20 END IF END IF dof = job == justf .OR. job == both dog = job == justg .OR. job == both cts%time = cts%time - cpusecs() IF (dof) cts%fevals = cts%fevals + 1 IF (dog) cts%gevals = cts%gevals + 1 !-----First compute required function and/or gradient values. docompare = PRESENT(test) .AND. dervs == comparetest IF (docompare) THEN zerr = test%worst zserr = test%accumsum zdcnt = test%accumnumber zindex = test%index zgcnt = test%gradcnt END IF rememberbad = ok EXPENSE: DO kk = 1, ev%expense ! Repeat if required to simulate expensive call. IF (dervs==analytic .OR. .NOT. dog) THEN calls = 0 ELSE calls = SIZE(x) ! extra calls to FName needed. END IF job_return = job ! First compute f(x) and/or g(x). CALL fname(x,fval,g,job_return) SELECT CASE (job_return) CASE (limit,abort,nof,nog,noforg) rememberbad = job_return CYCLE END SELECT IF (job==justg) THEN IF (scale/=0) ft = fx ELSE ft = fval END IF NCALLS: DO count = 1, calls !Do extra calls, if required: function values only. tt = x(count) ! save x h = SQRT(machtol) IF (tt/=zero) THEN ! finite difference step h = h*ABS(tt) END IF x(count) = tt + h ! Compute f( x + h * e[count] ) job_return = justf CALL fname(x,fval,g,job_return) x(count) = tt ! restore x SELECT CASE (job_return) CASE (limit,abort,nof,nog,noforg) rememberbad = job_return CYCLE EXPENSE END SELECT DOCOMPAR: IF (docompare) THEN IF (ev%tracedervtest) THEN CALL indent(id,ev%evaltraceunit,level) WRITE (ev%evaltraceunit,90030) g(count), count, (fval-ft)/h END IF ! Estimate error, and leave computed analytic gradients in g. ! Use f at x + a * e[count], for a = h and sqrt(h). rh = SQRT(h) x(count) = tt + rh job_return = justf CALL fname(x,fval2,g,job_return) x(count) = tt SELECT CASE (job_return) CASE (limit,abort,nof,nog,noforg) rememberbad = job_return CYCLE EXPENSE END SELECT IF (ABS(fval2-ft)>c100*machtol*ABS(ft)) THEN terr = (fval-ft-h*g(count))/(fval2-ft-rh*g(count)) IF (tt>one) terr = terr/tt terr = MAX(MIN(one,ABS(terr)),machtol) ! in [MachTol,1] zserr = zserr - LOG10(terr) ! Number of figures agreement. zdcnt = zdcnt + 1 IF (ev%tracedervtest) THEN CALL indent(id,ev%evaltraceunit,level) WRITE (ev%evaltraceunit,90020) terr, -LOG10(terr) END IF IF (terr>ABS(zerr)) THEN zindex = count zgcnt = cts%gevals zerr = SIGN(terr,zerr) END IF ELSE zerr = -ABS(zerr) ! Flag case with excessive cancellation. IF (ev%tracedervtest) THEN CALL indent(id,ev%evaltraceunit,level) WRITE (ev%evaltraceunit,90010) END IF END IF ELSE ! DOCOMPAR Estimate gradients with forward finite difference g(count) = (fval-ft)/h END IF DOCOMPAR END DO NCALLS IF (ev%scalef/=0) THEN ! Do scaling: define Fv and scale. CALL f_scale(ft,ev%scalef,fv,scale,job) ELSE fv = ft scale = one END IF IF (dof) THEN ! Revise function and gradient as necessary. fx = fv END IF IF (dog .AND. scale/=one) THEN g = scale*g END IF END DO EXPENSE IF (docompare) THEN test%worst = zerr test%index = zindex test%gradcnt = zgcnt test%accumnumber = zdcnt test%accumsum = zserr test%average = zserr/zdcnt