00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019 c
00020 c Least-squares fitting by singular-value decomposition.
00021 c From "Numerical Recipes," section 14.3, pp. 515-519, with
00022 c minor modifications.
00023 c
00024 c Perform a fit to the NDATA points Y(X), with weight 1/SIG, returning
00025 c the MA coefficients A of the fitting function, in terms of the basis
00026 c functions defined by the user-supplied function FUNCS. For now, the
00027 c input arrays U, V, and W, of dimension (MP,NP), (NP,NP), and NP,
00028 c respectively, are for workspace. We require MP >= NDATA, NP >= MA.
00029 c
00030
00031 c
00032 c Polynomial basis functions:
00033 c --------------------------
00034 c
00035 subroutine poly(x,f,n)
00036 save
00037 dimension f(n)
00038 c
00039 f(1) = 1.
00040 do 10 i=2,n
00041 10 f(i) = x*f(i-1)
00042 c
00043 return
00044 end
00045
00046 c
00047 c Trigonometric basis function1s:
00048 c -----------------------------
00049 c
00050 subroutine trig(x,f,n)
00051 save
00052 dimension f(n)
00053 c
00054 do 10 i=1,n,2
00055 10 f(i) = cos(((i-1)/2)*x)
00056 do 20 i=2,n-1,2
00057 20 f(i) = sin((i/2)*x)
00058 c
00059 return
00060 end
00061
00062
00063
00064 SUBROUTINE SVDFIT(X,Y,SIG,NDATA,A,MA,U,V,W,B,MP,NP,CHISQ,FUNCS)
00065 save
00066 c
00067 PARAMETER(MMAX=21,TOL=1.E-5)
00068 DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),V(NP,NP),
00069 * U(MP,NP),W(NP),B(MP),AFUNC(MMAX)
00070 c
00071 external funcs
00072 c
00073 DO 95000 I=1,NDATA
00074 CALL FUNCS(X(I),AFUNC,MA)
00075 TMP=1./SIG(I)
00076 DO 95001 J=1,MA
00077 U(I,J)=AFUNC(J)*TMP
00078 95001 CONTINUE
00079 B(I)=Y(I)*TMP
00080 95000 CONTINUE
00081 c
00082 c Pad U with blank rows, if necessary (see NR, p. 60)...
00083 c
00084 do 100 i=ndata+1,ma
00085 do 100 j=1,ma
00086 100 u(i,j) = 0.
00087 c
00088 CALL SVDCMP(U,max(ma,NDATA),MA,MP,NP,W,V)
00089 WMAX=0.
00090 DO 95002 J=1,MA
00091 IF(W(J).GT.WMAX)WMAX=W(J)
00092 95002 CONTINUE
00093 THRESH=TOL*WMAX
00094 DO 95003 J=1,MA
00095 IF(W(J).LT.THRESH)W(J)=0.
00096 95003 CONTINUE
00097 CALL SVBKSB(U,W,V,NDATA,MA,MP,NP,B,A)
00098 CHISQ=0.
00099 DO 95004 I=1,NDATA
00100 CALL FUNCS(X(I),AFUNC,MA)
00101 SUM=0.
00102 DO 95005 J=1,MA
00103 SUM=SUM+A(J)*AFUNC(J)
00104 95005 CONTINUE
00105 CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
00106 95004 CONTINUE
00107 RETURN
00108 END
00109
00110
00111
00112 SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V)
00113 save
00114 c
00115 PARAMETER (NMAX=1000)
00116 DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX)
00117 c
00118 G=0.0
00119 SCALE=0.0
00120 ANORM=0.0
00121 DO 95000 I=1,N
00122 L=I+1
00123 RV1(I)=SCALE*G
00124 G=0.0
00125 S=0.0
00126 SCALE=0.0
00127 IF (I.LE.M) THEN
00128 DO 95001 K=I,M
00129 SCALE=SCALE+ABS(A(K,I))
00130 95001 CONTINUE
00131 IF (SCALE.NE.0.0) THEN
00132 DO 95002 K=I,M
00133 A(K,I)=A(K,I)/SCALE
00134 S=S+A(K,I)*A(K,I)
00135 95002 CONTINUE
00136 F=A(I,I)
00137 G=-SIGN(SQRT(S),F)
00138 H=F*G-S
00139 A(I,I)=F-G
00140 IF (I.NE.N) THEN
00141 DO 95003 J=L,N
00142 S=0.0
00143 DO 95004 K=I,M
00144 S=S+A(K,I)*A(K,J
00145 + )
00146 95004 CONTINUE
00147 F=S/H
00148 DO 95005 K=I,M
00149 A(K,J)=A(K,J)+F*
00150 + A(K,I)
00151 95005 CONTINUE
00152 95003 CONTINUE
00153 ENDIF
00154 DO 95006 K= I,M
00155 A(K,I)=SCALE*A(K,I)
00156 95006 CONTINUE
00157 ENDIF
00158 ENDIF
00159 W(I)=SCALE *G
00160 G=0.0
00161 S=0.0
00162 SCALE=0.0
00163 IF ((I.LE.M).AND.(I.NE.N)) THEN
00164 DO 95007 K=L,N
00165 SCALE=SCALE+ABS(A(I,K))
00166 95007 CONTINUE
00167 IF (SCALE.NE.0.0) THEN
00168 DO 95008 K=L,N
00169 A(I,K)=A(I,K)/SCALE
00170 S=S+A(I,K)*A(I,K)
00171 95008 CONTINUE
00172 F=A(I,L)
00173 G=-SIGN(SQRT(S),F)
00174 H=F*G-S
00175 A(I,L)=F-G
00176 DO 95009 K=L,N
00177 RV1(K)=A(I,K)/H
00178 95009 CONTINUE
00179 IF (I.NE.M) THEN
00180 DO 95010 J=L,M
00181 S=0.0
00182 DO 95011 K=L,N
00183 S=S+A(J,K)*A(I,K
00184 + )
00185 95011 CONTINUE
00186 DO 95012 K=L,N
00187 A(J,K)=A(J,K)+S*
00188 + RV1(K)
00189 95012 CONTINUE
00190 95010 CONTINUE
00191 ENDIF
00192 DO 95013 K=L,N
00193 A(I,K)=SCALE*A(I,K)
00194 95013 CONTINUE
00195 ENDIF
00196 ENDIF
00197 ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
00198 95000 CONTINUE
00199 DO 95014 I=N,1,-1
00200 IF (I.LT.N) THEN
00201 IF (G.NE.0.0) THEN
00202 DO 95015 J=L,N
00203 V(J,I)=(A(I,J)/A(I,L))/G
00204 95015 CONTINUE
00205 DO 95016 J=L,N
00206 S=0.0
00207 DO 95017 K=L,N
00208 S=S+A(I,K)*V(K,J)
00209 95017 CONTINUE
00210 DO 95018 K=L,N
00211 V(K,J)=V(K,J)+S*V(K,I)
00212 95018 CONTINUE
00213 95016 CONTINUE
00214 ENDIF
00215 DO 95019 J=L,N
00216 V(I,J)=0.0
00217 V(J,I)=0.0
00218 95019 CONTINUE
00219 ENDIF
00220 V(I,I)=1.0
00221 G=RV1(I)
00222 L=I
00223 95014 CONTINUE
00224 DO 95020 I=N,1,-1
00225 L=I+1
00226 G=W(I)
00227 IF (I.LT.N) THEN
00228 DO 95021 J=L,N
00229 A(I,J)=0.0
00230 95021 CONTINUE
00231 ENDIF
00232 IF (G.NE.0.0) THEN
00233 G=1.0/G
00234 IF (I.NE.N) THEN
00235 DO 95022 J=L,N
00236 S=0.0
00237 DO 95023 K=L,M
00238 S=S+A(K,I)*A(K,J)
00239 95023 CONTINUE
00240 F=(S/A(I,I))*G
00241 DO 95024 K=I,M
00242 A(K,J)=A(K,J)+F*A(K,I)
00243 95024 CONTINUE
00244 95022 CONTINUE
00245 ENDIF
00246 DO 95025 J=I,M
00247 A(J,I)=A(J,I)*G
00248 95025 CONTINUE
00249 ELSE
00250 DO 95026 J= I,M
00251 A(J,I)=0.0
00252 95026 CONTINUE
00253 ENDIF
00254 A(I,I)=A(I,I)+1.0
00255 95020 CONTINUE
00256 DO 95027 K=N,1,-1
00257 DO 95028 ITS=1,30
00258 DO 95029 L=K,1,-1
00259 NM=L-1
00260 IF ((ABS(RV1(L))+ANORM).EQ.ANORM) GO TO
00261 + 2
00262 IF ((ABS(W(NM))+ANORM).EQ.ANORM) GO TO
00263 + 1
00264 95029 CONTINUE
00265 1 C=0.0
00266 S=1.0
00267 DO 95030 I=L,K
00268 F=S*RV1(I)
00269 IF ((ABS(F)+ANORM).NE.ANORM) THEN
00270 G=W(I)
00271 H=SQRT(F*F+G*G)
00272 W(I)=H
00273 H=1.0/H
00274 C= (G*H)
00275 S=-(F*H)
00276 DO 95031 J=1,M
00277 Y=A(J,NM)
00278 Z=A(J,I)
00279 A(J,NM)=(Y*C)+(Z*S)
00280 A(J,I)=-(Y*S)+(Z*C)
00281 95031 CONTINUE
00282 ENDIF
00283 95030 CONTINUE
00284 2 Z=W(K)
00285 IF (L.EQ.K) THEN
00286 IF (Z.LT.0.0) THEN
00287 W(K)=-Z
00288 DO 95032 J=1,N
00289 V(J,K)=-V(J,K)
00290 95032 CONTINUE
00291 ENDIF
00292 GO TO 3
00293 ENDIF
00294 IF (ITS.EQ.30)
00295 & stop 'No convergence in 30 iterations'
00296 X=W(L)
00297 NM=K-1
00298 Y=W(NM)
00299 G=RV1(NM)
00300 H=RV1(K)
00301 F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
00302 G=SQRT(F*F+1.0)
00303 F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X
00304 C=1.0
00305 S=1.0
00306 DO 95033 J=L,NM
00307 I=J+1
00308 G=RV1(I)
00309 Y=W(I)
00310 H=S*G
00311 G=C*G
00312 Z=SQRT(F*F+H*H)
00313 RV1(J)=Z
00314 C=F/Z
00315 S=H/Z
00316 F= (X*C)+(G*S)
00317 G=-(X*S)+(G*C)
00318 H=Y*S
00319 Y=Y*C
00320 DO 95034 NM=1,N
00321 X=V(NM,J)
00322 Z=V(NM,I)
00323 V(NM,J)= (X*C)+(Z*S)
00324 V(NM,I)=-(X*S)+(Z*C)
00325 95034 CONTINUE
00326 Z=SQRT(F*F+H*H)
00327 W(J)=Z
00328 IF (Z.NE.0.0) THEN
00329 Z=1.0/Z
00330 C=F*Z
00331 S=H*Z
00332 ENDIF
00333 F= (C*G)+(S*Y)
00334 X=-(S*G)+(C*Y)
00335 DO 95035 NM=1,M
00336 Y=A(NM,J)
00337 Z=A(NM,I)
00338 A(NM,J)= (Y*C)+(Z*S)
00339 A(NM,I)=-(Y*S)+(Z*C)
00340 95035 CONTINUE
00341 95033 CONTINUE
00342 RV1(L)=0.0
00343 RV1(K)=F
00344 W(K)=X
00345 95028 CONTINUE
00346 3 CONTINUE
00347 95027 CONTINUE
00348 RETURN
00349 END
00350
00351
00352
00353 SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X)
00354 save
00355 c
00356 PARAMETER (NMAX=1000)
00357 DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(NMAX)
00358 c
00359 DO 95000 J=1,N
00360 S=0.
00361 IF(W(J).NE.0.)THEN
00362 DO 95001 I=1,M
00363 S=S+U(I,J)*B(I)
00364 95001 CONTINUE
00365 S=S/W(J)
00366 ENDIF
00367 TMP(J)=S
00368 95000 CONTINUE
00369 DO 95002 J=1,N
00370 S=0.
00371 DO 95003 JJ=1,N
00372 S=S+V(J,JJ)*TMP(JJ)
00373 95003 CONTINUE
00374 X(J)=S
00375 95002 CONTINUE
00376 RETURN
00377 END
