1c program DRDNLSGU 2c>> 2001-05-24 DRDNLSGU Krogh Minor change for making .f90 version. 3c>> 1997-06-18 DRDNLSGU Krogh Changes to improve C portability. 4c>> 1994-11-02 DRDNLSGU Krogh Changes to use M77CON 5c>> 1994-09-14 DRDNLSGU CLL Set IV(OUTLEV) = 0 for comparing output. 6c>> 1992-04-13 CLL Rename and reorder common block [D/S]KEY. 7c>> 1992-02-03 CLL @ JPL 8c>> 1990-07-02 CLL @ JPL 9c>> 1990-06-27 CLL @ JPL 10c>> 1990-04-05 CLL @ JPL 11c>> 1990-03-29 CLL @ JPL 12c Demo driver for DNLSGU. A variant of the nonlinear LS code NL2SOL. 13c DNLSGU solves the "separable" problem. 14c DNLSGU requires values of the function and the Jacobian matrix. 15c Note: MDER is the number of ones in the array IND(). 16c ------------------------------------------------------------------ 17c--D replaces "?": DR?NLSGU, ?NLSGU, ?CALCA, ?CALCB, ?IVSET, ?KEY 18c ------------------------------------------------------------------ 19 external DCALCA, DCALCB 20 integer ITERM,IVAR,LIV,LV,MA,MB,MDATA,MDER,MLEN,NA,NB,NDATA 21 parameter(MDATA = 30, MA = 2, MB = 5) 22 parameter(MDER = 4) 23 parameter(LIV = 115 + MA + MB + 2*MDER) 24 parameter(MLEN = (MB+MA)*(MDATA+MB+MA+1)) 25 parameter(LV = 105 + MDATA*(MB+MDER+3) + 26 * MLEN + (MB*(MB+3))/2 + MA*(2*MA+17)) 27 integer IND(MB+1,MA), IV(LIV) 28 integer F, COVPRT, OUTLEV, SOLPRT, STATPR, X0PRT 29 parameter(F=10) 30 parameter(COVPRT=14, OUTLEV=19, SOLPRT=22, STATPR=23, X0PRT=24) 31 double precision ALF(MA), BET(MB), DOF, V(LV), YDATA(MDATA) 32 data ((IND(ITERM,IVAR),IVAR = 1,2),ITERM = 1,6)/ 33 * 0, 0, 34 * 1, 0, 35 * 1, 0, 36 * 0, 1, 37 * 0, 1, 38 * 0, 0/ 39 data YDATA / 40 * 1.700641d0, 1.793512d0, 1.838309d0, 1.838416d0, 1.792204d0, 41 * 1.700501d0, 1.579804d0, 1.426268d0, 1.260724d0, 1.084901d0, 42 * 0.917094d0, 0.761920d0, 0.627304d0, 0.522146d0, 0.446645d0, 43 * 0.404920d0, 0.392033d0, 0.409622d0, 0.453045d0, 0.510765d0, 44 * 0.584554d0, 0.663109d0, 0.747613d0, 0.829439d0, 0.908496d0, 45 * 0.983178d0, 1.051046d0, 1.114072d0, 1.171746d0, 1.227823d0/ 46c ------------------------------------------------------------------ 47 NDATA = MDATA 48 NA = MA 49 NB = MB 50 ALF(1) = 5.0d0 51 ALF(2) = 10.0d0 52 IV(1) = 0 53 54 print'(1x,a)', 55 * 'Program DRDNLSGU.. Demo driver for DNLSGU.', 56 * ' A variant of NL2SOL.', 57 * ' DNLSGU handles the Separable problem.', 58 * ' DNLSGU requires values of the function and the Jacobian.', 59 * ' ', 60 * 'Sample problem is a nonlinear curve fit to data.', 61 * 'Model function is B1 + B2 * cos(A1*t) + B3 * sin(A1*t) +', 62 * ' B4 * cos(A2*t) + B5 * sin(A2*t) + Noise', 63 * 'Data generated using', 64 * '(A1, A2, B1, ..., B5) = (6, 9, 1, 0.5, 0.4, 0.2, 0.1)', 65 * 'and Gaussian noise with mean 0 and', 66 * 'sample standard deviation 0.001', 67 * ' ' 68 69 call DIVSET(1, IV, LIV, LV, V) 70 IV( X0PRT) = 1 71 IV(OUTLEV) = 0 72 IV(STATPR) = 1 73 IV(SOLPRT) = 1 74 IV(COVPRT) = 1 75 76 call DNLSGU(NDATA, NA, NB, ALF, BET, YDATA, DCALCA, DCALCB, 77 * IND, NB+1, IV, LIV, LV, V) 78 79 DOF = max(NDATA - NA - NB, 1) 80 print'(1x/1x,a,g12.4)', 81 * 'SIGFAC: sqrt((2 * V(F))/DOF) =', 82 * sqrt(2.0d0 * V(F)/DOF) 83 stop 84 end 85c ================================================================== 86 subroutine DCALCA(NDATA, NA, NB, ALF, NCOUNT, PHI) 87c Test case for separable nonlinear least squares computation. 88c Computes MDATA x NB matrix PHI as a function of the 89c nonlinear parameters ALF(). 90c For J .le. NB the (I,J) term of PHI is the coefficient of the 91c linear coefficient B(J) in row I of the model. 92c In this example the model does not have a term that is not 93c multiplied by a linear coefficient. If such a term is present 94c then PHI must have an (NB+1)st column to hold this term. 95c This code Illustrates saving results in common between DCALCA and 96c DCALCB to avoid recalculation of common subexpressions. 97c ------------------------------------------------------------------ 98 common/DKEY/C1, C2, S1, S2, KEY 99 save /DKEY/ 100 integer I, KEY, MDATA, NA, NB, NCOUNT, NDATA 101 parameter(MDATA = 30) 102 double precision ALF(NA), C1(MDATA), C2(MDATA) 103 double precision DEL, PHI(NDATA,NB) 104 double precision S1(MDATA), S2(MDATA), T 105c ------------------------------------------------------------------ 106 T = 0.0D0 107 DEL = 1.0D0 / 29.0D0 108 KEY = NCOUNT 109 do 10 I = 1,NDATA 110 C1(I) = cos(ALF(1)*T) 111 S1(I) = sin(ALF(1)*T) 112 C2(I) = cos(ALF(2)*T) 113 S2(I) = sin(ALF(2)*T) 114 PHI(I,1) = 1.0D0 115 PHI(I,2) = C1(I) 116 PHI(I,3) = S1(I) 117 PHI(I,4) = C2(I) 118 PHI(I,5) = S2(I) 119 T = T + DEL 120 10 continue 121 return 122 end 123c ================================================================== 124 subroutine DCALCB(NDATA, NA, NB, ALF, NCOUNT, DER) 125c Test case for separable nonlinear least squares computation. 126c Computes the NDATA x NDER matrix DER. Here NDER is the number of 127c ones in the indicator array IND(). The columns of DER correspond 128c to nonzero entries of IND() traversed columnwise. 129c In this example NDER is 4 and the correspondence is as follows: 130c Col of DER: 1 2 3 4 131c Element of IND(): (2,1) (3,1) (4,2) (5,2) 132c In this example Row I of DER will be set to contain the values 133c of the four partial derivatives: 134c Partial of PHI(I,2) with respect to ALP(1) 135c Partial of PHI(I,3) with respect to ALP(1) 136c Partial of PHI(I,4) with respect to ALP(2) 137c Partial of PHI(I,5) with respect to ALP(2) 138c ------------------------------------------------------------------ 139 common/DKEY/C1, C2, S1, S2, KEY 140 save /DKEY/ 141 integer I, KEY, MDATA, NCOUNT, NDATA, NA, NB, NDER 142 parameter(MDATA = 30, NDER = 4) 143 double precision ALF(NA), C1(MDATA), C2(MDATA) 144 double precision DEL, DER(NDATA,NDER), S1(MDATA), S2(MDATA), T 145c ------------------------------------------------------------------ 146 T = 0.0D0 147 DEL = 1.0D0 / 29.0D0 148 if(NCOUNT .eq. KEY) then 149 do 10 I = 1,NDATA 150 DER(I,1) = -S1(I)*T 151 DER(I,2) = C1(I)*T 152 DER(I,3) = -S2(I)*T 153 DER(I,4) = C2(I)*T 154 T = T + DEL 155 10 continue 156 else 157 do 20 I = 1,NDATA 158 DER(I,1) = -sin(ALF(1)*T)*T 159 DER(I,2) = cos(ALF(1)*T)*T 160 DER(I,3) = -sin(ALF(2)*T)*T 161 DER(I,4) = cos(ALF(2)*T)*T 162 T = T + DEL 163 20 continue 164 endif 165 return 166 end 167