712 lines
18 KiB
Fortran
712 lines
18 KiB
Fortran
c****************************************************************
|
|
subroutine svdvecfit(i_mp,i_rd,i_fp,r_vecin,r_vobs,r_cov,
|
|
+ i_np,r_a,r_at2,r_u,r_v,r_w,r_chisq,l_chisq)
|
|
|
|
c****************************************************************
|
|
c**
|
|
c** FILE NAME: svdvecfit.f
|
|
c**
|
|
c** DATE WRITTEN: 01/02/95
|
|
c**
|
|
c** PROGRAMMER: Scott Hensley
|
|
c**
|
|
c** FUNCTIONAL DESCRIPTION: This routine does a least squares fit
|
|
c** to a vector valued observation least squares problem.
|
|
c**
|
|
c** ROUTINES CALLED: gaussj,svbksb,svdcmp,funcs
|
|
c**
|
|
c** NOTES: funcs is a user supplied function giving the jacobian
|
|
c** of the observation parameters wrt to fit parameters. This routine
|
|
c** is a generalization of Numerical Recipes svdfit. Note that this
|
|
c** routine can also be used in a nonlinear least squares procedure
|
|
c** by iterating properly.
|
|
c**
|
|
c** Solves the least problem
|
|
c**
|
|
c** T -1 -1 T -1
|
|
c** A = (AMAT COV AMAT) (AMAT COV )VOBS
|
|
c**
|
|
c** where AMAT is the jacobain of the observations vs parameters,
|
|
c** COV is the covriance matrix of observations
|
|
c** and VOBS is the vector of observations.
|
|
c**
|
|
c** r_a should be passed in with current best estimate of values
|
|
c**
|
|
c** UPDATE LOG:
|
|
c**
|
|
c** 4/17/95 - Reversed order of r_vecin, r_vobs, and r_cov SJS
|
|
c** revmoved r_vt, cleaned up parameter list
|
|
c**
|
|
c*****************************************************************
|
|
|
|
implicit none
|
|
|
|
c PARAMETERS:
|
|
integer I_NPE !number of parameters to estimate = i_np
|
|
integer I_RDE !number of observations per point = i_rd
|
|
real*8 R_TOL,R_LAMBDA
|
|
parameter(I_NPE=7)
|
|
parameter(I_RDE=2)
|
|
parameter(R_TOL=1.0d-20)
|
|
parameter (R_LAMBDA=1.d0)
|
|
|
|
c INPUT VARIABLES:
|
|
integer i_mp !number of input points
|
|
integer i_rd !number of observations each point
|
|
integer i_fp !number of input parameters to func
|
|
integer i_np !number of parameters to solve for
|
|
|
|
real*8 r_vecin(i_fp,i_mp) !vector values for func
|
|
real*8 r_vobs(i_rd,i_mp) !vector of observations
|
|
real*8 r_cov(i_rd,i_rd,i_mp) !covariance matrix of observation
|
|
real*8 r_chisq(i_rd,0:i_mp) !chisq for solution and fit vs observation
|
|
real*8 r_a(i_np) !solution to least squares
|
|
!for each point
|
|
logical l_chisq !evaluate the chisq for this fit
|
|
|
|
c OUTPUT VARIABLES:
|
|
real*8 r_at2(i_np) !delta to add to previous solution
|
|
real*8 r_u(i_np,i_np) !svd matrix, orthogonal matrix
|
|
real*8 r_v(i_np,i_np) !svd matrix, orthogonal matrix
|
|
real*8 r_w(i_np) !svd matrix, diagonal matrix
|
|
|
|
c LOCAL VARIABLES:
|
|
integer i,j,k,i_pts
|
|
real*8 r_covtemp(I_RDE,I_RDE)
|
|
real*8 r_am(I_NPE,I_RDE)
|
|
real*8 r_amat(I_RDE,I_NPE)
|
|
real*8 r_ptot(I_NPE)
|
|
real*8 r_wmax,r_thres,r_b(I_RDE,1),r_chird(I_RDE)
|
|
|
|
integer i_paramest(I_NPE),i_usedata(I_RDE)
|
|
common/funcom3/i_paramest,i_usedata
|
|
|
|
c DATA STATEMENTS:
|
|
|
|
C FUNCTION STATEMENTS:
|
|
|
|
c PROCESSING STEPS:
|
|
|
|
c init some arrays
|
|
|
|
c write(*,*) ' '
|
|
c write(*,*) 'Inside SVDVECFIT'
|
|
c write(*,*) ' '
|
|
|
|
if (i_rd .ne. I_RDE) then
|
|
write(*,*) 'ERROR - i_rd not equal to I_RDE in SVDVECFIT'
|
|
stop
|
|
end if
|
|
if (i_np .ne. I_NPE) then
|
|
write(*,*) 'ERROR - i_np not equal to I_NPE in SVDVECFIT'
|
|
stop
|
|
end if
|
|
|
|
do i=1,i_np
|
|
do j=1,i_np
|
|
r_u(i,j) = 0.0
|
|
enddo
|
|
r_ptot(i) = 0.0
|
|
enddo
|
|
|
|
c loop over the input points
|
|
|
|
do i_pts=1,i_mp
|
|
|
|
c write(*,*) 'i_pts = ',i_pts
|
|
|
|
c invert the covariance matrix of the observation
|
|
|
|
do i=1,i_rd
|
|
do j=1,i_rd
|
|
r_covtemp(i,j) = r_cov(i,j,i_pts)
|
|
enddo
|
|
enddo
|
|
|
|
call gaussj(r_covtemp,i_rd,i_rd,r_b,1,1)
|
|
|
|
c get the required jacobian matrix
|
|
|
|
call funcs(i_pts,i_rd,i_fp,r_vecin(1,i_pts),i_np,r_a,r_amat)
|
|
|
|
c do i=1,i_rd
|
|
c do j=1,i_np
|
|
c write(*,*) 'i,j,r_amat = ',i,j,r_amat(i,j)
|
|
c enddo
|
|
c enddo
|
|
|
|
c multiply amat transpose by the inverse cov matrix
|
|
|
|
do i=1,i_np
|
|
do j=1,i_rd
|
|
r_am(i,j) = 0.0
|
|
do k=1,i_rd
|
|
r_am(i,j) = r_am(i,j) + r_amat(k,i)*r_covtemp(k,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
c do i=1,i_np
|
|
c do j=1,i_rd
|
|
c write(*,*) 'i,j,r_am = ',i,j,r_am(i,j)
|
|
c enddo
|
|
c enddo
|
|
|
|
c multiply am by amat
|
|
|
|
do i=1,i_np
|
|
do j=1,i_np
|
|
do k=1,i_rd
|
|
r_u(i,j) = r_u(i,j) + r_am(i,k)*r_amat(k,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
c multilpy am by vobs
|
|
|
|
|
|
c write(*,*) 'r_vobs,i_pts = ',i_pts,r_vobs(1,i_pts),r_vobs(2,i_pts)
|
|
do i=1,i_np
|
|
do k=1,i_rd
|
|
r_ptot(i) = r_ptot(i) + r_am(i,k)*r_vobs(k,i_pts)
|
|
enddo
|
|
enddo
|
|
|
|
enddo !i_pts
|
|
|
|
c find the SVD of the r_u matrix
|
|
|
|
c do i=1,i_np
|
|
c do j=1,i_np
|
|
c write(*,*) 'i,j,r_u = ',i,j,r_u(i,j)
|
|
c enddo
|
|
c enddo
|
|
|
|
call svdcmp(r_u,i_np,i_np,i_np,i_np,r_w,r_v)
|
|
|
|
c do i=1,i_np
|
|
c do j=1,i_np
|
|
c write(*,*) 'i,j,r_u,r_v = ',i,j,r_u(i,j),r_v(i,j)
|
|
c enddo
|
|
c enddo
|
|
|
|
c do i=1,i_np
|
|
c write(*,*) 'w = ',i,r_w(i)
|
|
c enddo
|
|
|
|
c kill off all the singular values
|
|
|
|
r_wmax = 0.0
|
|
do i=1,i_np
|
|
if(r_w(i) .gt. r_wmax)then
|
|
r_wmax = r_w(i)
|
|
endif
|
|
enddo
|
|
r_thres = r_wmax*R_TOL
|
|
c write(*,*) 'r_thres = ',r_thres
|
|
|
|
do i=1,i_np
|
|
if(r_w(i) .lt. r_thres)then
|
|
r_w(i) = 0.0
|
|
endif
|
|
enddo
|
|
|
|
c do i=1,i_np
|
|
c write(*,*) 'w = ',i,r_w(i)
|
|
c enddo
|
|
|
|
c use the svbksb routine to solve for the desired parameters
|
|
|
|
call svbksb(r_u,r_w,r_v,i_np,i_np,i_np,i_np,r_ptot,r_at2)
|
|
|
|
c update the r_a vector
|
|
|
|
do i=1,i_np
|
|
r_at2(i) = -r_at2(i)*i_paramest(i)
|
|
r_a(i) = r_at2(i)/R_LAMBDA + r_a(i)
|
|
c write(*,*) 'a=',i,r_a(i),r_at2(i)
|
|
enddo
|
|
|
|
c evaluate the chisq array (linearized version)
|
|
|
|
if(l_chisq)then
|
|
|
|
c loop over data points
|
|
|
|
|
|
do i=1,i_rd
|
|
r_chird(i) = 0.
|
|
enddo
|
|
r_chisq(1,0) = 0.0
|
|
do i=1,i_mp
|
|
|
|
call funcs(i,i_rd,i_fp,r_vecin(1,i),i_np,r_a,r_amat)
|
|
|
|
do j=1,i_rd
|
|
r_chisq(j,i) = 0.0
|
|
do k=1,i_np
|
|
r_chisq(j,i) = r_chisq(j,i) + r_amat(j,k)*r_at2(k)
|
|
enddo
|
|
c write(*,*) 'r_chisq = ',i,j,r_chisq(j,i),r_vobs(j,i)
|
|
r_chisq(j,i) = r_covtemp(j,j)*(r_chisq(j,i) -
|
|
+ r_vobs(j,i))**2
|
|
r_chisq(1,0) = r_chisq(1,0) + r_chisq(j,i)
|
|
r_chird(j) = r_chird(j) + r_chisq(j,i)
|
|
enddo
|
|
|
|
enddo !i_pts loop for chisq
|
|
|
|
r_chisq(1,0) = sqrt(r_chisq(1,0)/(2.*i_mp))
|
|
write(*,*) 'r_chisq = ',r_chisq(1,0),sqrt(r_chird(1)/i_mp),sqrt(r_chird(2)/i_mp)
|
|
|
|
endif
|
|
|
|
end
|
|
|
|
c******************************************************************************
|
|
|
|
SUBROUTINE gaussj(a,n,np,b,m,mp)
|
|
INTEGER m,mp,n,np,NMAX
|
|
REAL*8 a(np,np),b(np,mp)
|
|
PARAMETER (NMAX=50)
|
|
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
|
|
REAL*8 big,dum,pivinv
|
|
do 11 j=1,n
|
|
ipiv(j)=0
|
|
11 continue
|
|
do 22 i=1,n
|
|
big=0.
|
|
do 13 j=1,n
|
|
if(ipiv(j).ne.1)then
|
|
do 12 k=1,n
|
|
if (ipiv(k).eq.0) then
|
|
if (abs(a(j,k)).ge.big)then
|
|
big=abs(a(j,k))
|
|
irow=j
|
|
icol=k
|
|
endif
|
|
else if (ipiv(k).gt.1) then
|
|
pause 'singular matrix in gaussj'
|
|
endif
|
|
12 continue
|
|
endif
|
|
13 continue
|
|
ipiv(icol)=ipiv(icol)+1
|
|
if (irow.ne.icol) then
|
|
do 14 l=1,n
|
|
dum=a(irow,l)
|
|
a(irow,l)=a(icol,l)
|
|
a(icol,l)=dum
|
|
14 continue
|
|
do 15 l=1,m
|
|
dum=b(irow,l)
|
|
b(irow,l)=b(icol,l)
|
|
b(icol,l)=dum
|
|
15 continue
|
|
endif
|
|
indxr(i)=irow
|
|
indxc(i)=icol
|
|
if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
|
|
pivinv=1./a(icol,icol)
|
|
a(icol,icol)=1.
|
|
do 16 l=1,n
|
|
a(icol,l)=a(icol,l)*pivinv
|
|
16 continue
|
|
do 17 l=1,m
|
|
b(icol,l)=b(icol,l)*pivinv
|
|
17 continue
|
|
do 21 ll=1,n
|
|
if(ll.ne.icol)then
|
|
dum=a(ll,icol)
|
|
a(ll,icol)=0.
|
|
do 18 l=1,n
|
|
a(ll,l)=a(ll,l)-a(icol,l)*dum
|
|
18 continue
|
|
do 19 l=1,m
|
|
b(ll,l)=b(ll,l)-b(icol,l)*dum
|
|
19 continue
|
|
endif
|
|
21 continue
|
|
22 continue
|
|
do 24 l=n,1,-1
|
|
if(indxr(l).ne.indxc(l))then
|
|
do 23 k=1,n
|
|
dum=a(k,indxr(l))
|
|
a(k,indxr(l))=a(k,indxc(l))
|
|
a(k,indxc(l))=dum
|
|
23 continue
|
|
endif
|
|
24 continue
|
|
return
|
|
END
|
|
|
|
SUBROUTINE svdcmp(a,m,n,mp,np,w,v)
|
|
INTEGER m,mp,n,np,NMAX
|
|
REAL*8 a(mp,np),v(np,np),w(np)
|
|
PARAMETER (NMAX=500)
|
|
INTEGER i,its,j,jj,k,l,nm
|
|
REAL*8 anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythag
|
|
real*8 r_one
|
|
|
|
g=0.0
|
|
r_one = 1.d0
|
|
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
|
|
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
|
|
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
|
|
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
|
|
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
|
|
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
|
|
do 39 i=min(m,n),1,-1
|
|
l=i+1
|
|
g=w(i)
|
|
do 33 j=l,n
|
|
a(i,j)=0.0
|
|
33 continue
|
|
if(g.ne.0.0)then
|
|
g=1.0/g
|
|
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
|
|
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
|
|
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) goto 2
|
|
if((abs(w(nm))+anorm).eq.anorm) goto 1
|
|
41 continue
|
|
1 c=0.0
|
|
s=1.0
|
|
do 43 i=l,k
|
|
f=s*rv1(i)
|
|
rv1(i)=c*rv1(i)
|
|
if((abs(f)+anorm).eq.anorm) goto 2
|
|
g=w(i)
|
|
h=pythag(f,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
|
|
43 continue
|
|
2 z=w(k)
|
|
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
|
|
goto 3
|
|
endif
|
|
if(its.eq.30) pause 'no convergence in svdcmp'
|
|
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=pythag(f,r_one)
|
|
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
|
|
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=pythag(f,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 jj=1,n
|
|
x=v(jj,j)
|
|
z=v(jj,i)
|
|
v(jj,j)= (x*c)+(z*s)
|
|
v(jj,i)=-(x*s)+(z*c)
|
|
45 continue
|
|
z=pythag(f,h)
|
|
w(j)=z
|
|
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 jj=1,m
|
|
y=a(jj,j)
|
|
z=a(jj,i)
|
|
a(jj,j)= (y*c)+(z*s)
|
|
a(jj,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
|
|
return
|
|
END
|
|
|
|
REAL*8 FUNCTION pythag(a,b)
|
|
REAL*8 a,b
|
|
REAL*8 absa,absb
|
|
absa=abs(a)
|
|
absb=abs(b)
|
|
if(absa.gt.absb)then
|
|
pythag=absa*sqrt(1.d0+(absb/absa)**2)
|
|
else
|
|
if(absb.eq.0.)then
|
|
pythag=0.
|
|
else
|
|
pythag=absb*sqrt(1.d0+(absa/absb)**2)
|
|
endif
|
|
endif
|
|
return
|
|
END
|
|
|
|
SUBROUTINE svbksb(u,w,v,m,n,mp,np,b,x)
|
|
INTEGER m,mp,n,np,NMAX
|
|
REAL*8 b(mp),u(mp,np),v(np,np),w(np),x(np)
|
|
PARAMETER (NMAX=500)
|
|
INTEGER i,j,jj
|
|
REAL*8 s,tmp(NMAX)
|
|
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
|
|
s=s/w(j)
|
|
endif
|
|
tmp(j)=s
|
|
12 continue
|
|
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
|
|
return
|
|
END
|
|
|
|
SUBROUTINE svdvar(v,ma,np,w,cvm,ncvm)
|
|
INTEGER ma,ncvm,np,MMAX
|
|
REAL*8 cvm(ncvm,ncvm),v(np,np),w(np)
|
|
PARAMETER (MMAX=20)
|
|
INTEGER i,j,k
|
|
REAL*8 sum,wti(MMAX)
|
|
do 11 i=1,ma
|
|
wti(i)=0.
|
|
if(w(i).ne.0.) wti(i)=1.d0/(w(i)*w(i))
|
|
11 continue
|
|
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 Modify Numerical Recipes program moment.f to compute only
|
|
c standard deviation and allow double precision
|
|
SUBROUTINE moment(data,p,sdev)
|
|
Implicit None
|
|
INTEGER p
|
|
REAL*8 adev,ave,curt,sdev,skew,var,data(p)
|
|
INTEGER j
|
|
REAL*8 t,s,ep
|
|
if(p.le.1)pause 'p must be at least 2 in moment'
|
|
s=0.0d0
|
|
do 11 j=1,p
|
|
s=s+data(j)
|
|
11 continue
|
|
ave=s/p
|
|
adev=0.0d0
|
|
var=0.0d0
|
|
skew=0.0d0
|
|
curt=0.0d0
|
|
ep=0.
|
|
do 12 j=1,p
|
|
s=data(j)-ave
|
|
t=s*s
|
|
var=var+t
|
|
12 continue
|
|
adev=adev/p
|
|
var=(var-ep**2/p)/(p-1)
|
|
sdev=sqrt(var)
|
|
return
|
|
END
|
|
|
|
c This program is used to find the rotation matrix from the affine matrix
|
|
SUBROUTINE qrdcmp(a,n,np,c,d,sing)
|
|
INTEGER n,np
|
|
REAL*8 a(np,np),c(n),d(n)
|
|
LOGICAL sing
|
|
INTEGER i,j,k
|
|
REAL*8 scale,sigma,sum,tau
|
|
sing=.false.
|
|
scale=0.
|
|
do 17 k=1,n-1
|
|
do 11 i=k,n
|
|
scale=max(scale,abs(a(i,k)))
|
|
11 continue
|
|
if(scale.eq.0.)then
|
|
sing=.true.
|
|
c(k)=0.
|
|
d(k)=0.
|
|
else
|
|
do 12 i=k,n
|
|
a(i,k)=a(i,k)/scale
|
|
12 continue
|
|
sum=0.
|
|
do 13 i=k,n
|
|
sum=sum+a(i,k)**2
|
|
13 continue
|
|
sigma=sign(sqrt(sum),a(k,k))
|
|
a(k,k)=a(k,k)+sigma
|
|
c(k)=sigma*a(k,k)
|
|
d(k)=-scale*sigma
|
|
do 16 j=k+1,n
|
|
sum=0.
|
|
do 14 i=k,n
|
|
sum=sum+a(i,k)*a(i,j)
|
|
14 continue
|
|
tau=sum/c(k)
|
|
do 15 i=k,n
|
|
a(i,j)=a(i,j)-tau*a(i,k)
|
|
15 continue
|
|
16 continue
|
|
endif
|
|
17 continue
|
|
d(n)=a(n,n)
|
|
if(d(n).eq.0.)sing=.true.
|
|
return
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software $23#1yR.3Z9.
|