]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:20 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE MINSQ (M,N,F,X,E,IPRINT,NFUN ,NW,W,COV,XSTEP) | |
11 | DIMENSION F(M),X(N),E(N),W(NW) | |
12 | IFLAG=1 | |
13 | ESCALE=1. | |
14 | MAXFUN=NFUN | |
15 | MPLUSN=M+N | |
16 | KST=N+MPLUSN | |
17 | NPLUS=N+1 | |
18 | KINV=NPLUS*(MPLUSN+1) | |
19 | KSTORE=KINV-MPLUSN-1 | |
20 | CALL FCN1(M,N,F,X,IFLAG,NW,W,KEND,KEND2N) | |
21 | IFLAG=4 | |
22 | NN=N+N | |
23 | K=NN | |
24 | DO 1 I=1,M | |
25 | K=K+1 | |
26 | W(K)=F(I) | |
27 | 1 CONTINUE | |
28 | IINV=2 | |
29 | K=KST | |
30 | I=1 | |
31 | 2 X(I)=X(I)+E(I) | |
32 | CALL FCN1(M,N,F,X,IFLAG,NW,W,KEND,KEND2N) | |
33 | X(I)=X(I)-E(I) | |
34 | DO 3 J=1,N | |
35 | K=K+1 | |
36 | W(K)=0. | |
37 | W(J)=0. | |
38 | 3 CONTINUE | |
39 | SUM=0. | |
40 | KK=NN | |
41 | DO 4 J=1,M | |
42 | KK=KK+1 | |
43 | F(J)=F(J)-W(KK) | |
44 | SUM=SUM+F(J)*F(J) | |
45 | 4 CONTINUE | |
46 | IF (SUM) 5,5,6 | |
47 | 5 WRITE(6,7)I | |
48 | 7 FORMAT(5X,'E(',I3,') UNREASONABLY SMALL') | |
49 | RETURN | |
50 | 6 SUM=1./SQRT (SUM) | |
51 | J=K-N+I | |
52 | W(J)=E(I)*SUM | |
53 | DO 9 J=1,M | |
54 | K=K+1 | |
55 | W(K)=F(J)*SUM | |
56 | KK=NN+J | |
57 | DO 11 II=1,I | |
58 | KK=KK+MPLUSN | |
59 | W(II)=W(II)+W(KK)*W(K) | |
60 | 11 CONTINUE | |
61 | 9 CONTINUE | |
62 | ILESS=I-1 | |
63 | IGAMAX=N+I-1 | |
64 | INCINV=N-ILESS | |
65 | INCINP=INCINV+1 | |
66 | IF (ILESS) 13,13,14 | |
67 | 13 W(KINV)=1. | |
68 | GO TO 15 | |
69 | 14 B=1. | |
70 | DO 16 J=NPLUS,IGAMAX | |
71 | W(J)=0. | |
72 | 16 CONTINUE | |
73 | KK=KINV | |
74 | DO 17 II=1,ILESS | |
75 | IIP=II+N | |
76 | W(IIP)=W(IIP)+W(KK)*W(II) | |
77 | JL=II+1 | |
78 | IF (JL-ILESS) 18,18,19 | |
79 | 18 DO 20 JJ=JL,ILESS | |
80 | KK=KK+1 | |
81 | JJP=JJ+N | |
82 | W(IIP)=W(IIP)+W(KK)*W(JJ) | |
83 | W(JJP)=W(JJP)+W(KK)*W(II) | |
84 | 20 CONTINUE | |
85 | 19 B=B-W(II)*W(IIP) | |
86 | KK=KK+INCINP | |
87 | 17 CONTINUE | |
88 | B=1./B | |
89 | KK=KINV | |
90 | DO 21 II=NPLUS,IGAMAX | |
91 | BB=-B*W(II) | |
92 | DO 22 JJ=II,IGAMAX | |
93 | W(KK)=W(KK)-BB*W(JJ) | |
94 | KK=KK+1 | |
95 | 22 CONTINUE | |
96 | W(KK)=BB | |
97 | KK=KK+INCINV | |
98 | 21 CONTINUE | |
99 | W(KK)=B | |
100 | 15 GO TO (27,24),IINV | |
101 | 24 I=I+1 | |
102 | IF (I-N) 2,2,25 | |
103 | 25 IINV=1 | |
104 | FF=0. | |
105 | KL=NN | |
106 | DO 26 I=1,M | |
107 | KL=KL+1 | |
108 | F(I)=W(KL) | |
109 | FF=FF+F(I)*F(I) | |
110 | 26 CONTINUE | |
111 | ICONT = 1 | |
112 | ISS=1 | |
113 | MC=N+1 | |
114 | IPP=IPRINT*(IPRINT-1) | |
115 | ITC=0 | |
116 | IPS=1 | |
117 | IPC=0 | |
118 | 27 IPC=IPC-IPRINT | |
119 | IF (IPC) 28,29,29 | |
120 | 28 WRITE(6,30)ITC,MC,FF | |
121 | 30 FORMAT(//5X,'ITERATION',I4,I9,' CALLS OF FCN ',5X,'F=',E24.14) | |
122 | WRITE(6,31)(X(I),I=1,N) | |
123 | 31 FORMAT(5X,'VARIABLES',/(5E24.14)) | |
124 | WRITE(6,32)(F(I),I=1,M) | |
125 | 32 FORMAT(5X,'FUNCTIONS',/(5E24.14)) | |
126 | IPC=IPP | |
127 | GO TO (29,33),IPS | |
128 | 29 GO TO (34,35),ICONT | |
129 | 35 IF (CHANGE-1.) 10,10,36 | |
130 | 10 CONTINUE | |
131 | C*UL 37 WRITE(6,38) | |
132 | WRITE(6,38) | |
133 | 38 FORMAT(//5X,' FINAL VALUES OF FUNCTIONS AND VARIABLES') | |
134 | IPS=2 | |
135 | GO TO 28 | |
136 | 33 IFLAG=3 | |
137 | CALL FCN1(M,N,F,X,IFLAG,NW,W,KEND,KEND2N) | |
138 | IF(COV.EQ.0.) RETURN | |
139 | IWC1 = 2*N+N*(N+M) | |
140 | IWC=IWC1 | |
141 | DO 91 I=1,N | |
142 | DO 91 J=1,I | |
143 | IWC=IWC+1 | |
144 | W(IWC)=0. | |
145 | DO 90 MR=1,N | |
146 | JCOV=N+J+MR*(N+M) | |
147 | ICOV=N+I+MR*(N+M) | |
148 | 90 W(IWC)=W(IWC)+W(ICOV)*W(JCOV) | |
149 | 91 CONTINUE | |
150 | WRITE(6,1000) | |
151 | 1000 FORMAT(///40X,'VARIANCE-COVARIANCE MATRIX'///) | |
152 | IB=IWC1 | |
153 | DO 92 I=1,N | |
154 | IA=IB+1 | |
155 | IB=IB+I | |
156 | 92 WRITE(6,1001)(W(J),J=IA,IB) | |
157 | 1001 FORMAT(10X,5G20.6//) | |
158 | RETURN | |
159 | 36 ICONT=1 | |
160 | 34 ITC=ITC+1 | |
161 | K=N | |
162 | KK=KST | |
163 | DO 39 I=1,N | |
164 | K=K+1 | |
165 | W(K)=0. | |
166 | KK=KK+N | |
167 | W(I)=0. | |
168 | DO 40 J=1,M | |
169 | KK=KK+1 | |
170 | W(I)=W(I)+W(KK)*F(J) | |
171 | 40 CONTINUE | |
172 | 39 CONTINUE | |
173 | DM=0. | |
174 | K=KINV | |
175 | DO 41 II=1,N | |
176 | IIP=II+N | |
177 | W(IIP)=W(IIP)+W(K)*W(II) | |
178 | JL=II+1 | |
179 | IF (JL-N) 42,42,43 | |
180 | 42 DO 44 JJ=JL,N | |
181 | JJP=JJ+N | |
182 | K=K+1 | |
183 | W(IIP)=W(IIP)+W(K)*W(JJ) | |
184 | W(JJP)=W(JJP)+W(K)*W(II) | |
185 | 44 CONTINUE | |
186 | K=K+1 | |
187 | 43 IF (DM-ABS (W(II)*W(IIP))) 45,41,41 | |
188 | 45 DM=ABS (W(II)*W(IIP)) | |
189 | KL=II | |
190 | 41 CONTINUE | |
191 | II=N+MPLUSN*KL | |
192 | CHANGE=0. | |
193 | DO 46 I=1,N | |
194 | JL=N+I | |
195 | W(I)=0. | |
196 | DO 47 J=NPLUS,NN | |
197 | JL=JL+MPLUSN | |
198 | W(I)=W(I)+W(J)*W(JL) | |
199 | 47 CONTINUE | |
200 | II=II+1 | |
201 | W(II)=W(JL) | |
202 | W(JL)=X(I) | |
203 | IF (ABS (E(I)*CHANGE)-ABS (W(I))) 48,48,46 | |
204 | 48 CHANGE=ABS (W(I)/E(I)) | |
205 | 46 CONTINUE | |
206 | DO 49 I=1,M | |
207 | II=II+1 | |
208 | JL=JL+1 | |
209 | W(II)=W(JL) | |
210 | W(JL)=F(I) | |
211 | 49 CONTINUE | |
212 | FC=FF | |
213 | ACC=0.1/CHANGE | |
214 | IT=3 | |
215 | XC=0. | |
216 | XL=0. | |
217 | IS=3 | |
218 | IF (CHANGE-1.) 50,50,51 | |
219 | 50 ICONT=2 | |
220 | 51 CALL VD01A (IT,XC,FC,20,ACC,0.1,XSTEP) | |
221 | GO TO (52,53,53,53),IT | |
222 | 52 MC=MC+1 | |
223 | IF (MC-MAXFUN) 54,54,55 | |
224 | 55 WRITE(6,56)MAXFUN | |
225 | 56 FORMAT(5X,I6,' CALLS OF FCN') | |
226 | ISS=2 | |
227 | GO TO 53 | |
228 | 54 XL=XC-XL | |
229 | DO 57 J=1,N | |
230 | X(J)=X(J)+XL*W(J) | |
231 | 57 CONTINUE | |
232 | XL=XC | |
233 | CALL FCN1(M,N,F,X,IFLAG,NW,W,KEND,KEND2N) | |
234 | FC=0. | |
235 | DO 58 J=1,M | |
236 | FC=FC+F(J)*F(J) | |
237 | 58 CONTINUE | |
238 | GO TO (59,59,60),IS | |
239 | 60 K=N | |
240 | IF (FC-FF) 61,51,62 | |
241 | 61 IS=2 | |
242 | FMIN=FC | |
243 | FSEC=FF | |
244 | GO TO 63 | |
245 | 62 IS=1 | |
246 | FMIN=FF | |
247 | FSEC=FC | |
248 | GO TO 63 | |
249 | 59 IF (FC-FSEC) 64,51,51 | |
250 | 64 K=KSTORE | |
251 | GO TO (75,74),IS | |
252 | 75 K=N | |
253 | 74 IF (FC-FMIN) 65,51,66 | |
254 | 66 FSEC=FC | |
255 | GO TO 63 | |
256 | 65 IS=3-IS | |
257 | FSEC=FMIN | |
258 | FMIN=FC | |
259 | 63 DO 67 J=1,N | |
260 | K=K+1 | |
261 | W(K)=X(J) | |
262 | 67 CONTINUE | |
263 | DO 68 J=1,M | |
264 | K=K+1 | |
265 | W(K)=F(J) | |
266 | 68 CONTINUE | |
267 | GO TO 51 | |
268 | 53 K=KSTORE | |
269 | KK=N | |
270 | GO TO (69,70,69),IS | |
271 | 70 K=N | |
272 | KK=KSTORE | |
273 | 69 SUM=0. | |
274 | DM=0. | |
275 | JJ=KSTORE | |
276 | DO 71 J=1,N | |
277 | K=K+1 | |
278 | KK=KK+1 | |
279 | JJ=JJ+1 | |
280 | X(J)=W(K) | |
281 | W(JJ)=W(K)-W(KK) | |
282 | 71 CONTINUE | |
283 | DO 72 J=1,M | |
284 | K=K+1 | |
285 | KK=KK+1 | |
286 | JJ=JJ+1 | |
287 | F(J)=W(K) | |
288 | W(JJ)=W(K)-W(KK) | |
289 | SUM=SUM+W(JJ)*W(JJ) | |
290 | DM=DM+F(J)*W(JJ) | |
291 | 72 CONTINUE | |
292 | GO TO (73,10),ISS | |
293 | 73 J=KINV | |
294 | KK=NPLUS-KL | |
295 | DO 76 I=1,KL | |
296 | K=J+KL-I | |
297 | J=K+KK | |
298 | W(I)=W(K) | |
299 | W(K)=W(J-1) | |
300 | 76 CONTINUE | |
301 | IF (KL-N) 77,78,78 | |
302 | 77 KL=KL+1 | |
303 | JJ=K | |
304 | DO 79 I=KL,N | |
305 | K=K+1 | |
306 | J=J+NPLUS-I | |
307 | W(I)=W(K) | |
308 | W(K)=W(J-1) | |
309 | 79 CONTINUE | |
310 | W(JJ)=W(K) | |
311 | B=1./W(KL-1) | |
312 | W(KL-1)=W(N) | |
313 | GO TO 88 | |
314 | 78 B=1./W(N) | |
315 | 88 K=KINV | |
316 | DO 80 I=1,ILESS | |
317 | BB=B*W(I) | |
318 | DO 81 J=I,ILESS | |
319 | W(K)=W(K)-BB*W(J) | |
320 | K=K+1 | |
321 | 81 CONTINUE | |
322 | K=K+1 | |
323 | 80 CONTINUE | |
324 | IF (FMIN-FF) 82,83,83 | |
325 | 83 CHANGE=0. | |
326 | GO TO 84 | |
327 | 82 FF=FMIN | |
328 | CHANGE=ABS (XC)*CHANGE | |
329 | 84 XL=-DM/FMIN | |
330 | SUM=1./SQRT (SUM+DM*XL) | |
331 | K=KSTORE | |
332 | DO 85 I=1,N | |
333 | K=K+1 | |
334 | W(K)=SUM*W(K) | |
335 | W(I)=0. | |
336 | 85 CONTINUE | |
337 | DO 86 I=1,M | |
338 | K=K+1 | |
339 | W(K)=SUM*(W(K)+XL*F(I)) | |
340 | KK=NN+I | |
341 | DO 87 J=1,N | |
342 | KK=KK+MPLUSN | |
343 | W(J)=W(J)+W(KK)*W(K) | |
344 | 87 CONTINUE | |
345 | 86 CONTINUE | |
346 | GO TO 14 | |
347 | END |