+++ /dev/null
- subroutine fitsq(rms,x,y,nn,t,b,non_conv)
- implicit real*8 (a-h,o-z)
- include 'COMMON.IOUNITS'
-c x and y are the vectors of coordinates (dimensioned (3,n)) of the two
-c structures to be superimposed. nn is 3*n, where n is the number of
-c points. t and b are respectively the translation vector and the
-c rotation matrix that transforms the second set of coordinates to the
-c frame of the first set.
-c eta = machine-specific variable
-
- dimension x(3*nn),y(3*nn),t(3)
- dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3)
- logical non_conv
-c eta = z00100000
-c small=25.0*rmdcon(3)
-c small=25.0*eta
-c small=25.0*10.e-10
-c the following is a very lenient value for 'small'
- small = 0.0001D0
- non_conv=.false.
- fn=nn
- do 10 i=1,3
- xav(i)=0.0D0
- yav(i)=0.0D0
- do 10 j=1,3
- 10 b(j,i)=0.0D0
- nc=0
-c
- do 30 n=1,nn
- do 20 i=1,3
-c write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i)
- xav(i)=xav(i)+x(nc+i)/fn
- 20 yav(i)=yav(i)+y(nc+i)/fn
- 30 nc=nc+3
-c
- do i=1,3
- t(i)=yav(i)-xav(i)
- enddo
-
- rms=0.0d0
- do n=1,nn
- do i=1,3
- rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2
- enddo
- enddo
- rms=dabs(rms/fn)
-
-c write(iout,*)'xav = ',(xav(j),j=1,3)
-c write(iout,*)'yav = ',(yav(j),j=1,3)
-c write(iout,*)'t = ',(t(j),j=1,3)
-c write(iout,*)'rms=',rms
- if (rms.lt.small) return
-
-
- nc=0
- rms=0.0D0
- do 50 n=1,nn
- do 40 i=1,3
- rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
- do 40 j=1,3
- b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
- 40 c(j,i)=b(j,i)
- 50 nc=nc+3
- call sivade(b,q,r,d,non_conv)
- sn3=dsign(1.0d0,d)
- do 120 i=1,3
- do 120 j=1,3
- 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3)
- call mvvad(b,xav,yav,t)
- do 130 i=1,3
- do 130 j=1,3
- rms=rms+2.0*c(j,i)*b(j,i)
- 130 b(j,i)=-b(j,i)
- if (dabs(rms).gt.small) go to 140
-* write (6,301)
- return
- 140 if (rms.gt.0.0d0) go to 150
-c write (iout,303) rms
- rms=0.0d0
-* stop
-c 150 write (iout,302) dsqrt(rms)
- 150 continue
- return
- 301 format (5x,'rms deviation negligible')
- 302 format (5x,'rms deviation ',f14.6)
- 303 format (//,5x,'negative ms deviation - ',f14.6)
- end
-c
- subroutine sivade(x,q,r,dt,non_conv)
- implicit real*8(a-h,o-z)
-c computes q,e and r such that q(t)xr = diag(e)
- dimension x(3,3),q(3,3),r(3,3),e(3)
- dimension h(3,3),p(3,3),u(3,3),d(3)
- logical non_conv
-c eta = z00100000
-c write (2,*) "SIVADE"
- nit = 0
- small=25.0*10.d-10
-c small=25.0*eta
-c small=2.0*rmdcon(3)
- xnrm=0.0d0
- do 20 i=1,3
- do 10 j=1,3
- xnrm=xnrm+x(j,i)*x(j,i)
- u(j,i)=0.0d0
- r(j,i)=0.0d0
- 10 h(j,i)=0.0d0
- u(i,i)=1.0
- 20 r(i,i)=1.0
- xnrm=dsqrt(xnrm)
- do 110 n=1,2
- xmax=0.0d0
- do 30 j=n,3
- 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n))
- a=0.0d0
- do 40 j=n,3
- h(j,n)=x(j,n)/xmax
- 40 a=a+h(j,n)*h(j,n)
- a=dsqrt(a)
- den=a*(a+dabs(h(n,n)))
- d(n)=1.0/den
- h(n,n)=h(n,n)+dsign(a,h(n,n))
- do 70 i=n,3
- s=0.0d0
- do 50 j=n,3
- 50 s=s+h(j,n)*x(j,i)
- s=d(n)*s
- do 60 j=n,3
- 60 x(j,i)=x(j,i)-s*h(j,n)
- 70 continue
- if (n.gt.1) go to 110
- xmax=dmax1(dabs(x(1,2)),dabs(x(1,3)))
- h(2,3)=x(1,2)/xmax
- h(3,3)=x(1,3)/xmax
- a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
- den=a*(a+dabs(h(2,3)))
- d(3)=1.0/den
- h(2,3)=h(2,3)+sign(a,h(2,3))
- do 100 i=1,3
- s=0.0d0
- do 80 j=2,3
- 80 s=s+h(j,3)*x(i,j)
- s=d(3)*s
- do 90 j=2,3
- 90 x(i,j)=x(i,j)-s*h(j,3)
- 100 continue
- 110 continue
- do 130 i=1,3
- do 120 j=1,3
- 120 p(j,i)=-d(1)*h(j,1)*h(i,1)
- 130 p(i,i)=1.0+p(i,i)
- do 140 i=2,3
- do 140 j=2,3
- u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2)
- 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3)
- call mmmul(p,u,q)
- 150 np=1
- nq=1
- nit=nit+1
-c write (2,*) "nit",nit," e",(x(i,i),i=1,3)
- if (nit.gt.10000) then
- print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
- non_conv=.true.
- return
- endif
- if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160
- x(2,3)=0.0d0
- nq=nq+1
- 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180
- x(1,2)=0.0d0
- if (x(2,3).ne.0.0d0) go to 170
- nq=nq+1
- go to 180
- 170 np=np+1
- 180 if (nq.eq.3) go to 310
- npq=4-np-nq
-c write (2,*) "np",np," npq",npq
- if (np.gt.npq) go to 230
- n0=0
- do 220 n=np,npq
- nn=n+np-1
-c write (2,*) "nn",nn
- if (dabs(x(nn,nn)).gt.small*xnrm) go to 220
- x(nn,nn)=0.0d0
- if (x(nn,nn+1).eq.0.0d0) go to 220
- n0=n0+1
-c write (2,*) "nn",nn
- go to (190,210,220),nn
- 190 do 200 j=2,3
- 200 call givns(x,q,1,j)
- go to 220
- 210 call givns(x,q,2,3)
- 220 continue
-c write (2,*) "nn",nn," np",np," nq",nq," n0",n0
-c write (2,*) "x",(x(i,i),i=1,3)
- if (n0.ne.0) go to 150
- 230 nn=3-nq
- a=x(nn,nn)*x(nn,nn)
- if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn)
- b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1)
- c=x(nn,nn)*x(nn,nn+1)
- dd=0.5*(a-b)
- xn2=c*c
- rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd))
- y=x(np,np)*x(np,np)-rt
- z=x(np,np)*x(np,np+1)
- do 300 n=np,nn
-c write (2,*) "n",n," a",a," b",b," c",c," y",y," z",z
- if (dabs(y).lt.dabs(z)) go to 240
- t=z/y
- c=1.0/dsqrt(1.0d0+t*t)
- s=c*t
- go to 250
- 240 t=y/z
- s=1.0/dsqrt(1.0d0+t*t)
- c=s*t
- 250 do 260 j=1,3
- v=x(j,n)
- w=x(j,n+1)
- x(j,n)=c*v+s*w
- x(j,n+1)=-s*v+c*w
- a=r(j,n)
- b=r(j,n+1)
- r(j,n)=c*a+s*b
- 260 r(j,n+1)=-s*a+c*b
- y=x(n,n)
- z=x(n+1,n)
- if (dabs(y).lt.dabs(z)) go to 270
- t=z/y
- c=1.0/dsqrt(1.0+t*t)
- s=c*t
- go to 280
- 270 t=y/z
- s=1.0/dsqrt(1.0+t*t)
- c=s*t
- 280 do 290 j=1,3
- v=x(n,j)
- w=x(n+1,j)
- a=q(j,n)
- b=q(j,n+1)
- x(n,j)=c*v+s*w
- x(n+1,j)=-s*v+c*w
- q(j,n)=c*a+s*b
- 290 q(j,n+1)=-s*a+c*b
- if (n.ge.nn) go to 300
- y=x(n,n+1)
- z=x(n,n+2)
- 300 continue
- go to 150
- 310 do 320 i=1,3
- 320 e(i)=x(i,i)
- nit=0
- 330 n0=0
- nit=nit+1
- if (nit.gt.10000) then
- print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
- non_conv=.true.
- return
- endif
-c write (2,*) "e",(e(i),i=1,3)
- do 360 i=1,3
- if (e(i).ge.0.0d0) go to 350
- e(i)=-e(i)
- do 340 j=1,3
- 340 q(j,i)=-q(j,i)
- 350 if (i.eq.1) go to 360
- if (dabs(e(i)).lt.dabs(e(i-1))) go to 360
- call switch(i,1,q,r,e)
- n0=n0+1
- 360 continue
- if (n0.ne.0) go to 330
-c write (2,*) "e",(e(i),i=1,3)
- if (dabs(e(3)).gt.small*xnrm) go to 370
- e(3)=0.0d0
- if (dabs(e(2)).gt.small*xnrm) go to 370
- e(2)=0.0d0
- 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3))
-c write (2,*) "nit",nit
-c write (2,501) (e(i),i=1,3)
- return
- 501 format (/,5x,'singular values - ',3e15.5)
- end
- subroutine givns(a,b,m,n)
- implicit real*8 (a-h,o-z)
- dimension a(3,3),b(3,3)
- if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10
- t=a(n,n)/a(m,n)
- s=1.0/dsqrt(1.0+t*t)
- c=s*t
- go to 20
- 10 t=a(m,n)/a(n,n)
- c=1.0/dsqrt(1.0+t*t)
- s=c*t
- 20 do 30 j=1,3
- v=a(m,j)
- w=a(n,j)
- x=b(j,m)
- y=b(j,n)
- a(m,j)=c*v-s*w
- a(n,j)=s*v+c*w
- b(j,m)=c*x-s*y
- 30 b(j,n)=s*x+c*y
- return
- end
- subroutine switch(n,m,u,v,d)
- implicit real*8 (a-h,o-z)
- dimension u(3,3),v(3,3),d(3)
- do 10 i=1,3
- tem=u(i,n)
- u(i,n)=u(i,n-1)
- u(i,n-1)=tem
- if (m.eq.0) go to 10
- tem=v(i,n)
- v(i,n)=v(i,n-1)
- v(i,n-1)=tem
- 10 continue
- tem=d(n)
- d(n)=d(n-1)
- d(n-1)=tem
- return
- end
- subroutine mvvad(b,xav,yav,t)
- implicit real*8 (a-h,o-z)
- dimension b(3,3),xav(3),yav(3),t(3)
-c dimension a(3,3),b(3),c(3),d(3)
-c do 10 j=1,3
-c d(j)=c(j)
-c do 10 i=1,3
-c 10 d(j)=d(j)+a(j,i)*b(i)
- do 10 j=1,3
- t(j)=yav(j)
- do 10 i=1,3
- 10 t(j)=t(j)+b(j,i)*xav(i)
- return
- end
- double precision function det (a,b,c)
- implicit real*8 (a-h,o-z)
- dimension a(3),b(3),c(3)
- det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3))
- 1 +a(3)*(b(1)*c(2)-b(2)*c(1))
- return
- end
- subroutine mmmul(a,b,c)
- implicit real*8 (a-h,o-z)
- dimension a(3,3),b(3,3),c(3,3)
- do 10 i=1,3
- do 10 j=1,3
- c(i,j)=0.0d0
- do 10 k=1,3
- 10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
- return
- end
- subroutine matvec(uvec,tmat,pvec,nback)
- implicit real*8 (a-h,o-z)
- real*8 tmat(3,3),uvec(3,nback), pvec(3,nback)
-c
- do 2 j=1,nback
- do 1 i=1,3
- uvec(i,j) = 0.0d0
- do 1 k=1,3
- 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)
- 2 continue
- return
- end
projects / unres.git / blobdiff