SUBROUTINE UCTPG(PNO, M,N, X,TAU0,DX0, XMIN,FXMIN) ************************************************************************ * Get parameters for unconstrained minimization test problem. * See H.B. Nielsen: "UCTP - Test Problems for Unconstrained * Optimization", Report IMM-REP-2000-18 * Version 00.11.27 ************************************************************************ IMPLICIT NONE C Parameters INTEGER PNO,M,N DOUBLE PRECISION X(N),TAU0,DX0,XMIN(N),FXMIN C Common variables INTEGER PBNO, MM DOUBLE PRECISION PBDATA(2000) COMMON /PROB/ PBNO, MM, PBDATA C Local variables INTEGER I,J,MN,MNA, MEYER(16) DOUBLE PRECISION BARD(15),KANDO(22),OSB(33),EFIT(45),XWAT(27), " XBD(12),XCHQ(34),XOSB(10),XEFIT(8), " DTAU0(22),DDX0(22), A,BB,T,Y INTRINSIC ABS, DBLE, MIN, EXP, COS, SIN C BLAS functions DOUBLE PRECISION DDOT, DNRM2 C If BLAS is available, then remove the * in col. 1 of the next line * EXTERNAL DCOPY, DDOT, DGEMV, DNRM2, DSYRK C DATA DTAU0/3*1D-8, 2*1D0, 1D-8, 1D0, 1D-8, 1D0, 1D-6, & 2*1D-8, 1D0, 1D-3, 2*1D0, 1D-8, 2*1D-3, 2*1D0, 1D-3/ DATA DDX0 /3*1D1, 6*1D0, 1D2, & 2*1D0, 5D-2, 5D-1, 2*1D0, 1D-1, 3*1D0, 1D2, 1D0/ DATA BARD /0.14D0, 0.18D0, 0.22D0, 0.25D0, 0.29D0, & 0.32D0, 0.35D0, 0.39D0, 0.37D0, 0.58D0, & 0.73D0, 0.96D0, 1.34D0, 2.10D0, 4.39D0/ DATA KANDO /0.1957D0, 0.1947D0, 0.1735D0, 0.1600D0, 0.0844D0, & 0.0627D0, 0.0456D0, 0.0342D0, 0.0323D0, 0.0235D0, 0.0246D0, & 4.0000D0, 2.0000D0, 1.0000D0, 0.5000D0, 0.2500D0, & 0.1670D0, 0.1250D0, 0.1000D0, 0.0833D0, 0.0714D0, 0.0625D0/ DATA MEYER /34780, 28610, 23650, 19630, 16370, 13720, 11540, & 9744, 8261, 7030, 6005, 5147, 4427, 3820, 3307, 2872/ DATA OSB /8.44D-1, 9.08D-1, 9.32D-1, 9.36D-1, 9.25D-1, 9.08D-1, & 8.81D-1, 8.50D-1, 8.18D-1, 7.84D-1, 7.51D-1, 7.18D-1, 6.85D-1, & 6.58D-1, 6.28D-1, 6.03D-1, 5.80D-1, 5.58D-1, 5.38D-1, 5.22D-1, & 5.06D-1, 4.90D-1, 4.78D-1, 4.67D-1, 4.57D-1, 4.48D-1, 4.38D-1, & 4.31D-1, 4.24D-1, 4.20D-1, 4.14D-1, 4.11D-1, 4.06D-1/ DATA EFIT / 0.905420D-1, 1.245690D-1, 1.793670D-1, 1.956540D-1, & 2.697070D-1, 2.860270D-1, 2.898920D-1, 3.174750D-1, 3.081910D-1, & 3.369950D-1, 3.483710D-1, 3.213370D-1, 2.994230D-1, 3.389720D-1, & 3.047630D-1, 2.889030D-1, 3.008200D-1, 3.039740D-1, 2.839870D-1, & 2.620780D-1, 2.815930D-1, 2.675310D-1, 2.189260D-1, 2.255720D-1, & 2.005940D-1, 1.973750D-1, 1.824400D-1, 1.838920D-1, 1.522850D-1, & 1.740280D-1, 1.508740D-1, 1.262200D-1, 1.262660D-1, 1.063840D-1, & 1.189230D-1, 0.918680D-1, 1.289260D-1, 1.192730D-1, 1.159970D-1, & 1.058310D-1, 0.752610D-1, 0.683870D-1, 0.908230D-1, 0.852050D-1, & 0.672030D-1/ DATA XWAT /-1.57251D-2, 1.01243D0,-2.32992D-1, 1.26043D0, & -1.51373D0, 9.92996D-1, -1.53070D-5, 9.99790D-1, & 1.47640D-2, 1.46342D-1,1.00082D0, -2.61773D0, 4.10440D0, & -3.14361D0, 1.05263D0, -6.65800D-9, 1.00000D0, & -5.62881D-4, 3.47798D-1,-1.56525D-1,1.05178D0, -3.24414D0, & 7.28248D0, -1.02647D1, 9.06875D0,-4.53913D0, 1.01161D0/ DATA XBD / 0.77781D0, 1.87187D0, 1.15301D0, -0.67095D0, & -0.18950D0, 3.45424D0, 1.32570D0, -1.336679D0, & -11.5944D0, 13.2036D0, -0.40339D0, 0.236779D0/ DATA XCHQ / 4.31528D-2, 1.93091D-1, 2.66329D-1, 5.00000D-1, & 5.00000D-1, 7.33671D-1, 8.06909D-1, 9.56847D-1, & 7.97431D-2, 1.78355D-1, 3.10120D-1, 4.41205D-1, & 5.58795D-1, 6.89880D-1, 8.21645D-1, 9.20257D-1, & 4.42053D-2, 1.99491D-1, 2.35619D-1, 4.16047D-1, 5.00000D-1, & 5.83953D-1, 7.64381D-1, 8.00509D-1, 9.55795D-1, & 1.01891D-1, 1.89483D-1, 2.94482D-1, 3.97099D-1, 5.00000D-1, & 6.02901D-1, 7.05518D-1, 8.10517D-1, 8.98109D-1 / DATA XOSB / .37541D0, 1.93585D0, -1.46469D0, .1287D-1, & .2212D-1, .5D0, 1.5D0, -1D0, 1D-2, 2D-2 / DATA XEFIT / -4D0, -5D0, 4D0, -4D0, -1D0, -2D0, 1D0, -1D0/ C PBNO = PNO MN = M * N TAU0 = DTAU0(PNO) DX0 = DDX0(PNO) IF (PNO .EQ. 1) THEN C ... Linear function. Full rank CALL DCONST(-2D0/DBLE(M), MN, PBDATA) CALL DCONST(-1D0, M, PBDATA(MN+1)) I = 1 DO 10 J = 1,N PBDATA(I) = PBDATA(I) + 1D0 I = I + M + 1 10 CONTINUE C ... Solution CALL DCONST(-1D0, N, XMIN) FXMIN = DBLE(M - N) / 2D0 C ... Starting point CALL DCONST(1D0, N, X) ELSEIF (PNO .EQ. 2) THEN C ... Linear function. Rank 1 MN = 0 DO 30 J = 1, N DO 20 I = 1, M 20 PBDATA(MN+I) = DBLE(I*J) 30 MN = MN + M CALL DCONST(-1D0, M, PBDATA(MN+1)) C ... Solution CALL DCONST(0D0, N, XMIN) XMIN(1) = 3D0 / DBLE(2*M + 1) FXMIN = DBLE(M*(M - 1)) / DBLE(8*M + 4) C ... Starting point CALL DCONST(1D0, N, X) ELSEIF (PNO .EQ. 3) THEN C ... Linear function. Rank 1. Zero columns and rows CALL DCONST(0D0, M*N, PBDATA) MN = M+1 DO 50 J = 2, N-2 DO 40 I = 1, M-2 40 PBDATA(MN+I) = DBLE(I*J) 50 MN = MN + M CALL DCONST(-1D0, M, PBDATA(M*N+1)) C ... Solution CALL DCONST(0D0, N, XMIN) XMIN(2) = 1.5D0 / DBLE(2*M -3) FXMIN = DBLE(M*(M + 3) - 6) / DBLE(8*M - 12) C ... Starting point CALL DCONST(1D0, N, X) ELSEIF (PNO .EQ. 4) THEN C ... Rosenbrock C ... Solution CALL DCONST(1D0, N, XMIN) FXMIN = 0D0 C ... Starting point X(1) = -1.2D0 X(2) = 1D0 ELSEIF (PNO .EQ. 5) THEN C ... Helical Valley C ... Solution CALL DCONST(0D0, N, XMIN) XMIN(1) = 1D0 FXMIN = 0D0 C ... Starting point CALL DCONST(0D0, N, X) X(1) = -1D0 ELSEIF (PNO .EQ. 6) THEN C ... Powell singular C ... Solution CALL DCONST(0D0, 4, XMIN) FXMIN = 0D0 C ... Starting point X(1) = 3D0 X(2) = -1D0 X(3) = 0D0 X(4) = 1D0 ELSEIF (PNO .EQ. 7) THEN C ... Freudenstein and Roth C ... Solution XMIN(1) = 11.4128D0 XMIN(2) = -0.896805D0 FXMIN = 24.4921D0 C ... Starting point X(1) = .5D0 X(2) = -2D0 ELSEIF (PNO .EQ. 8) THEN C ... Bard M = 15 CALL DCOPY(M, BARD,1, PBDATA,1) DO 60 I = 1, M PBDATA(15+I) = DBLE(I) PBDATA(30+I) = DBLE(16-I) PBDATA(45+I) = DBLE(MIN(I, 16-I)) 60 CONTINUE C ... Solution XMIN(1) = 8.2411D-2 XMIN(2) = 1.133036D0 XMIN(3) = 2.343695D0 FXMIN = 4.10744D-3 C ... Starting point CALL DCONST(1D0, N,X) ELSEIF (PNO .EQ. 9) THEN C ... Kowalik and Osborne M = 11 CALL DCOPY(2*M, KANDO,1, PBDATA,1) C ... Solution XMIN(1) = .192807D0 XMIN(2) = .191282D0 XMIN(3) = .123057D0 XMIN(4) = .136062D0 FXMIN = 1.53753D-4 C ... Starting point X(1) = .25D0 X(2) = .39D0 X(3) = .415D0 X(4) = .39D0 ELSEIF (PNO .EQ. 10) THEN C ... Meyer DO 70 I = 1, M PBDATA(I) = DBLE(45 + 5*I) PBDATA(M+I) = DBLE(MEYER(I)) 70 CONTINUE C ... Solution XMIN(1) = 5.60964D-3 XMIN(2) = 6.18135D3 XMIN(3) = 3.45224D2 FXMIN = 43.9729D0 C ... Starting point X(1) = 2D-2 X(2) = 4D3 X(3) = 2.5D2 ELSEIF (PNO .EQ. 11) THEN C ... Watson MNA = M*N CALL DCONST(1D0, M, PBDATA) CALL DCONST(1D0, M, PBDATA(MNA+1)) DO 80 I = 1, M 80 PBDATA(M+I) = DBLE(I) / 29D0 MN = M DO 100 J = 3, N MN = MN + M MNA = MNA + M DO 90 I = 1, M BB = PBDATA(MN-M+I) PBDATA(MNA+I) = DBLE(J-1) * BB PBDATA(MN+I) = PBDATA(M+I) * BB 90 CONTINUE 100 CONTINUE C ... Solution IF (N .EQ. 6) THEN CALL DCOPY(6, XWAT,1, XMIN,1) FXMIN = 1.143835D-3 ELSEIF (N .EQ. 9) THEN CALL DCOPY(9, XWAT(7),1, XMIN,1) FXMIN = 6.998801D-7 ELSEIF (N .EQ. 12) THEN CALL DCOPY(12, XWAT(16),1, XMIN,1) FXMIN = 2.361196D-10 ENDIF C ... Starting point CALL DCONST(0D0, N, X) ELSEIF (PNO .EQ. 12) THEN C ... Box 3D DO 110 I = 1, M PBDATA(I) = DBLE(I) / 10D0 PBDATA(M+I) = EXP(-PBDATA(I)) - EXP(-DBLE(I)) 110 CONTINUE C ... Solution XMIN(1) = 1D0 XMIN(2) = 1D1 XMIN(3) = 1D0 FXMIN = 0D0 C ... Starting point X(1) = 0D0 X(2) = 1D1 X(3) = 2D1 ELSEIF (PNO .EQ. 13) THEN C ... Jennrich and Sampson DO 120 I = 1, M PBDATA(I) = DBLE(I) PBDATA(M+I) = DBLE(2 + 2*I) 120 CONTINUE C ... Solution IF (M .EQ. 5) THEN CALL DCONST(.378468D0, 2, XMIN) FXMIN = 4.8879031D0 ELSEIF (M .EQ. 10) THEN CALL DCONST(.257825D0, 2, XMIN) FXMIN = 6.21811D1 ELSEIF (M .EQ. 20) THEN CALL DCONST(.165191D0, 2, XMIN) FXMIN = 7.2474D2 ENDIF C ... Starting point X(1) = .3D0 X(2) = .4D0 ELSEIF (PNO .EQ. 14) THEN C ... Brown and Dennis J = 1 DO 130 I = 1, M PBDATA(J) = DBLE(I)/5D0 PBDATA(J+1) = EXP(PBDATA(J)) PBDATA(J+2) = SIN(PBDATA(J)) PBDATA(J+3) = COS(PBDATA(J)) J = J + 4 130 CONTINUE C ... Solution IF (M .EQ. 5) THEN CALL DCOPY(4, XBD,1, XMIN,1) FXMIN = 9.08309D-5 ELSEIF (M .EQ. 10) THEN CALL DCOPY(4, XBD(5),1, XMIN,1) FXMIN = 7.21613D-1 ELSEIF (M .EQ. 20) THEN CALL DCOPY(4, XBD(9),1, XMIN,1) FXMIN = 4.29111D4 ENDIF C ... Starting point X(1) = 25D0 X(2) = 5D0 X(3) = -5D0 X(4) = -1D0 ELSEIF (PNO .EQ. 15) THEN C ... Chebyquad CALL DCONST(0D0, M, PBDATA) DO 140 I = 2, M, 2 140 PBDATA(I) = 1D0 / DBLE(I**2 - 1) C ... Solution IF ((N .EQ. 8) .AND. (M .EQ. 8)) THEN CALL DCOPY(8, XCHQ,1, XMIN,1) FXMIN = 1.75844D-3 ELSEIF ((N .EQ. 8) .AND. (M .EQ. 16)) THEN CALL DCOPY(8, XCHQ(9),1, XMIN,1) FXMIN = 2.94780D-2 ELSEIF ((N .EQ. 9) .AND. (M .EQ. 9)) THEN CALL DCOPY(9, XCHQ(17),1, XMIN,1) FXMIN = 0D0 ELSEIF ((N .EQ. 9) .AND. (M .EQ. 18)) THEN CALL DCOPY(9, XCHQ(26),1, XMIN,1) FXMIN = 3.55274D-2 ENDIF C ... Starting point DO 150 I = 1, N 150 X(I) = DBLE(I) / DBLE(N+1) DX0 = 1D0 / DBLE(N + 1) ELSEIF (PNO .EQ. 16) THEN C ... Brown almost linear CALL DCONST(DBLE(N+1), N-1, PBDATA) CALL DCONST(1D0, (N-1)*N+1, PBDATA(N)) DO 160 I = N+1, N*N, N 160 PBDATA(I) = 2D0 C ... Solution CALL DCONST(1D0, N, XMIN) FXMIN = 0D0 C ... Starting point CALL DCONST(.5D0, N, X) ELSEIF (PNO .EQ. 17) THEN C ... Osborne 1 DO 170 I = 1, 33 170 PBDATA(I) = DBLE(10*(I - 1)) CALL DCOPY(33, OSB,1, PBDATA(34),1) C ... Solution CALL DCOPY(N, XOSB,1, XMIN,1) FXMIN = 2.73245D-5 C ... Starting point CALL DCOPY(N, XOSB(6),1, X,1) ELSEIF (PNO .EQ. 18) THEN C ... Exp fit, 4 parameters DO 180 I = 1, M 180 PBDATA(I) = .02D0 * DBLE(I) CALL DCOPY(M, EFIT,1, PBDATA(M+1),1) C ... Solution CALL DCOPY(N, XEFIT,1, XMIN,1) FXMIN = 5D-3 C ... Starting point CALL DCOPY(N, XEFIT(N+1),1, X,1) ELSEIF (PNO .EQ. 19) THEN C ... Exp fit, 2 parameters DO 190 I = 1, M 190 PBDATA(I) = .02D0 * DBLE(I) CALL DCOPY(M, EFIT,1, PBDATA(M+1),1) C ... Solution CALL DCOPY(N, XEFIT,1, XMIN,1) FXMIN = 5D-3 C ... Starting point CALL DCOPY(N, XEFIT(5),1, X,1) ELSEIF (PNO .EQ. 20) THEN C ... Modified Meyer DO 200 I = 1, M PBDATA(I) = 1D-2*DBLE(45 + 5*I) PBDATA(M+I) = 1D-3*DBLE(MEYER(I)) 200 CONTINUE C ... Solution XMIN(1) = 2.48178D0 XMIN(2) = 6.18135D0 XMIN(3) = 3.45224D0 FXMIN = 4.39729D-5 C ... Starting point X(1) = 8.85D0 X(2) = 4D0 X(3) = 2.5D0 ELSEIF (PNO .EQ. 21) THEN C ... Separated Meyer DO 210 I = 1, M PBDATA(I) = DBLE(45 + 5*I) PBDATA(M+I) = DBLE(MEYER(I)) 210 CONTINUE C ... Solution XMIN(1) = 6.18135D3 XMIN(2) = 3.45224D2 FXMIN = 43.9729D0 C ... Starting point X(1) = 4D3 X(2) = 2.5D2 ELSEIF (PNO .EQ. 22) THEN C ... Exp and squares A = 1D0 DO 220 I = 1, N 220 A = A + 1D0/DBLE(I)**2 C ... Newton iteration to solve -A*exp(-y) + y = 0 Y = .75D0 230 T = Y BB = A * EXP(-Y) Y = Y - (Y - BB)/(1D0 + BB) IF (ABS(Y - T) .GT. 1D-12) GOTO 230 C ... Solution A = EXP(-Y) FXMIN = 0D0 DO 240 I = 1, N BB = DBLE(I)**2 XMIN(I) = A / BB FXMIN = FXMIN + BB*XMIN(I)**2 240 CONTINUE FXMIN = A + .5D0 * FXMIN C ... Starting point CALL DCONST(0D0, N, X) ENDIF C ... Store numer of functions in PROB MM = M RETURN END C SUBROUTINE UCTPV(N,X,DF,F) ************************************************************************ * Unconstrained minimization test problem. Scalar-gradient formulation * See H.B. Nielsen: "UCTP - Test Problems for Unconstrained * Optimization", Report IMM-REP-2000-18 * Version 00.11.27 ************************************************************************ IMPLICIT NONE C Parameters INTEGER N DOUBLE PRECISION X(N),DF(N),F C Common variables INTEGER PBNO, MM DOUBLE PRECISION PBDATA(2000) COMMON /PROB/ PBNO, MM, PBDATA C Local variables INTEGER J DOUBLE PRECISION W(2000), E,Y,XX INTRINSIC DBLE, EXP C BLAS functions DOUBLE PRECISION DNRM2 C If BLAS is available, then remove the * in col. 1 of the next line * EXTERNAL DGEMV, DNRM2 C IF (PBNO .LE. 21) THEN C ... Basically a least squares problem. Use UCTPV2 IF (MM*(N+1) .GT. 2000) THEN PRINT*, 'WORKSPACE W(2000) SHOULD BE INCREASED TO W(', & MM*(N+1),')' RETURN ENDIF CALL UCTPV2(N,MM,X,W(MM+1),W) F = .5D0 * DNRM2(MM, W, 1)**2 CALL DCONST(0D0, N, DF) CALL DGEMV('T', MM,N, 1D0,W(MM+1),MM, W,1, 0D0,DF,1) ELSEIF (PBNO .EQ. 22) THEN C ... Exp and squares Y = 0 DO 10 J = 1, N 10 Y = Y + X(J) E = EXP(-Y) F = 0 DO 20 J = 1, N XX = DBLE(J)**2 * X(J) F = F + XX * X(J) DF(J) = XX - E 20 CONTINUE F = E + .5D0 * F ENDIF RETURN END C SUBROUTINE UCTPV2(N,M,X,DF,F) ************************************************************************ * Unconstrained minimization test problem. Vector-Jacobian formulation * See H.B. Nielsen: "UCTP - Test Problems for Unconstrained * Optimization", Report IMM-REP-2000-18 * Version 00.11.27 ************************************************************************ IMPLICIT NONE C Parameters INTEGER N,M DOUBLE PRECISION X(N),DF(M,N),F(M) C Common variables INTEGER PBNO, MM DOUBLE PRECISION PBDATA(2000) COMMON /PROB/ PBNO, MM, PBDATA C Local variables INTEGER I,J,MN,MNA DOUBLE PRECISION PI4,T,NX,C,S5,S10,D3,D4,X1,X2,U,II,NN,T0,T1,T2, " D0,D1,D2,A(2,2),BY(2),BD(2,2),C1,C2,EE,EI INTRINSIC ATAN, SQRT, EXP C BLAS functions DOUBLE PRECISION DDOT, DNRM2 C If BLAS is available, then remove the * in col. 1 of the next line * EXTERNAL DCOPY, DDOT, DGEMV, DNRM2, DSYRK C IF (PBNO .LE. 3) THEN C ... Linear function CALL DCOPY(M*N, PBDATA,1, DF,1) CALL DCOPY(M, PBDATA(M*N+1),1, F,1) CALL DGEMV('N', M,N, 1D0,DF,M, X,1, 1D0,F,1) ELSEIF (PBNO .EQ. 4) THEN C ... Rosenbrock F(1) = 1D1*(X(2) - X(1)**2) F(2) = 1D0 - X(1) DF(1,1) = -2D1*X(1) DF(1,2) = 1D1 DF(2,1) = -1D0 DF(2,2) = 0D0 ELSEIF (PBNO .EQ. 5) THEN C ... Helical Valley PI4 = ATAN(1D0) T = ATAN(X(2)/X(1)) / (8D0 * PI4) IF (X(1) .LT. 0D0) T = T + .5D0 NX = DNRM2(2, X,1) F(1) = 1D1 * (X(3) - 1D1*T) F(2) = 1D1*(NX - 1D0) F(3) = X(3) C = -12.5D0 / PI4 / NX**2 DF(1,1) = -C*X(2) DF(1,2) = C*X(1) DF(1,3) = 1D1 C = 1D1 / NX DF(2,1) = C*X(1) DF(2,2) = C*X(2) DF(2,3) = 0D0 DF(3,1) = 0D0 DF(3,2) = 0D0 DF(3,3) = 1D0 ELSEIF (PBNO .EQ. 6) THEN C ... Powell singular S5 = SQRT(5D0) S10 = SQRT(1D1) D3 = X(2) - 2D0*X(3) D4 = X(1) - X(4) F(1) = x(1) + 1D1*x(2) F(2) = s5*(x(3) - x(4)) F(3) = D3**2 F(4) = S10 * D4**2 CALL DCONST(0D0, 16, DF) DF(1,1) = 1D0 DF(1,2) = 1D1 DF(2,3) = S5 DF(2,4) = -S5 DF(3,2) = 2D0*D3 DF(3,3) = -4D0*D3 DF(4,1) = 2D0 * S10 * D4 DF(4,4) = -DF(4,1) ELSEIF (PBNO .EQ. 7) THEN C ... Freudenstein and Roth X1 = X(1) X2 = X(2) F(1) = X1 - X2*(2D0 - X2*(5D0 - X2)) - 13D0 F(2) = X1 - X2*(14D0 - X2*(1D0 + X2)) - 29D0 DF(1,1) = 1D0 DF(1,2) = -2D0 + x2*(1D1 - 3D0*x2) DF(2,1) = 1D0 DF(2,2) = -14D0 + x2*(2D0 + 3D0*x2) ELSEIF (PBNO .EQ. 8) THEN C ... Bard DO 10 I = 1, M C = X(2)*PBDATA(30+I) + X(3)*PBDATA(45+I) F(I) = PBDATA(I) - (X(1) + PBDATA(15+I)/C) D3 = PBDATA(15+I) / C**2 DF(I,1) = -1D0 DF(I,2) = PBDATA(30+I) * D3 DF(I,3) = PBDATA(45+I) * D3 10 CONTINUE ELSEIF (PBNO .EQ. 9) THEN C ... Kowalik and Osborne DO 20 I = 1, M U = PBDATA(11+I) D3 = U + X(3) D4 = U*D3 + X(4) C = U / D4 F(I) = PBDATA(I) - C*X(1)*(U + X(2)) DF(I,1) = -C * (U + X(2)) DF(I,2) = -C * X(1) DF(I,4) = C * X(1) * (U + X(2)) / D4 DF(I,3) = U * DF(I,4) 20 CONTINUE ELSEIF (PBNO .EQ. 10) THEN C ... Meyer DO 30 I = 1, M C = PBDATA(I) + X(3) U = EXP(X(2)/C) F(I) = X(1) * U - PBDATA(M+I) DF(I,1) = U DF(I,2) = X(1)/C*U DF(I,3) = -DF(I,2)*X(2)/C 30 CONTINUE ELSEIF (PBNO .EQ. 11) THEN C ... Watson CALL DCONST(0D0, M*N, DF) DO 50 I = 1, 29 MNA = M * N MN = M C = DDOT(N, PBDATA(I),M, X,1) F(I) = DDOT(N-1, PBDATA(MNA+I),M, X(2),1) - C**2 - 1D0 DF(I,1) = -2D0*PBDATA(I)*C DO 40 J = 2, N DF(I,J) = PBDATA(MNA+I) - 2D0*PBDATA(MN+I)*C MN = MN + M MNA = MNA + M 40 CONTINUE 50 CONTINUE F(30) = X(1) DF(30,1) = 1D0 F(31) = X(2) - X(1)**2 - 1D0 DF(31,1) = -2D0 * X(1) DF(31,2) = 1D0 ELSEIF (PBNO .EQ. 12) THEN C ... Box 3D DO 60 I = 1, M X1 = EXP(-X(1)*PBDATA(I)) X2 = EXP(-X(2)*PBDATA(I)) F(I) = X1 - X2 - PBDATA(M+I)*X(3) DF(I,1) = -PBDATA(I)*X1 DF(I,2) = PBDATA(I)*X2 DF(I,3) = -PBDATA(M+I) 60 CONTINUE ELSEIF (PBNO .EQ. 13) THEN C ... Jennrich and Sampson DO 70 I = 1, M X1 = EXP(X(1)*PBDATA(I)) X2 = EXP(X(2)*PBDATA(I)) F(I) = PBDATA(M+I) - X1 - X2 DF(I,1) = -PBDATA(I)*X1 DF(I,2) = -PBDATA(I)*X2 70 CONTINUE ELSEIF (PBNO .EQ. 14) THEN C ... Brown and Dennis J = 1 DO 80 I = 1, M X1 = X(1) + X(2)*PBDATA(J) - PBDATA(J+1) X2 = X(3) + X(4)*PBDATA(J+2) - PBDATA(J+3) F(I) = X1**2 + X2**2 DF(I,1) = 2D0*X1 DF(I,2) = 2D0*PBDATA(J)*X1 DF(I,3) = 2D0*X2 DF(I,4) = 2D0*X2*PBDATA(J+2) J = J + 4 80 CONTINUE ELSEIF (PBNO .EQ. 15) THEN C ... Chebyquad NN = DBLE(N) CALL DCOPY(M, PBDATA,1, F,1) DO 100 J = 1, N C = 2D0 * (2D0*X(J) - 1D0) T1 = 1D0 / NN T2 = (.5D0 * C) / NN D1 = 0D0 D2 = 2D0 / NN F(1) = F(1) + T2 DF(1,J) = D2 DO 90 I = 2, M D0 = D1 D1 = D2 D2 = 4D0*T2 + C*D2 - D0 DF(I,J) = D2 T0 = T1 T1 = T2 T2 = C*T2 - T0 F(I) = F(I) + T2 90 CONTINUE 100 CONTINUE ELSEIF (PBNO .EQ. 16) THEN C ... Brown almost linear CALL DCOPY(N, PBDATA,1, F,1) CALL DGEMV('N', M-1,N, 1D0,PBDATA(N+1),N-1, X,1, -1D0,F,1) U = 1D0 DO 110 I = 1, N 110 U = U * X(I) F(N) = U - 1D0 MN = N + 1 DO 120 J = 1, N CALL DCOPY(N-1, PBDATA(MN),1, DF(1,J),1) MN = MN + N - 1 DF(N,J) = U / X(J) 120 CONTINUE ELSEIF (PBNO .EQ. 17) THEN C ... Osborne 1 DO 130 I = 1, M T = PBDATA(I) D1 = EXP(-X(4)*T) D2 = EXP(-X(5)*T) F(I) = PBDATA(M+I) - (X(1) + D1*X(2) + D2*X(3)) DF(I,1) = -1D0 DF(I,2) = -D1 DF(I,3) = -D2 DF(I,4) = T*X(2)*D1 DF(I,5) = T*X(3)*D2 130 CONTINUE ELSEIF (PBNO .EQ. 18) THEN C ... Exp fit, 4 parameters DO 140 I = 1, M T = PBDATA(I) D1 = -EXP(X(1)*T) D2 = -EXP(X(2)*T) F(I) = PBDATA(M+I) + D1*X(3) + D2*X(4) DF(I,1) = T*X(3)*D1 DF(I,2) = T*X(4)*D2 DF(I,3) = D1 DF(I,4) = D2 140 CONTINUE ELSEIF (PBNO .EQ. 19) THEN C ... Exp fit, 2 parameters. Model DO 150 I = 1, M T = PBDATA(I) DF(I,1) = EXP(X(1)*T) DF(I,2) = EXP(X(2)*T) 150 CONTINUE C ... Set up and solve normal eqsuations CALL DCOPY(M, PBDATA(M+1),1, F,1) CALL DSYRK('U','T',2,M, 1D0,DF,M, 0D0,A,2) CALL DGEMV('T', M,2, 1D0,DF,M, F,1, 0D0, BY,1) A(2,1) = A(1,2)/A(1,1) A(2,2) = A(2,2) - A(2,1)*A(1,2) BY(2) = (BY(2) - A(2,1)*BY(1))/A(2,2) BY(1) = (BY(1) - A(1,2)*BY(2))/A(1,1) CALL DGEMV('N', M,2, -1D0,DF,M, BY,1, 1D0, F,1) C ... Gradients CALL DCONST(0D0, 4, BD) C1 = BY(1) C2 = BY(2) DO 170 I = 1, M T = PBDATA(I) DO 160 J = 1, 2 BD(J,J) = BD(J,J) + T * F(I) * DF(I,J) BD(J,1) = BD(J,1) - T * C1 * DF(I,J) * DF(I,1) BD(J,2) = BD(J,2) - T * C2 * DF(I,J) * DF(I,2) 160 CONTINUE 170 CONTINUE C ... Solve for c' DO 180 J = 1, 2 BD(2,J) = (BD(2,J) - A(2,1)*BD(1,J))/A(2,2) BD(1,J) = (BD(1,J) - A(1,2)*BD(2,J))/A(1,1) 180 CONTINUE C ... Store Jacobian DO 190 I = 1, M T = PBDATA(I) X1 = DF(I,1) X2 = DF(I,2) DF(I,1) = -T * C1 * X1 - X1 * BD(1,1) - X2 * BD(2,1) DF(I,2) = -T * C2 * X2 - X1 * BD(1,2) - X2 * BD(2,2) 190 CONTINUE ELSEIF (PBNO .EQ. 20) THEN C ... Modified Meyer DO 200 I = 1, M C = .1D0 * (PBDATA(I) + X(3)) U = EXP(X(2)/C - 13D0) F(I) = X(1) * U - PBDATA(M+I) DF(I,1) = U DF(I,2) = X(1)/C*U DF(I,3) = -.1D0 * DF(I,2) * X(2) / C 200 CONTINUE ELSEIF (PBNO .EQ. 21) THEN C ... Separated Meyer C ... Compute coefficient EE = 0D0 C2 = 0D0 DO 210 I = 1, M EI = EXP(X(1)/(PBDATA(I) + X(2))) EE = EE + EI**2 C2 = C2 + EI * PBDATA(M+I) 210 CONTINUE C1 = C2 / EE C ... Partial derivatives BY(1) = 0D0 BY(2) = 0D0 DO 220 I = 1, M D1 = PBDATA(I) + X(2) EI = EXP(X(1)/D1) D2 = EI * (PBDATA(M+I) - 2D0*C1*EI)/D1 BY(1) = BY(1) + D2 BY(2) = BY(2) - X(1)*D2/D1 220 CONTINUE BY(1) = BY(1)/EE BY(2) = BY(2)/EE C ... Function and Jacobian DO 230 I = 1, M D1 = PBDATA(I) + X(2) D2 = C1/D1 EI = EXP(X(1)/D1) F(I) = C1*EI - PBDATA(M+I) DF(I,1) = EI*(BY(1) + D2) DF(I,2) = EI*(BY(2) - X(1)*D2/D1) 230 CONTINUE ENDIF RETURN END C SUBROUTINE DCONST(CONST, N, S) ************************************************************************ * S(i) := const , i = 1,...,n * Hans Bruun Nielsen, Numerisk Institut, DTH. 92.10.05 / 00.09.22 ************************************************************************ INTEGER N, I,M,MP1 DOUBLE PRECISION CONST,S(*) C IF (N .LE. 0) RETURN M = MOD(N, 5) IF (M .EQ. 0) GOTO 20 DO 10 I = 1, M 10 S(I) = CONST 20 MP1 = M + 1 DO 30 I = MP1, N, 5 S(I) = CONST S(I+1) = CONST S(I+2) = CONST S(I+3) = CONST S(I+4) = CONST 30 CONTINUE RETURN ****************************** end of DCONST ************************ END