C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C R G D F I T - FITTING OF RARE-GAS ATOM POSITION IN THE DIMER BETWEEN C A MOLECULE AND A RARE GAS ATOM C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Program requires only A,B,C of the dimer and of the acceptor molecule, C and the masses of molecule and of the Rg atom to determine C the position of the Rg atom in the inertial axes of the molecule. C C see: C C 1/ Suenram, Fraser, Lovas, J.Chem.Phys. 89,6141(1988) - method C 2/ Kisiel, Fowler, Legon, J.Chem.Phys. 95,2283(1991) - application C of this program C C C Coordinates of the Rg atom are fitted to rotational constants of the C complex by least squares using singular value decomposition from C 'Numerical Recipes'. C C Since three Cartesians are fitted to three rotational constants of C the complex their reliability is critically dependent on C the magnitude of inertial vibration-rotation contributions - any C subtraction of such contributions will increase reliability C considerably C C C version 5.XII.2001 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c 1990: first version c 1998: modification to use SVD C 5.12.01: overhaul for adding to PROSPE C_____________________________________________________________________________ C C INPUT IS FROM USER SPECIFIED FILE CONTAINING THE FOLLOWING: C C - descriptive comment in the first line followed by a spacing line C - mass of the molecule, mass of the Rg atom (in amu) C - rotational constants (in MHz) of the molecule, defining the Z,X,Y C cartesian axes relative to which the polar coordinates of the Rg C are defined C For consistency with the SFL paper these should be specified in C the order A,B,C in which case the ZX plane of the molecule C will be the reference plane. C If it is desired for the Z axis to be perpendicular to the AB molecular C plane then the order should be C,A,B for which the XY plane will be C the reference plane. C - rotational constants A, B, C (in MHz) of the Rg...molecule dimer C - initial polar coordinates R, Theta, Phi of Rg in the inertial C frame of the molecule (in Angstroms and degrees) C C Only the first set of declarations found in the data is read in, C remainder of the file can be used for storage C C OPUTPUT IS WRITTEN TO FILE RGDFIT.OUT C C----------------------------------------------------------------------------- C C NOTE: The whole singular value decomposition package is in single C precision as adapted from program FIT. To save on changes C dimensionality declared in FIT is preserved and buffering C between REAL*8 internal working of M/P (as developed with C normal least squares) and SVD is provided. C C----------------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) EXTERNAL FDATA COMMON /TEN/MU,IMON/HBLK/T,S COMMON /BLKCF/DERSVD CHARACTER FILNAM*30,coment*79 REAL*8 MMON1,MMON2,MU,CONMON(3),CONCMP(3),IMON(3),ICMP(3), . T(3,3),S(3,3),CORD(3),DCORD(3),DMOFI(3),DER(3,3), . MICALC(3),DDCORD(3),CONCLC(3) REAL*4 XSVD(200),SIGSVD(200),YSVD(200),ASVD(20),USVD(200,20), . VSVD(20,20),WSVD(20),DERSVD(3,3),CHISQ EQUIVALENCE (R,CORD(1)),(THETA,CORD(2)),(PHI,CORD(3)) C DATA PI,TWOPI/3.141592654D0,6.283185307D0/ CONST=505379.07D0 CONV=57.29577951D0 NITER=0 DCORD(1)=0.0001 DCORD(2)=0.0001 DCORD(3)=0.0001 C WRITE(*,33) 33 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | R G D F I T - Fit of the position of a Rare-Gas ', * 'atom in a', * T79,'|'/ * ' | dimer with an acceptor molecule', * T79,'|'/ * ' |',76(1H_),'|'/' version 5.XII.2001',T64,'Zbigniew KISIEL'/) 1 FORMAT(1X,A) C C C...Read the data C 3 WRITE(*,2)'Name of file containing the data: ' 2 FORMAT(1X/1X,A\) READ(*,'(A)',ERR=3)FILNAM OPEN(3,FILE=FILNAM,status='old',ERR=3) READ(3,'(a/)')coment READ(3,*)MMON1,MMON2 READ(3,*)CONMON(1),CONMON(2),CONMON(3) READ(3,*)CONCMP(1),CONCMP(2),CONCMP(3) READ(3,*)CORD(1),CORD(2),CORD(3) CLOSE(3) CORD(2)=CORD(2)/CONV CORD(3)=CORD(3)/CONV C DO 4 I=1,3 IMON(I)=CONST/CONMON(I) ICMP(I)=CONST/CONCMP(I) CONCLC(I)=CONCMP(I) 4 CONTINUE MU=(MMON1*MMON2)/(MMON1+MMON2) C C...................................................................... C C FITTING LOOP C C...Set up and diagonalise the inertial tensor C open(4,file='rgdfit.out',status='unknown') WRITE(4,33) WRITE(*,1)' ' C 9 CALL NORMAL(1,THETA) CALL NORMAL(2,PHI) CALL MOFI(R,THETA,PHI) WRITE(*,21)R,THETA*CONV,PHI*CONV WRITE(4,21)R,THETA*CONV,PHI*CONV 21 FORMAT(1X,' R, Theta, PHI = ',3F12.5) WRITE(*,5)(CONST/T(I,I),I=1,3) WRITE(4,5)(CONST/T(I,I),I=1,3) 5 FORMAT(1X,'Calculated A,B,C = ',3F12.5) DO 6 I=1,3 MICALC(I)=T(I,I) DMOFI(I)=ICMP(I)-MICALC(I) 6 CONTINUE C C...Derivatives C DO 7 I=1,3 CORD(I)=CORD(I)+DCORD(I) CALL MOFI(R,THETA,PHI) CORD(I)=CORD(I)-DCORD(I) DO 8 J=1,3 DER(J,I)=(T(J,J)-MICALC(J))/DCORD(I) 8 CONTINUE 7 CONTINUE C NPTS=3 NCOEFF=3 DO 50 I=1,3 YSVD(I)=DMOFI(I) XSVD(I)=I SIGSVD(I)=1.D0 DO 50 J=1,3 DERSVD(I,J)=DER(I,J) 50 CONTINUE C C...Evaluate new coordinates C CALL SVDFIT(XSVD,YSVD,SIGSVD,NPTS,ASVD,NCOEFF,USVD, * VSVD,WSVD,200,20,CHISQ,FDATA) C DO 10 I=1,3 DDCORD(I)=ASVD(I) IF(NITER.GT.20)DDCORD(I)=DDCORD(I)/NITER 10 CONTINUE if(dabs(ddcord(2)).gt.0.5) . ddcord(2)=dsign(1.0d0,ddcord(2))*0.5d0 if(dabs(ddcord(3)).gt.0.5) . ddcord(3)=dsign(1.0d0,ddcord(3))*0.5d0 NITER=NITER+1 WRITE(*,11)NITER,DDCORD(1),(DDCORD(I)*CONV,I=2,3) WRITE(4,11)NITER,DDCORD(1),(DDCORD(I)*CONV,I=2,3) 11 FORMAT(1X/' Iteration no.',I3/ . 1X,'dR, dTheta, dPhi = ',3F12.5) DO 12 I=1,3 CORD(I)=CORD(I)+DDCORD(I) 12 CONTINUE C C...convergence test and exit C IF(NITER.EQ.1)GOTO 9 SUM=0.0 DO 22 I=1,3 SUM=SUM+DABS(CONCLC(I)-CONST/MICALC(I)) CONCLC(I)=CONST/MICALC(I) 22 CONTINUE IF(SUM.GT.0.001D0)GOTO 9 C C...................................................................... C C...output of final results C WRITE(*,31)COMENT WRITE(4,31)COMENT 31 FORMAT(1X/1x,78(1h-)/1x,a/1x,78(1H-) * //' THE FOUR SOLUTIONS FOR RARE-GAS ATOM POSITION ARE:'// * 22X,'R',7X,'THETA',5X,'PHI',10X,'X',9X,'Y',9X,'Z'/) C CALL NORMAL(1,THETA) CALL NORMAL(2,PHI) IF(PI-THETA.LT.THETA)THETA=PI-THETA IF(PI-PHI.LT.PHI)PHI=PI-PHI IF(TWOPI-PHI.LT.PHI)PHI=TWOPI-PHI TH=THETA*CONV PS=PHI*CONV CALL XYZAR(R,TH,PS,X,Y,Z) WRITE(*,30)' R,T0,P0:',R,TH,PS,X,Y,Z WRITE(4,30)' R,T0,P0:',R,TH,PS,X,Y,Z 30 FORMAT(1X,A,T17,F10.5,2F10.3,2X,3F10.5) C PHI1=PHI+PI CALL NORMAL(2,PHI1) PH1=PHI1*CONV CALL XYZAR(R,TH,PH1,X,Y,Z) WRITE(*,30)' R,T0,PI+P0:',R,TH,PH1,X,Y,Z WRITE(4,30)' R,T0,PI+P0:',R,TH,PH1,X,Y,Z C THETA1=PI-THETA CALL NORMAL(1,THETA1) PHI2=TWOPI-PHI CALL NORMAL(2,PHI2) TH1=THETA1*CONV PH2=PHI2*CONV CALL XYZAR(R,TH1,PH2,X,Y,Z) WRITE(*,30)'R,PI-T0,2*PI-P0:',R,TH1,PH2,X,Y,Z WRITE(4,30)'R,PI-T0,2*PI-P0:',R,TH1,PH2,X,Y,Z C PHI3=PI-PHI CALL NORMAL(2,PHI3) PH3=PHI3*CONV CALL XYZAR(R,TH1,PH3,X,Y,Z) WRITE(*,30)' R,PI-T0,PI-P0:',R,TH1,PH3,X,Y,Z WRITE(4,30)' R,PI-T0,PI-P0:',R,TH1,PH3,X,Y,Z WRITE(*,1)' ' WRITE(4,1)' ' C WRITE(*,35)' OBSERVED RG-MOLECULE CONSTANTS:',(CONCMP(I),I=1,3) WRITE(*,35)'CALCULATED RG-MOLECULE CONSTANTS:',(CONCLC(I),I=1,3) WRITE(*,35)' Chi Squared = SUM(I.o-I.c)**2 :',CHISQ WRITE(*,1)' ' WRITE(*,35)'Z,X,Y AXIS CONSTANTS OF MOLECULE:',(CONMON(I),I=1,3) WRITE(4,35)' OBSERVED RG-MOLECULE CONSTANTS:',(CONCMP(I),I=1,3) WRITE(4,35)'CALCULATED RG-MOLECULE CONSTANTS:',(CONCLC(I),I=1,3) WRITE(4,35)' Chi Squared = SUM(I.o-I.c)**2 :',CHISQ WRITE(4,1)' ' WRITE(4,35)'Z,X,Y AXIS CONSTANTS OF MOLECULE:',(CONMON(I),I=1,3) 35 FORMAT(1X,A,5X,3F12.3) WRITE(*,'(1x/1x,78(1H-)/)') WRITE(4,'(1x/1x,78(1H-)/)') c close(4) C STOP END C C_____________________________________________________________________________ C SUBROUTINE NORMAL(I,ANGLE) C C Normalisation of an angle (in radians) to range 0 to PI for I=1 and C 0 to 2*PI for I=2 C IMPLICIT REAL*8 (A-H,O-Z) DATA PI,TWOPI/3.141592654D0,6.283185307D0/ C IF(I.EQ.1)THEN A=DABS(ANGLE) IF(A.LE.PI)RETURN A=A-DINT(A/TWOPI)*TWOPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.-PI)ANGLE=ANGLE+TWOPI IF(ANGLE.GT.PI)ANGLE=ANGLE-TWOPI IF(ANGLE.LT.0.D0)ANGLE=-ANGLE RETURN ELSE IF(I.EQ.2)THEN A=DABS(ANGLE) A=A-DINT(A/TWOPI)*TWOPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.0.D0)ANGLE=ANGLE+TWOPI RETURN ENDIF C RETURN END C C_____________________________________________________________________________ C SUBROUTINE XYZAR(R,T,P,X,Y,Z) IMPLICIT REAL*8 (A-H,O-Z) CONV=57.29577951D0 C C...Cartesian coordinates of Ar from its polar coordinates C THETA=T/CONV PHI=P/CONV Z=R*DCOS(THETA) X=R*DSIN(THETA)*DCOS(PHI) Y=R*DSIN(THETA)*DSIN(PHI) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE INV IMPLICIT REAL*8(A-H,O-Z) COMMON /HBLK/A(3,3),S(3,3) c c...Inversion of a 3x3 general array c s(1,1)=a(2,2)*a(3,3)-a(3,2)*a(2,3) s(1,2)=a(3,2)*a(1,3)-a(1,2)*a(3,3) s(1,3)=a(1,2)*a(2,3)-a(2,2)*a(1,3) s(2,1)=a(2,3)*a(3,1)-a(2,1)*a(3,3) s(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3) s(2,3)=a(2,1)*a(1,3)-a(1,1)*a(2,3) s(3,1)=a(2,1)*a(3,2)-a(2,2)*a(3,1) s(3,2)=a(3,1)*a(1,2)-a(1,1)*a(3,2) s(3,3)=a(1,1)*a(2,2)-a(2,1)*a(1,2) c det=s(1,1)*a(1,1)+s(1,2)*a(2,1)+s(1,3)*a(3,1) do 1 i=1,3 do 1 j=1,3 a(i,j)=s(i,j)/det 1 continue c RETURN END C C_____________________________________________________________________________ C SUBROUTINE MOFI(R,THETA,PHI) IMPLICIT REAL*8(A-H,O-Z) COMMON /TEN/MU,IMON/HBLK/T(3,3),S(3,3) REAL*8 MU,IMON(3) C c ... inertial tensor of the complex c RSQ=R*R T(1,1)=IMON(1)+MU*RSQ*DSIN(THETA)**2 T(2,2)=IMON(2)+MU*RSQ*(DCOS(THETA)**2 . +(DSIN(THETA)*DSIN(PHI))**2) T(3,3)=IMON(3)+MU*RSQ*(DCOS(THETA)**2 . +(DSIN(THETA)*DCOS(PHI))**2) T(1,2)=-MU*RSQ*DCOS(THETA)*DSIN(THETA)*DCOS(PHI) T(2,1)=T(1,2) T(1,3)=-MU*RSQ*DCOS(THETA)*DSIN(THETA)*DSIN(PHI) T(3,1)=T(1,3) T(2,3)=-MU*RSQ*DSIN(THETA)**2*DSIN(PHI)*DCOS(PHI) T(3,2)=T(2,3) c c ... moments of inertia of the complex c CALL JACKS(3) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE JACKS(N) c c JACOBI type matrix diagonalisation routine for real, symmetric matrices c c H - array to be diagonalised (contains the eigenvalues on the diagonal c on output c S - contains the eigenvectors on output (in order of eigenvalues) c N - dimension of the array to be diagonalised contained in the top c left hand corner of H c c This version sorts eigenvalues and eigenvectors c IMPLICIT REAL*8 (A-H,O-Z) COMMON /HBLK/H(3,3),S(3,3)/SORTCC/WK(3),IPT(3) REAL*8 D(3) c c ... initialisation and setup of limits c ss=0.7071067811865475d0 n1=n-1 sum=0.d0 do 20 i=1,n1 s(i,i)=1.d0 k=i+1 do 20 m=k,n sum=sum+h(i,m)**2 s(i,m)=0.d0 s(m,i)=0.d0 20 continue s(n,n)=1.d0 unu=dsqrt(2.d0*sum) ene=n c c ... convergence criterion c unufi=1.d-12 c c c ... iteration sequence - determination of SIN and COS for rotation c 54 unu=unu/ene 55 ind=0 do 89 m=2,n l=m-1 do 89 i=1,l if(h(i,m))60,89,60 60 q=dabs(h(i,m)) if(q-unu)89,62,62 62 ind=1 if(h(i,i)-h(m,m))68,64,68 64 si=dsign(ss,h(i,m)) co=ss goto 77 68 x1=-h(i,m)*2.d0/(h(i,i)-h(m,m)) x1a=dabs(x1) if(1.d30-x1a)66,66,71 66 si=dsign(ss,x1) co=ss goto 77 71 if(1.d-30-x1a)72,89,89 72 t=datan(x1)*0.5d0 si=dsin(t) co=dcos(t) c c ... transformation loop c 77 do 79 j=1,n hji=h(j,i) sji=s(j,i) h(j,i)=hji*co-h(j,m)*si h(j,m)=hji*si+h(j,m)*co s(j,i)=sji*co-s(j,m)*si s(j,m)=sji*si+s(j,m)*co if(j.eq.i.or.j.eq.m)goto 79 h(i,j)=h(j,i) h(m,j)=h(j,m) 79 continue c h(i,i)=co*h(i,i)-si*h(m,i) h(m,m)=co*h(m,m)+si*h(i,m) h(i,m)=0.d0 h(m,i)=0.d0 89 continue if(ind)92,91,55 91 if(unu-unufi)92,92,54 c c ... sorting of eigenvalues and eigenvectors c 92 DO 115 I=1,N D(I)=H(I,I) DO 115 J=1,N H(I,J)=S(I,J) 115 CONTINUE C DO 122 I=1,N IPT(I)=I WK(I)=D(I) 122 CONTINUE C CALL SORTC(1,N) C DO 117 II=1,N I=IPT(II) DO 117 J=1,N S(J,II)=H(J,I) 117 CONTINUE DO 118 II=1,N I=IPT(II) H(II,II)=D(I) 118 CONTINUE C return end C C____________________________________________________________________________ C SUBROUTINE SORTC(N,M) IMPLICIT REAL*8(A-H,O-Z) COMMON /SORTCC/WK(3),IPT(3) C C ... This routine sorts the quantities in vector WK in ascending order C of magnitude and also accordingly rearranges vector IPT of pointers C to original positions of sorted quantities C DO 101 I=N,M-1 J=I 106 J=J+1 IF(WK(J)-WK(I))103,104,104 103 EE=WK(I) WK(I)=WK(J) WK(J)=EE K=IPT(I) IPT(I)=IPT(J) IPT(J)=K 104 IF(J.EQ.M)GOTO 101 GOTO 106 101 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE FDATA(X,TERMS,NTERMS) C C User supplied routine for general SVD fit. This sets up the linear eq. C for a given value of X (for point X(i),Y(i)) ie. determines the NTERMS C coefficients of the NTERMS parameters to be fitted and places them in C TERMS. C In the present case X is just the equation number and the coefficients C are only to be transferred from the appropriate row of COEFF read in C from the data file. C COMMON /BLKCF/DERSVD(3,3) DIMENSION TERMS(NTERMS) C N=NINT(X) DO 1 I=1,NTERMS TERMS(I)=DERSVD(N,I) 1 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SVDFIT(X,Y,SIG,NDATA,A,MA,U,V,W,MP,NP,CHISQ,FUNCS) C C Given a set of NDATA points X(I),Y(I) with individual standard deviations C SIG(I) use CHISQ minimisation to determine the MA coefficients A of the C fitting function \-- C y = > A . AFUNC(x) C /-- C Here we solve the fitting equations using singular value decomposition C of the NDATA by MA matrix. Arrays U,W,V provide workspace on input, on C output they define the singular value decomposition, and can be used to C obtain the covariance matrix. MP,NP are the physical dimensions of the C matrices U,V,W, as indicated below. It is necessary that MP.ge.NDATA, C NP.ge.MA. The program returns values for the MA parameters A, and C CHISQ. The user supplies a subroutine FUNCS(X,AFUNC,MA) that returns C the MA basis functions evaluated at x=X in the array AFUNC. C C C...NMAX,MMAX - maximum expected NDATA and MA C PARAMETER(NMAX=200,MMAX=20,TOL=1.E-5) DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),V(NP,NP), * U(MP,NP),W(NP),B(NMAX),AFUNC(MMAX) C C...Accumulate coefficients of the fitting matrix C DO 12 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) TMP=1./SIG(I) DO 11 J=1,MA U(I,J)=AFUNC(J)*TMP 11 CONTINUE B(I)=Y(I)*TMP 12 CONTINUE C C...Singular value decomposition C CALL SVDCMP(U,NDATA,MA,MP,NP,W,V) C C...Edit the singular values, given TOL from the parameter statement, C between here ..... C WMAX=0. DO 13 J=1,MA IF(W(J).GT.WMAX)WMAX=W(J) 13 CONTINUE THRESH=TOL*WMAX DO 14 J=1,MA IF(W(J).LT.THRESH)W(J)=0. 14 CONTINUE C C...and here ......... C CALL SVBKSB(U,W,V,NDATA,MA,MP,NP,B,A) C C...Evaulate chi-square C CHISQ=0. DO 16 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) SUM=0. DO 15 J=1,MA SUM=SUM+A(J)*AFUNC(J) 15 CONTINUE CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2 16 CONTINUE RETURN END C C_____________________________________________________________________________ C SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V) C C Given a matrix A with logical dimensions M by N and physical dimensions C MP by NP, this routine computes its singular value decomposition: C C A = U W V' C C The matrix U replaces A on output. The diagonal matrix of singular C values W is output as vector W. The matrix V (not the transpose V') is C output as V. M must be greater or equal to N; if it is smaller then A C should be filled up to square with zero rows. C C C...NMAX = maximum anticipated value of N C PARAMETER (NMAX=200) DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX) C C...Housholder reduction to bidiagonal form C G=0.0 SCALE=0.0 ANORM=0.0 DO 25 I=1,N L=I+1 RV1(I)=SCALE*G G=0.0 S=0.0 SCALE=0.0 IF (I.LE.M) THEN DO 11 K=I,M SCALE=SCALE+ABS(A(K,I)) 11 CONTINUE IF (SCALE.NE.0.0) THEN DO 12 K=I,M A(K,I)=A(K,I)/SCALE S=S+A(K,I)*A(K,I) 12 CONTINUE F=A(I,I) G=-SIGN(SQRT(S),F) H=F*G-S A(I,I)=F-G IF (I.NE.N) THEN DO 15 J=L,N S=0.0 DO 13 K=I,M S=S+A(K,I)*A(K,J) 13 CONTINUE F=S/H DO 14 K=I,M A(K,J)=A(K,J)+F*A(K,I) 14 CONTINUE 15 CONTINUE ENDIF DO 16 K= I,M A(K,I)=SCALE*A(K,I) 16 CONTINUE ENDIF ENDIF W(I)=SCALE *G G=0.0 S=0.0 SCALE=0.0 IF ((I.LE.M).AND.(I.NE.N)) THEN DO 17 K=L,N SCALE=SCALE+ABS(A(I,K)) 17 CONTINUE IF (SCALE.NE.0.0) THEN DO 18 K=L,N A(I,K)=A(I,K)/SCALE S=S+A(I,K)*A(I,K) 18 CONTINUE F=A(I,L) G=-SIGN(SQRT(S),F) H=F*G-S A(I,L)=F-G DO 19 K=L,N RV1(K)=A(I,K)/H 19 CONTINUE IF (I.NE.M) THEN DO 23 J=L,M S=0.0 DO 21 K=L,N S=S+A(J,K)*A(I,K) 21 CONTINUE DO 22 K=L,N A(J,K)=A(J,K)+S*RV1(K) 22 CONTINUE 23 CONTINUE ENDIF DO 24 K=L,N A(I,K)=SCALE*A(I,K) 24 CONTINUE ENDIF ENDIF ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I)))) 25 CONTINUE C C...Accumulation of right-hand transformations (double division in loop over C J to avoid possible underflow C DO 32 I=N,1,-1 IF (I.LT.N) THEN IF (G.NE.0.0) THEN DO 26 J=L,N V(J,I)=(A(I,J)/A(I,L))/G 26 CONTINUE DO 29 J=L,N S=0.0 DO 27 K=L,N S=S+A(I,K)*V(K,J) 27 CONTINUE DO 28 K=L,N V(K,J)=V(K,J)+S*V(K,I) 28 CONTINUE 29 CONTINUE ENDIF DO 31 J=L,N V(I,J)=0.0 V(J,I)=0.0 31 CONTINUE ENDIF V(I,I)=1.0 G=RV1(I) L=I 32 CONTINUE C C...Accumulation of left-hand trnsformations C DO 39 I=N,1,-1 L=I+1 G=W(I) IF (I.LT.N) THEN DO 33 J=L,N A(I,J)=0.0 33 CONTINUE ENDIF IF (G.NE.0.0) THEN G=1.0/G IF (I.NE.N) THEN DO 36 J=L,N S=0.0 DO 34 K=L,M S=S+A(K,I)*A(K,J) 34 CONTINUE F=(S/A(I,I))*G DO 35 K=I,M A(K,J)=A(K,J)+F*A(K,I) 35 CONTINUE 36 CONTINUE ENDIF DO 37 J=I,M A(J,I)=A(J,I)*G 37 CONTINUE ELSE DO 38 J= I,M A(J,I)=0.0 38 CONTINUE ENDIF A(I,I)=A(I,I)+1.0 39 CONTINUE C C...Diagonalisation of the bidiagonal form (K loops over singular values, C ITS loops over allowed iterations, L loop tests for splitting - RV1(1) C is always zero) C DO 49 K=N,1,-1 DO 48 ITS=1,30 DO 41 L=K,1,-1 NM=L-1 IF ((ABS(RV1(L))+ANORM).EQ.ANORM) GO TO 2 IF ((ABS(W(NM))+ANORM).EQ.ANORM) GO TO 1 41 CONTINUE C C...cancellation of RV1(L) if L>1 1 C=0.0 S=1.0 DO 43 I=L,K F=S*RV1(I) IF ((ABS(F)+ANORM).NE.ANORM) THEN G=W(I) H=SQRT(F*F+G*G) W(I)=H H=1.0/H C= (G*H) S=-(F*H) DO 42 J=1,M Y=A(J,NM) Z=A(J,I) A(J,NM)=(Y*C)+(Z*S) A(J,I)=-(Y*S)+(Z*C) 42 CONTINUE ENDIF 43 CONTINUE 2 Z=W(K) C C...if L=K then convergence and singular value made nonnegative IF (L.EQ.K) THEN IF (Z.LT.0.0) THEN W(K)=-Z DO 44 J=1,N V(J,K)=-V(J,K) 44 CONTINUE ENDIF GO TO 3 ENDIF IF (ITS.EQ.30) PAUSE 'No convergence in 30 iterations' C C...shift from bottom 2x2 minor X=W(L) NM=K-1 Y=W(NM) G=RV1(NM) H=RV1(K) F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y) G=SQRT(F*F+1.0) F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X C C...next QR transformation C=1.0 S=1.0 DO 47 J=L,NM I=J+1 G=RV1(I) Y=W(I) H=S*G G=C*G Z=SQRT(F*F+H*H) RV1(J)=Z C=F/Z S=H/Z F= (X*C)+(G*S) G=-(X*S)+(G*C) H=Y*S Y=Y*C DO 45 NM=1,N X=V(NM,J) Z=V(NM,I) V(NM,J)= (X*C)+(Z*S) V(NM,I)=-(X*S)+(Z*C) 45 CONTINUE Z=SQRT(F*F+H*H) W(J)=Z C C...rotation can be arbitrary if Z=0 IF (Z.NE.0.0) THEN Z=1.0/Z C=F*Z S=H*Z ENDIF F= (C*G)+(S*Y) X=-(S*G)+(C*Y) DO 46 NM=1,M Y=A(NM,J) Z=A(NM,I) A(NM,J)= (Y*C)+(Z*S) A(NM,I)=-(Y*S)+(Z*C) 46 CONTINUE 47 CONTINUE RV1(L)=0.0 RV1(K)=F W(K)=X 48 CONTINUE 3 CONTINUE 49 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X) C C Solves A X = B for a vector X, where A is specified by the arrys U,W,V C as returned by SVDCMP. M and N are the logical dimensions of A, and will C be equal for square matrices. MP and NP are the physical dimensions C of A. B is the outpur right-hand side. X is the output solution vector. C No input quantities are destroyed, so the routine may be called C sequentially with different B's C C C...NMAX is the maximum anticipated value of N C PARAMETER (NMAX=200) DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(NMAX) C C...calculate U' B DO 12 J=1,N S=0. IF(W(J).NE.0.)THEN DO 11 I=1,M S=S+U(I,J)*B(I) 11 CONTINUE C C...divide by W.j S=S/W(J) ENDIF TMP(J)=S 12 CONTINUE C C...matrix multiply by V to get the answer X DO 14 J=1,N S=0. DO 13 JJ=1,N S=S+V(J,JJ)*TMP(JJ) 13 CONTINUE X(J)=S 14 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SVDVAR(V,MA,NP,W,CVM,NCVM) C C To evaluate the covariance matrix CVM of the fit for MA parameters C obtained SVDFIT, call this routine with matrices V,W as returned from C SVDFIT. NP,NCVM give the physical dimensions of V,W,CVM as indicated C below. C C C...NMAX is the maximum number of fit parameters C PARAMETER (MMAX=20) DIMENSION V(NP,NP),W(NP),CVM(NCVM,NCVM),WTI(MMAX) DO 11 I=1,MA WTI(I)=0. IF(W(I).NE.0.) WTI(I)=1./(W(I)*W(I)) 11 CONTINUE C C...Sum contributions to covariance matrix C DO 14 I=1,MA DO 13 J=1,I SUM=0. DO 12 K=1,MA SUM=SUM+V(I,K)*V(J,K)*WTI(K) 12 CONTINUE CVM(I,J)=SUM CVM(J,I)=SUM 13 CONTINUE 14 CONTINUE RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________