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 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 crc 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 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 eta = z00100000 nit = 0 small=25.0*10.e-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 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 if (np.gt.npq) go to 230 n0=0 do 220 n=np,npq nn=n+np-1 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 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 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 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 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 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)) * write (1,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