| 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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame