c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c program add3.f April 23, 1993
c
c Sept.,1995 declared plev,df,xf as real*4
c June,1995 input/output changed
c
c this version uses LAPACK for matrix inversion
c
c Combines three rotations independently estimated from datasets with
c unequal kappas.
c
c Critical value for the confidence region construction is calculated
c based on estimates developed by S. Johansen
c
c
C COMPILE WITH NUMLIB.F , numlib2.f and dtrans.f
C
c add three rotations and calculate the covariance matrix of the
c product rotation
c dhat=chat.bhat.ahat
c
c cov(tdhat)~(ba)**t*cov(tchat)*(ba) + (a**t)*cov(tbhat)*a + cov(tahat)
C
c
c
c
program add3
implicit double precision (a-h,o-z)
parameter (ndatmx=500,ndtmx2=2*ndatmx)
parameter (msect = 100,msig=2*msect*msect*msect+7*msect+6,msig2=
& 2*msect+3,mwork=2*msig2)
character filnam*70,name(3)*15,dum*1,datdum*30,datnam(3)*15
dimension a(3,3),cova(3,3),b(3,3),covb(3,3),c(3,3),covc(3,3),
& w(3,3),xa(3),xb(3),xc(3),ba(3,3),bat(3,3)
dimension axis(3),h(3),etai(msect,3,2),comp(3,3),
& qhati(4),eta(msect,3,3), num(2),wt(2,msect,20)
dimension u(2,msect,20,3),x2(ndatmx,ndatmx),xd(3),at(3,3),
& x(ndatmx,ndatmx),xt(ndatmx,ndatmx),npts(2,msect),axm(3,3),
& tempo(ndatmx,ndatmx),xtsx(ndatmx,ndatmx),covd(3,3)
dimension xtsxinv(ndatmx,ndatmx),pone(ndatmx,ndatmx),
& p11(ndatmx,ndatmx),p22(ndatmx,ndatmx)
dimension p11sq(ndatmx,ndatmx),p22sq(ndatmx,ndatmx),
& sigma(2,msect,3,3),p33(ndatmx,ndatmx),p33sq(ndatmx,ndatmx)
REAL*4 PLEV,df,XF
external R1
Common nsect,sigma,qhati,eta,etai
DATA NDAT/0/,NSECT/0/,rfact/4.0528473E07/,rfact2/6366.1977/,
& nsig/4/,maxfn/1000/,fmin/5000./,fmax/0./
1010 format(A)
1012 format(A29,A15)
1013 format(1x,A29,A15)
c
itour=1
write(*,*)'to estimate the composition A3*A2*A1'
1 write (*,1000) itour
1000 format ('Enter name of file for fit of rotation A',i1,': ',$)
read (*,'(a)') name(itour)
open (unit=10,file=name(itour),status='old',err=1)
write (*,'(a)') name(itour)
c read two first comment lines
read (10,1012) datdum,datnam(itour)
write (*,1013) datdum,datnam(itour)
read (10,'(a)') filnam
write (*,'(a)') filnam
c read alat,alon,angle
read(10,'(a)') filnam
write (*,'(a)') filnam
read(10,*) xc
write (*,'(3x,3f8.2)') xc
c read conf. level
read(10,'(a)') filnam
write (*,'(a)') filnam
read (10,*) plev,d1,d2
write(*,'(3f6.2)') plev,d1,d2
c read kappa estimated, degrees of freedom
read(10,'(a)') filnam
write (*,'(a)') filnam
read(10,*) hatkapc,dfc,xchi
write(*,*) hatkapc,dfc,xchi
if (dfc .lt. fmin) fmin=dfc
if (dfc .gt. fmax) fmax=dfc
c read number of points, sections ...
read(10,'(a)') dum
read(10,*) nrowx1,LSc
c read rotation matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
do 3 i=1,3
read(10,*) (c(i,j),j=1,3)
3 write (*,*) (c(i,j),j=1,3)
c read covariance matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
do 5 i=1,3
read (10,*) (covc(i,j),j=1,3)
5 write (*,*) (covc(i,j),j=1,3)
write(*,*) ' '
close(10)
c
c
write(*,1001)
1001 format(' enter data file name for this rotation ',$)
read(*,'(a)')filnam
OPEN(10,FILE=FILNAM,STATUS='UNKNOWN')
REWIND 10
READ(10,*) nsect
NDAT=0
num(iside)=0
NUM(1)=0
NUM(2)=0
do 51 I1 = 1,2
do 51 I2 = 1,nsect
npts(i1,i2) = 0
do 51 i3 = 1,3
do 51 i4 = 1,3
51 sigma(i1,i2,i3,i4) = 0.
do 200 ktimes = 1,nrowx1
100 READ(10,*,err=100) ISIDE,ISECT,ALAT,ALONG,SD
if(iside.gt.2) goto 100
npts(iside,isect) = npts(iside,isect) + 1
CALL TRANS1(ALAT,ALONG,AXIS)
n = npts(iside,isect)
wt(iside,isect,n) = 1.0/sd
do 105 i = 1,3
105 u(iside,isect,n ,i) = axis(i)
NUM(ISIDE)=Num(iside) + 1
do 110 j1 = 1,3
do 110 j2=1,3
110 sigma(iside,isect,j1,j2)=sigma(iside,isect,j1,j2)
& + axis(j1)*axis(j2)/(sd**2)
200 continue
201 close (10)
c
c setting up the eta vectors
c
250 alat = xc(1)
along = xc(2)
rho = xc(3)
CALL TRANS3(ALAT,ALONG,RHO,QHATI)
DO 207 I=1,3
207 H(I)=0.
RMIN=R1(H)*RFACT
c
C
C ETA IS THE MATRIX M(ETA) OF PAPER; ETAI IS THE MATRIX O-SUB I.
C ETA AND ETAI ARE SET BY R1.
C
c
c construct estimate of separate fit design matrix X, beginning
c with first 3 columns.
c
nrowx = nrowx1
ncolx = 3 + 2*nsect
do 505 k = 1, nrowx
do 505 j = 1, ncolx
505 x2(k,j) = 0.
nrow = 1
do 610 isect = 1, nsect
do 570 k = 1,3
do 570 j = 1,3
570 axm(k,j) = 0.
do 580 k = 1,3
do 580 j = 1,3
do 580 m = 1,3
580 axm(k,j) = axm(k,j) + c(k,m)*eta(isect,j,m)
do 600 np = 1,npts(2,isect)
do 595 k = 1,3
do 590 j = 1,3
590 x2(nrow,k) = x2(nrow,k) + u(2,isect,np,j)*axm(j,k)
595 x2(nrow,k) = x2(nrow,k)*wt(2,isect,np)*rfact2
600 nrow = nrow + 1
610 continue
nr = 1
nc = 4
do 650 isect = 1,nsect
do 620 k = 1,3
do 620 j = 1,2
620 axm(k,j) = 0.
do 630 k = 1,3
do 630 j = 1,2
do 630 m = 1,3
630 axm(k,j) = axm(k,j) + c(k,m)*etai(isect,m,j)
do 645 np = 1,npts(2,isect)
do 640 k = 1,3
x2(nr,nc) = x2(nr,nc) + u(2,isect,np,k)*axm(k,1)
640 x2(nr,nc+1) = x2(nr,nc+1) + u(2,isect,np,k)*axm(k,2)
x2(nr,nc) = wt(2,isect,np)*x2(nr,nc)*rfact2
x2(nr,nc+1) = wt(2,isect,np)*x2(nr,nc+1)*rfact2
645 nr = nr + 1
650 nc = nc + 2
nc = 4
do 680 isect = 1,nsect
do 670 np = 1,npts(1,isect)
do 660 k = 1,3
x2(nr,nc) = x2(nr,nc) + u(1,isect,np,k)*etai(isect,k,1)
660 x2(nr,nc+1) = x2(nr,nc+1) +
& u(1,isect,np,k)*etai(isect,k,2)
x2(nr,nc) = x2(nr,nc)*wt(1,isect,np)*rfact2
x2(nr,nc+1) = x2(nr,nc+1)*wt(1,isect,np)*rfact2
670 nr = nr + 1
680 nc = nc + 2
do 690 i = 1,nrowx
do 690 j = 1,ncolx
690 x2(i,j) = x2(i,j)*hatkapc**0.5
close(10)
c
itour=itour+1
goto (20,35,160) itour-1
20 do 30 i=1,3
xa(i)=xc(i)
do 30 j=1,3
a(i,j)=c(i,j)
30 cova(i,j)=covc(i,j)
hatkapa=hatkapc
dfa=dfc
lsa = nsect
nrowxa = nrowx
ncolxa = ncolx
do 31 i = 1,nrowxa
do 31 j = 1,ncolxa
31 x(i,j) = x2(i,j)
goto 1
c
c Construct estimated design matrix XK from the separate fit matrices
c
c
35 do 42 i = 1,3
xb(i) = xc(i)
do 42 j = 1,3
b(i,j) = c(i,j)
42 covb(i,j) = covc(i,j)
hatkapb = hatkapc
dfb = dfc
lsb = nsect
nrowxb = nrowx
ncolxb = ncolx
nrab = nrowxa + nrowxb
ncab = ncolxa + ncolxb
do 45 i = 1,nrowxb
do 45 j = 1,ncolxb
45 x(nrowxa + i,ncolxa+j) = x2(i,j)
do 150 i = 1,nrowxa
do 150 j = ncolxa + 1, ncab
150 x(i,j) = 0.0
do 155 i = nrowxa+1,nrab
do 155 j = 1,ncolxa
155 x(i,j) = 0.0
go to 1
c
c
160 call mul(b,a,ba)
nsect = nsect + lsa+lsb
nrowxc = nrowx
ncolxc = ncolx
nrowx = nrab + nrowxc
if (nrowx .gt. ndatmx) then
write(*,*)'parameter ndatmx is exceeded. '
write(*,*)'there are ',nrowx,' data points.'
stop
endif
ncolx = ncab + ncolxc - 3
if (ncolx .gt. 2*msect+3) then
write(*,*)'STOP! msect exceeded.'
stop
endif
do 170 i = 1,nrowxc
do 170 j = 4,ncolxc
170 x(nrab + i, ncab+j-3) = x2(i,j)
do 175 i = 1,nrab
do 175 j = ncab+1,ncolx
175 x(i,j) = 0.0
do 180 i = nrab+1,nrowx
do 185 j = 4,ncolxa
185 x(i,j) = 0.0
do 180 j = ncolxa + 4, ncab
180 x(i,j) = 0.0
do 190 i = 1,nrowxc
do 190 j = 1,3
x(nrab+i,j) = 0.0
do 195 k = 1,3
195 x(nrab+i,j) = x(nrab+i,j)-x2(i,k)*ba(k,j)
190 continue
do 800 i = 1,nrowxc
do 800 j = 1,3
x(nrab+i,ncolxa+j) = 0.0
do 810 k = 1,3
810 x(nrab+i,ncolxa+j) = x(nrab+i,ncolxa+j) - x2(i,k)*b(k,j)
800 continue
c
c open(13,file='x.out',status='unknown')
c do 370 i = 1,nrowx
c write(13,*) (x(i,j),j = 1,ncolx)
c 370 write(13,*) ' '
c close (13)
c
c transform X to S^(-1/2)*X
c
do 303 i = 1,nrowxa
do 303 j = 1,ncolx
303 x(i,j)=x(i,j)*hatkapa**.5
do 305 i = nrowxa + 1,nrowxa+nrowxb
do 305 j = 1,ncolx
305 x(i,j)=x(i,j)*hatkapb**.5
do 307 i = nrowxa+nrowxb+1,nrowx
do 307 j = 1,ncolx
307 x(i,j) = x(i,j)*hatkapc**.5
c
c
c
c calculate P-one
c
write(6,*)'calculating projection matrix'
do 400 i = 1,nrowx
do 400 j = 1,ncolx
400 xt(j,i) = x(i,j)
call matmul(xt,x,xtsx,ncolx,nrowx,ncolx)
do 401 i = 1,ncolx
do 401 j = 1,ncolx
401 xtsxinv(i,j)=xtsx(i,j)
c
c open(13,file='xtsx',status='unknown')
c do 405 i = 1,ncolx
c 405 write(13,*)(xtsx(i,j),j=1,ncolx)
c close(13)
c
call dpotrf('U',ncolx,xtsxinv,ndatmx,info)
if (info.ne.0) write(6,*)'failure in LAPACK routine. info = ',info
call dpotri('U',ncolx,xtsxinv,ndatmx,info)
if (info.ne.0) write(6,*)'failure in LAPACK routine. info = ',info
do 407 i = 2,ncolx
do 407 j = 1,i-1
407 xtsxinv(i,j)=xtsxinv(j,i)
c
c
c call matmul(xtsxinv, xtsx, tempo,ncolx,ncolx,ncolx)
c open(13,file='Lapackproduct',status='unknown')
c do 408 i = 1,ncolx
c write(13,*)(tempo(i,k),k=1,ncolx)
c 408 write(13,*)' '
c close(13)
c
call matmul(xtsxinv,xt,tempo,ncolx,ncolx,nrowx)
call matmul(x,tempo,pone,nrowx,ncolx,nrowx)
c
c
c Calculation of correction terms A and C
c
do 801 i = 1,nrowxa
do 801 j = 1,nrowxa
801 p11(i,j)=pone(i,j)
do 805 i = 1,nrowxb
do 805 j = 1,nrowxb
805 p22(i,j)=pone(i+nrowxa,j+nrowxa)
do 807 i = 1,nrowxc
do 807 j = 1,nrowxc
807 p33(i,j)=pone(i+nrab,j+nrab)
call matmul(p11,p11,p11sq,nrowxa,nrowxa,nrowxa)
call matmul(p22,p22,p22sq,nrowxb,nrowxb,nrowxb)
call matmul(p33,p33,p33sq,nrowxc,nrowxc,nrowxc)
call trace(p11,nrowxa,trp11)
call trace(p11sq,nrowxa,trp11sq)
call trace(p22,nrowxb,trp22)
call trace(p22sq,nrowxb,trp22sq)
call trace(p33,nrowxc,trp33)
call trace(p33sq,nrowxc,trp33sq)
c
capa=(ncolxa+trp11sq-2*trp11)/dfa + (ncolxb+trp22sq-2*trp22)/dfb
& + (ncolxc+trp33sq-2*trp33)/dfc
capc = (ncolxa-trp11)**2/dfa + (ncolxb-trp22)**2/dfb
& + (ncolxc-trp33)**2/dfc - capa
write(6,*)'A, C',capa,capc
c
c calculate c and f, where Q is cF(3,f)
c
exq = 3 + 2*capa
vq = 6 + 14*capa + 2*capc
write (6,*) 'mean of Q = ', exq
write (6,*) 'variance of Q = ', vq
ex2 = exq*exq
write (6,*)
df = (12*vq + 2*ex2)/(3*vq - 2*ex2)
cfact = exq*(1 - 2/df)
write (6,*) 'deg. freedom ', df
write (6,*) ' c = ',cfact
call xidf(plev,3.0,df,xf,ier)
cxf = cfact*xf
write(6,*)'xf, critical value = c*f(3,df):',xf,cxf
c
DO 36 I=1,3
DO 36 J=1,3
COVA(I,J)=COVA(I,J)/HATKAPA
COVB(I,J)=COVB(I,J)/HATKAPB
36 covc(i,j)=covc(i,j)/hatkapc
c
c adjust hatkapd for use in conreg.f
HATKAPd = 3./cfact
C
c calculate covariance matrix for the composition
call transp(ba,bat)
call mul(bat,covc,w)
call mul(w,ba,covd)
call transp(a,at)
call mul(at,covb,w)
call mul(w,a,covb)
do 40 i=1,3
do 40 j=1,3
40 covd(i,j)=covd(i,j)+covb(i,j)+cova(i,j)
c calculate the inverse of covariance matrix H11.2
call inv(covd,w)
c calculate comp=c.b.a
call mul(c,ba,comp)
call trans9(comp,h,rho,errchk)
call trans6(h,xd(1),xd(2))
xd(3)=rho*57.29577951
c
50 write (*,*) ' Product rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xd
write (*,*) 'Rotation matrix'
do 55 i=1,3
55 write (*,*) (comp(i,j),j=1,3)
write (*,*) 'Covariance matrix'
do 60 i=1,3
60 write (*,*) (covd(i,j),j=1,3)
c
65 write(*,1065)
1065 format(/,' Enter name of the output file: ',$)
read(*,'(a)') filnam
open(unit=15,file=filnam,status='unknown',err=65)
write(*,'(a)') filnam
write(15,1060) name(3),name(2),name(1)
1060 format (' ','Results from adding ',a,a,a,/)
write(15,*) 'Fitted rotation--alat,along,rho: '
write(15,*) xd
write(15,*)'confidence level'
write(15,*) plev,1.0,1.0
write(15,*) '3/cfact, degrees of freedom, crit.value, flag'
c flag = 1 is output for use by the program addplus.f
write(15,'(f12.6,2x,f9.4,2x,f12.4,2x,f3.1)') hatkapd,df,cxf,1.0
write(15,*) 'Number of points, fmin, fmax, rank'
write(15,*) nrowx,fmin,fmax,ncolx
write(15,*) 'ahat: '
do 70 i=1,3
70 write(15,*) (comp(i,j),j=1,3)
write(15,*) 'covariance matrix'
do 75 i=1,3
75 write(15,*) (covd(i,j),j=1,3)
write(15,*) 'H11.2 matrix: '
do 80 i=1,3
80 write(15,*) (w(i,j),j=1,3)
close (15)
c
stop
end
c****************************************************************
c
subroutine inv(a,b)
c invert 3x3 matrix
double precision a(3,3), b(3,3),det
det=a(1,1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))-
& a(1,2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))+
& a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
if (det.eq.0.) goto 5
b(1,1)=a(2,2)*a(3,3)-a(3,2)*a(2,3)
b(1,2)=a(3,2)*a(1,3)-a(1,2)*a(3,3)
b(1,3)=a(1,2)*a(2,3)-a(2,2)*a(1,3)
b(2,1)=a(3,1)*a(2,3)-a(2,1)*a(3,3)
b(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3)
b(2,3)=a(2,1)*a(1,3)-a(1,1)*a(2,3)
b(3,1)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
b(3,2)=a(3,1)*a(1,2)-a(1,1)*a(3,2)
b(3,3)=a(1,1)*a(2,2)-a(2,1)*a(1,2)
do 1 i=1,3
do 1 j=1,3
1 b(i,j)=b(i,j)/det
return
5 write (*,*) 'Trouble: det=0'
stop
end
c****************************************************************
c
subroutine mul(a,b,c)
c multiply 3x3 matrices: c=a.b
double precision a(3,3),b(3,3),c(3,3)
do 5 i=1,3
do 5 j=1,3
c(i,j)=0.
do 5 k=1,3
5 c(i,j)=c(i,j)+a(i,k)*b(k,j)
return
end
c****************************************************************
c
subroutine transp(a,b)
c calculate transpose of a: b=a**t
double precision a(3,3),b(3,3)
do 5 i=1,3
do 5 j=1,3
5 b(i,j)=a(j,i)
return
end
c****************************************************************
C
C
FUNCTION R1(H)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER (MSECT=100)
DIMENSION H(3),SIGMA(2,MSECT,3,3),QHATI(4),ETA(MSECT,3,3),
& ETAI(MSECT,3,2),AHAT(3,3),SIG(3,3),D(3),Z(3,3),WK(6),QHAT(4)
COMMON NSECT,SIGMA,QHATI,ETA,ETAI
C
CALL TRANS2(H,QHATI,QHAT)
CALL TRANS4(QHAT,AHAT)
R1=0.
DO 100 I=1,NSECT
DO 110 J=1,3
DO 110 K=1,3
SIG(J,K)=SIGMA(1,I,J,K)
DO 110 K1=1,3
DO 110 K2=1,3
110 SIG(J,K)=SIG(J,K)+AHAT(K1,J)*SIGMA(2,I,K1,K2)*AHAT(K2,K)
CALL JACOBI(SIG,3,3,D,Z,3,WK,NROT)
IF (NROT.LT.0) WRITE(6,*) 'SUBROUTINE JACOBI(3)--NROT: ',NROT
ETA(I,1,1)=0.
ETA(I,2,1)=Z(3,1)
ETA(I,3,1)=-Z(2,1)
ETA(I,1,2)=-Z(3,1)
ETA(I,2,2)=0.
ETA(I,3,2)=Z(1,1)
ETA(I,1,3)=Z(2,1)
ETA(I,2,3)=-Z(1,1)
ETA(I,3,3)=0.
IF (ABS(Z(3,1)).GT.(.2)) THEN
DO 120 J=1,3
ETAI(I,J,1)=ETA(I,J,1)
120 ETAI(I,J,2)=ETA(I,J,2)
ELSE IF (ABS(Z(1,1)).GT.(.2)) THEN
DO 125 J=1,3
ETAI(I,J,1)=ETA(I,J,2)
125 ETAI(I,J,2)=ETA(I,J,3)
ELSE
DO 130 J=1,3
ETAI(I,J,1)=ETA(I,J,1)
130 ETAI(I,J,2)=ETA(I,J,3)
ENDIF
100 R1=R1+D(1)
RETURN
END
C
c
c
subroutine matmul(a,b,c,nra,nrb,ncb)
c
c this subroutine only works for matrices dimensioned ndatmx by ndatmx
c
implicit double precision (a-h,o-z)
parameter (ndatmx = 500)
dimension a(ndatmx, ndatmx),b(ndatmx,ndatmx),c(ndatmx,ndatmx)
do 10 i = 1,nra
do 10 j = 1, ncb
10 c(i,j) = 0.
do 21 i = 1, nra
do 21 j = 1,ncb
do 20 k = 1, nrb
20 c(i,j) = c(i,j) + a(i,k)*b(k,j)
21 continue
return
end
c
c
c
subroutine trace(a,ndima,t)
c
implicit double precision (a-h,o-z)
parameter (ndatmx = 500)
dimension a(ndatmx,ndatmx)
t = 0.
do 10 k = 1,ndima
10 t = t + a(k,k)
return
end