c***************teclib.f************************************************
c23456789012345678901234567890123456789012345678901234567890123456789012
c December, 1994
c
c Sept., 1995 modified critcalc, declaring plev,xf,df as real*4
c May, 1995 corrections to trimat (dimensioning)
c October, 1999 dimension statement in vside rewritten to limit
c lines to 72 characters
c
c subroutines used to set up design matrix for regression equations
c for combining rotation estimates with unequal concentration
c parameters
c
c
c CAUTION: the parameters ndatmx and msect must be set with the same
c values in all subroutines and calling programs
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
function r2(h)
c calculation of hellinger criterion r for a triple junction.
implicit double precision(a-h,o-z)
parameter (msect=70)
dimension h(6),sigma(3,msect,3,3),qhati(12),eta(msect,3,3),
& etai(msect,3,2),ahat(3,3,2),sig(3,3),d(3),z(3,3),wk(6),qhat(8)
common nsect,sigma,qhati,eta,etai
c
call trans2(h,qhati,qhat)
call trans4(qhat,ahat(1,1,1))
call trans2(h(4),qhati(5),qhat(5))
call trans4(qhat(5),ahat(1,1,2))
r2=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 ks=2,3
ks1=ks-1
do 110 k1=1,3
do 110 k2=1,3
110 sig(j,k)=sig(j,k)+ahat(k1,j,ks1)*sigma(ks,i,k1,k2)*ahat(k2,k,ks1)
call jacobi(sig,3,3,d,z,3,wk,nrot)
if (nrot.lt.0) then
write(6,*) 'subroutine jacobi(3)--nrot: ',nrot
write(6,*) 'vector h: ',h
endif
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 r2=r2+d(1)
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine qmulti(q)
implicit double precision(a-h,o-z)
dimension q(12),qbar(4)
do 10 i=2,4
10 qbar(i)=-q(i)
qbar(1)=q(1)
call qmult(q(5),qbar,q(9))
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Construction of the design matrix X from a triple-junction fit
c
subroutine trimat(ndat,ncolx,ahat,num,nsect1,nfsect1,
& nsect2,nfsect2,nsect3,nfsect3,x)
c
implicit double precision(a-h,o-z)
parameter(msect=70,ndatmx=400,msig2=2*msect+6)
dimension ahat(3,3,3),num(6),x(ndatmx,ndatmx),
& y(ndatmx,msig2),b(3,3),sigma(3,msect,3,3),qhati(12),
& eta(msect,3,3),etai(msect,3,2),wt(3,msect,ndatmx),
& u(3,msect,ndatmx,3)
common nsect, sigma, qhati, eta, etai, wt, npts,u
c
c
do 225 i = 1,ndat
do 225 j = 1,ncolx
225 x(i,j) = 0.
do 230 i = 1,3
do 230 j = 1,3
230 b(i,j)=ahat(i,j,1)
call vside(num(1),b,nsect1,nfsect1,2,y)
do 275 i = 1,num(1)
do 270 j = 1,3
270 x(i,j) = y(i,j)
do 275 j = 4,2*nsect1+3
275 x(i,j+3) = y(i,j)
call uside(num(2),nsect1,nfsect1,1,y)
do 285 i = 1,num(2)
do 285 j = 1,2*nsect1
285 x(i+num(1),j+6)=y(i,j)
call vside(num(5),b,nsect3,nfsect3,2,y)
moverow = num(1)+num(2)+num(3)+num(4)
movecol=6+2*nsect1+2*nsect2
do 295 i = 1,num(5)
do 290 j = 1,3
290 x(i+moverow,j)=y(i,j)
do 295 j = 4,3+2*nsect3
295 x(i+moverow,j+movecol-3)=y(i,j)
do 300 i = 1,3
do 300 j = 1,3
300 b(i,j) = ahat(i,j,2)
call vside(num(3),b,nsect2,nfsect2,3,y)
moverow = num(1)+num(2)
movecol = 6+2*nsect1
do 330 i = 1,num(3)
do 325 j = 1,3
325 x(i+moverow,j+3) = y(i,j)
do 330 j = 4,2*nsect2 +3
330 x(i+moverow,j+movecol-3)=y(i,j)
call uside(num(4),nsect2,nfsect2,1,y)
moverow = moverow + num(3)
do 335 i = 1,num(4)
do 335 j = 1,2*nsect2
335 x(i+moverow,j+movecol)=y(i,j)
call vside(num(6),b,nsect3,nfsect3,3,y)
moverow = moverow + num(4) + num(5)
movecol=6+2*nsect1+2*nsect2
do 345 i = 1,num(6)
do 340 j = 1,3
340 x(i+moverow,j+3)=y(i,j)
do 345 j = 4,2*nsect3+3
345 x(i+moverow,j+movecol-3)=y(i,j)
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c23456789012345678901234567890123456789012345678901234567890123456789012
c
subroutine vside(nrowx,b,numsec,nfsec,iside,y)
c
implicit double precision(a-h,o-z)
parameter(msect=70,ndatmx=400,msig2=2*msect+6)
dimension b(3,3),eta(msect,3,3),u(3,msect,ndatmx,3),
& y(ndatmx,msig2),axm(3,3),npts(3,msect),wt(3,msect,ndatmx),
& etai(msect,3,2)
dimension sigma(3,msect,3,3),qhati(12)
common nsect,sigma,qhati,eta,etai,wt,npts,u
data rfact2/6366.1977/
c
do 505 k = 1,nrowx
do 505 j = 1,3+numsec*2
505 y(k,j)=0.
nrow = 1
do 610 isect = nfsec,nfsec + numsec-1
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) + b(k,m)*eta(isect,j,m)
do 600 np = 1,npts(iside,isect)
do 595 k = 1,3
do 590 j = 1,3
590 y(nrow,k) = y(nrow,k)+u(iside,isect,np,j)*axm(j,k)
595 y(nrow,k)= y(nrow,k)*wt(iside,isect,np)*rfact2
600 nrow = nrow + 1
610 continue
nr=1
nc=4
do 650 isect = nfsec,nfsec+numsec-1
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)+b(k,m)*etai(isect,m,j)
do 645 np = 1,npts(iside,isect)
do 640 k = 1,3
y(nr,nc) = y(nr,nc) + u(iside,isect,np,k)*axm(k,1)
640 y(nr,nc+1) = y(nr,nc+1) + u(iside,isect,np,k)*axm(k,2)
y(nr,nc) = wt(iside,isect,np)*y(nr,nc)*rfact2
y(nr,nc+1) = wt(iside,isect,np)*y(nr,nc+1)*rfact2
645 nr = nr+1
650 nc = nc + 2
c
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine uside(nrowx,numsec,nfsec,iside,y)
c
implicit double precision(a-h,o-z)
parameter(msect=70,ndatmx=400,msig2 = 2*msect+6)
dimension u(3,msect,ndatmx,3),y(ndatmx,msig2),
& npts(3,msect),wt(3,msect,ndatmx),etai(msect,3,2),
& qhati(12),eta(msect,3,3),sigma(3,msect,3,3)
common nsect,sigma,qhati,eta,etai,wt,npts,u
data rfact2/6366.1977/
c
do 600 i = 1,nrowx
do 600 j = 1,2*numsec
600 y(i,j) = 0.
nr=1
nc=1
do 680 isect = nfsec, nfsec+ numsec-1
do 670 np = 1,npts(iside,isect)
do 660 k = 1,3
y(nr,nc)=y(nr,nc)+u(iside,isect,np,k)*etai(isect,k,1)
660 y(nr,nc+1) = y(nr,nc+1) + u(iside,isect,np,k)*etai(isect,k,2)
y(nr,nc) = wt(iside,isect,np)*y(nr,nc)*rfact2
y(nr,nc+1) = wt(iside,isect,np)*y(nr,nc+1)*rfact2
670 nr = nr + 1
680 nc = nc + 2
c
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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 = 400)
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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine trace(a,ndima,t)
c
implicit double precision (a-h,o-z)
parameter (ndatmx = 400)
dimension a(ndatmx,ndatmx)
t = 0.
do 10 k = 1,ndima
10 t = t + a(k,k)
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine invert(x,ndim)
c
c inverts matrix dimensioned ndatmx by ndatmx by overwriting x
c uses LAPACK routines
c
implicit double precision (a-h,o-z)
parameter(ndatmx=400)
dimension x(ndatmx,ndatmx)
c
write(6,*) 'calculating a matrix inverse'
call dpotrf('U',ndim,x,ndatmx,info)
if (info .ne. 0) write(6,*) 'failure in LAPACK routine. info = ',
& info
call dpotri('U',ndim,x,ndatmx,info)
if (info.ne.0) write (6,*) 'failure in LAPACK routine. info = ',
& info
do 30 i = 2,ndim
do 30 j = 1, i-1
30 x(i,j)=x(j,i)
c
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine projection(x,nrowx,ncolx,p)
c
c matrices must be dimensioned ndatmx by ndatmx in the calling program
c
implicit double precision(a-h,o-z)
parameter (ndatmx=400)
dimension x(ndatmx,ndatmx),xt(ndatmx,ndatmx),xtx(ndatmx,ndatmx),
& tempo(ndatmx,ndatmx),p(ndatmx,ndatmx)
c
c
write(6,*) 'calculating a projection matrix'
do 20 i = 1,nrowx
do 20 j = 1,ncolx
20 xt(j,i)=x(i,j)
call matmul(xt,x,xtx,ncolx,nrowx,ncolx)
c
c xtx is replaced with its inverse
call invert(xtx,ncolx)
call matmul(x,xtx,tempo,nrowx,ncolx,ncolx)
do 40 i = 1,ncolx
do 40 j =1,nrowx
40 xt(i,j)=x(j,i)
call matmul(tempo,xt,p,nrowx,ncolx,nrowx)
c
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine abc(p0,p1,a,b,c)
c
c calculates the correction terms A,B,C of Johansens'paper
c
implicit double precision (a-h,o-z)
parameter (ndatmx=400,msect=70)
dimension p0(ndatmx,ndatmx),p1(ndatmx,ndatmx),p01(ndatmx,ndatmx),
& p02(ndatmx,ndatmx),p03(ndatmx,ndatmx),p11(ndatmx,ndatmx),
& p12(ndatmx,ndatmx),p13(ndatmx,ndatmx),sigma(3,msect,3,3)
dimension qhati(12),eta(msect,3,3),etai(msect,3,2),npts(3,msect),
& wt(3,msect,ndatmx),u(3,msect,ndatmx,3)
common nsect,sigma,qhati,eta,etai,wt,npts,u,f1,f2,f3,n1,n2,n3
c
do 50 i = 1,n1
do 50 j = 1,n1
p11(i,j)=p1(i,j)
50 p01(i,j)=p0(i,j)
c
c
do 60 i = 1,n2
do 60 j = 1,n2
p12(i,j)=p1(i+n1,j+n1)
60 p02(i,j)=p0(i+n1,j+n1)
c
c
do 70 i = 1,n3
do 70 j = 1,n3
p13(i,j)=p1(i+n1+n2,j+n1+n2)
70 p03(i,j)=p0(i+n1+n2,j+n1+n2)
c
c call trace1(p01,p11,p11,n1,t)
a = t/f1
call trace1(p02,p12,p12,n2,t)
a = a + t/f2
call trace1(p03,p13,p13,n3,t)
a = a + t/f3
call trace1(p01,p11,p01,n1,t)
b = t/f1
call trace1(p02,p12,p02,n2,t)
b = b + t/f2
call trace1(p03,p13,p03,n3,t)
b = b + t/f3
call trace2(p01,p11,n1,t)
c=t/f1
call trace2(p02,p12,n2,t)
c = c + t/f2
call trace2(p03,p13,n3,t)
c = c + t/f3
write(6,*) 'a, b, c ',a,b,c
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine trace1(q1,q2,q3,ndim,t)
c
c calculates trace of the product (q1-q2)*(I-q3)
c
implicit double precision(a-h,o-z)
parameter (ndatmx=400)
dimension q1(ndatmx,ndatmx),q2(ndatmx,ndatmx),q3(ndatmx,ndatmx)
real iden(ndatmx,ndatmx)
c
do 20 i = 1,ndim
do 10 j = 1,ndim
10 iden(i,j)=0.
20 iden(i,i) = 1.0
t = 0.
do 30 i = 1,ndim
do 30 j = 1,ndim
30 t = t +(q1(i,j)-q2(i,j))*(iden(i,j)-q3(i,j))
c
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
subroutine trace2(q1,q2,ndim,t)
c
c calculates {trace(q1-q2)}**2 - trace{[q1-q2]**2}
c
implicit double precision (a-h,o-z)
parameter (ndatmx=400)
dimension q1(ndatmx,ndatmx),q2(ndatmx,ndatmx)
c
t = 0.
do 50 i = 1,ndim
do 50 j = 1,ndim
50 t = t + (q1(i,j)-q2(i,j))**2
t2 = 0.
do 60 i = 1,ndim
60 t2 = t2 + q1(i,i)-q2(i,i)
t = t2**2 - t
c
return
end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine critcalc(a,b,c,plev,cfact,df,cxf)
c
implicit double precision(a-h,o-z)
real*4 plev,df,xf
c
eq=3+2*a+2*b
vq=6+14*a+2*b+2*c
ex2=eq*eq
df = (12*vq+2*ex2)/(3*vq-2*ex2)
cfact = eq*(1-2/df)
write(6,*)' '
write(6,*) 'degrees of freedom ',df
write (6,*) 'multiple of F ',cfact
call xidf(plev,3.0,df,xf,ier)
cxf=cfact*xf
write (6,*) ' critical value cF(3,df)',cxf
write(6,*)' '
write(6,*)' '
c
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine inv(a,b)
c
c inverts a 3 by 3 matrix
c
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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine resultsout(xa,a,cov,covinv,plev,hatkapc,df,npts,
& nuside,nvside,cxf,fmin,fmax,rank)
c
implicit double precision(a-h,o-z)
dimension xa(3),a(3,3),cov(3,3),covinv(3,3)
character filnam*20
c
write(6,*) ' rotation from side ',nuside, ' to side ',nvside
write(6,*) ' '
write(6,*) ' latitude, longitude, angle'
write(*,'(3x,3f8.2)') xa
write(6,*) 'Rotation matrix'
do 55 i = 1,3
55 write(6,*)(a(i,j),j=1,3)
write(6,*) ' covariance matrix'
do 60 i = 1,3
60 write(6,*)(cov(i,j),j=1,3)
c
65 write(6,1065)
1065 format(/,' Enter name of the output file: ',$)
1000 format(a)
read(6,1000)filnam
open(unit=15,file=filnam,status='unknown')
write(15,*) 'Results from trikap, rotation from plate',nuside,
& ' to plate ',nvside
write(15,*)' '
write(15,*) ' fitted rotation--alat,along,rho: '
write(15,*) xa
write(15,*) 'plev'
write(15,*) plev,1.0,1.0
write(15,*) '3/cfact, degrees of freedom,crit.value,flag'
write(15,'(f12.6,2x,f9.4,2x,f9.4,2x,f3.1)') hatkapc,df,cxf,1.0
write(15,*) 'number of points,fmin,fmax,rank'
write(15,'(7x,i4,3f5.0)') npts,fmin,fmax,rank
write(15,*) 'ahat: '
do 70 i = 1,3
70 write(15,*)(a(i,j),j=1,3)
write(15,*) ' covariance matrix'
do 75 i = 1,3
75 write(15,*)(cov(i,j),j=1,3)
write(15,*) 'inverse covariance matrix'
do 80 i = 1,3
80 write(15,*) (covinv(i,j),j=1,3)
close(15)
return
c
end