1c program DRSNLAGU 2c>> 1997-06-18 DRSNLAGU Krogh Changes to improve C portability. 3c>> 1994-11-02 DRSNLAGU Krogh Changes to use M77CON 4c>> 1994-09-14 DRSNLAGU CLL Set IV(OUTLEV) = 0 for comparing output. 5c>> 1992-04-13 CLL Rename and reorder common block [D/S]KEY. 6c>> 1992-02-03 CLL @ JPL 7c>> 1990-07-02 CLL @ JPL 8c>> 1990-06-27 CLL @ JPL 9c>> 1990-06-14 CLL @ JPL 10c>> 1990-03-28 CLL @ JPL 11c Demo driver for SNLAGU. A variant of the nonlinear LS code NL2SOL. 12c SNLAGU requires values of the function and the Jacobian matrix. 13c ------------------------------------------------------------------ 14c--S replaces "?": DR?NLAGU, ?NLAGU, ?CALCR, ?CALCJ, ?IVSET, ?KEY 15c ------------------------------------------------------------------ 16 external SCALCR, SCALCJ 17 integer LIV, LV, MC, MDATA, NC, NDATA 18 parameter(MDATA = 30, MC = 7) 19 parameter(LIV = 82 + MC) 20 parameter(LV = 105 + MC*(MDATA + 2*MC + 17) + 2*MDATA) 21 integer IV(LIV) 22 integer F, COVPRT, OUTLEV, SOLPRT, STATPR, X0PRT 23 parameter(F=10) 24 parameter(COVPRT=14, OUTLEV=19, SOLPRT=22, STATPR=23, X0PRT=24) 25 real COEF(MC), DOF, V(LV) 26c ------------------------------------------------------------------ 27 NDATA = MDATA 28 NC = MC 29 COEF(1) = 5.0e0 30 COEF(2) = 10.0e0 31 COEF(3) = 0.5e0 32 COEF(4) = 0.5e0 33 COEF(5) = 0.5e0 34 COEF(6) = 0.5e0 35 COEF(7) = 0.5e0 36 IV(1) = 0 37 38 print'(1x,a)', 39 * 'Program DRSNLAGU.. Demo driver for SNLAGU.', 40 * ' A variant of NL2SOL.', 41 * ' SNLAGU requires values of the function and the Jacobian.', 42 * ' ', 43 * 'Sample problem is a nonlinear curve fit to data.', 44 * 'Model function is C3 + C4 * cos(C1*t) + C5 * sin(C1*t) +', 45 * ' C6 * cos(C2*t) + C7 * sin(C2*t) + Noise', 46 * 'Data generated using', 47 * '(C1, ..., C7) = (6, 9, 1, 0.5, 0.4, 0.2, 0.1)', 48 * 'and Gaussian noise with mean 0 and', 49 * 'sample standard deviation 0.001', 50 * ' ' 51 52 call SIVSET(1, IV, LIV, LV, V) 53 IV( X0PRT) = 1 54 IV(OUTLEV) = 0 55 IV(STATPR) = 1 56 IV(SOLPRT) = 1 57 IV(COVPRT) = 1 58 59 call SNLAGU(NDATA, NC, COEF, SCALCR, SCALCJ, IV, LIV, LV, V) 60 61 DOF = max(NDATA - NC, 1) 62 print'(1x/1x,a,g12.4)', 63 * 'SIGFAC: sqrt((2 * V(F))/DOF) =', 64 * sqrt(2.0e0 * V(F)/DOF) 65 stop 66 end 67c ================================================================== 68 subroutine SCALCR(NDATA, NC, C, NCOUNT, RVEC) 69c Function evaluation to test nonlinear least squares computation. 70c Illustrates saving results in common between SCALCR and SCALCJ 71c to avoid recalculation of common subexpressions. 72c ------------------------------------------------------------------ 73 common/SKEY/C1, C2, S1, S2, KEY 74 save /SKEY/ 75 integer I, KEY, MDATA, NCOUNT, NDATA, NC 76 parameter(MDATA = 30) 77 real C(NC), C1(MDATA), C2(MDATA), DEL, RVEC(NDATA) 78 real S1(MDATA), S2(MDATA), T, YDATA(MDATA) 79 data YDATA / 80 * 1.700641e0, 1.793512e0, 1.838309e0, 1.838416e0, 1.792204e0, 81 * 1.700501e0, 1.579804e0, 1.426268e0, 1.260724e0, 1.084901e0, 82 * 0.917094e0, 0.761920e0, 0.627304e0, 0.522146e0, 0.446645e0, 83 * 0.404920e0, 0.392033e0, 0.409622e0, 0.453045e0, 0.510765e0, 84 * 0.584554e0, 0.663109e0, 0.747613e0, 0.829439e0, 0.908496e0, 85 * 0.983178e0, 1.051046e0, 1.114072e0, 1.171746e0, 1.227823e0/ 86c ------------------------------------------------------------------ 87 T = 0.0E0 88 DEL = 1.0E0 / 29.0E0 89 KEY = NCOUNT 90 do 10 I = 1,NDATA 91 C1(I) = cos(C(1)*T) 92 S1(I) = sin(C(1)*T) 93 C2(I) = cos(C(2)*T) 94 S2(I) = sin(C(2)*T) 95 RVEC(I) = C(3) + C(4)*C1(I) + C(5)*S1(I) + 96 * C(6)*C2(I) + C(7)*S2(I) - 97 * YDATA(I) 98 T = T + DEL 99 10 continue 100 return 101 end 102c ================================================================== 103 subroutine SCALCJ(NDATA, NC, C, NCOUNT, AJAC) 104c Jacobian evaluation to test nonlinear least squares computation. 105c ------------------------------------------------------------------ 106 common/SKEY/C1, C2, S1, S2, KEY 107 save /SKEY/ 108 integer I, KEY, MDATA, NCOUNT, NDATA, NC 109 parameter(MDATA = 30) 110 real AJAC(NDATA,NC), C(NC), C1(MDATA), C2(MDATA) 111 real DEL, S1(MDATA), S2(MDATA), T 112c ------------------------------------------------------------------ 113 T = 0.0E0 114 DEL = 1.0E0 / 29.0E0 115 if(NCOUNT .eq. KEY) then 116 do 10 I = 1,NDATA 117 AJAC(I,1) = -C(4)*S1(I)*T + C(5)*C1(I)*T 118 AJAC(I,2) = -C(6)*S2(I)*T + C(7)*C2(I)*T 119 AJAC(I,3) = 1.0E0 120 AJAC(I,4) = C1(I) 121 AJAC(I,5) = S1(I) 122 AJAC(I,6) = C2(I) 123 AJAC(I,7) = S2(I) 124 T = T + DEL 125 10 continue 126 else 127 do 20 I = 1,NDATA 128 AJAC(I,3) = 1.0E0 129 AJAC(I,4) = cos(C(1)*T) 130 AJAC(I,5) = sin(C(1)*T) 131 AJAC(I,6) = cos(C(2)*T) 132 AJAC(I,7) = sin(C(2)*T) 133 AJAC(I,1) = -C(4)*AJAC(I,5)*T + C(5)*AJAC(I,4)*T 134 AJAC(I,2) = -C(6)*AJAC(I,7)*T + C(7)*AJAC(I,6)*T 135 T = T + DEL 136 20 continue 137 endif 138 return 139 end 140