- 积分
- 3268
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-8-2
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2015-11-20 16:58:13
|
显示全部楼层
FORTRAN程序:
program main
parameter(m=105,n=30,mnl=n,np=4,ll=m)
parameter(ks=0,undef=-999.9)
!-----input array
dimension f0(m,n),F(M,N)
!-----work array
dimension gvt(ll,mnl),rgvt(ll,np),cof(mnl,n),rcof(np,n)
!-----output arrays
dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
dimension rer(np,2),regvt(m,np),recof(np,n)
c-----read data
open(1,file='njrr12.txt')
do i=1,m
c read(1,101)(f0(i,j),J=1,N)
read(1,*)(f0(i,j),j=1,n)
end do
101 FORMAT(33F10.9)
close(1)
c-----Remove the terrain or missing value.
l=0
do j=1,m
if(f0(j,1).ne.undef)then
l=l+1
do k=1,n
f(l,k)=f0(j,k)
enddo
endif
enddo
write(*,*)'ll=',l
c-----call the subroutine
call reof(ll,n,mnl,np,f,ks,er,gvt,ecof,rer,rgvt,recof)
c-----Add the terrain or missing value.
l=0
do i=1,m
if(f0(i,1).ne.undef)then
l=l+1
do k=1,n
egvt(i,k)=gvt(l,k)
end do
do k=1,np
regvt(i,k)=rgvt(l,k)
end do
else
do k=1,n
egvt(i,k)=undef
end do
do k=1,np
regvt(i,k)=undef
end do
endif
enddo
!-----output the result
c-----output the error
ccccccccccccccccccccccccccccccc
open(99,file='bzh-l.txt')
write(99,*)((f(i,j),j=1,n),i=1,m)
close(99)
ccccccccccccccccccccccccccccccc
OPEN(10,file='er.dat')
! write(10,*)'error=',sum-tr
write(10,*)'EOF lanmda (eigenvalues) from big to small'
write(10,*)(er(i,1),i=1,mnl)
write(*,*)'lamda='
write(*,*)(er(i,1),i=1,mnl)
write(10,*)'EOF accumulated eigenvalues from big to small'
write(10,*)(er(i,2),i=1,mnl)
write(10,*)'EOF explained variances'
write(10,*)(er(i,3),i=1,mnl)
write(10,*)'EOF accumulated explained variances'
write(10,*)(er(i,4),i=1,mnl)
write(10,*)'REOF explained variances'
write(10,*)(rer(i,1),i=1,np)
write(10,*)'REOF accumulated explained variances'
write(10,*)(rer(i,2),i=1,np)
c-----output eigenvactors matrix of EOF
OPEN(11,FILE='egvt.dat',form='binary')
do i=1,m
write(11)(egvt(i,j),j=1,10)
end do
close(11)
OPEN(11,FILE='egvt.txt')
do i=1,m
write(11,102)(egvt(i,j),j=1,10)
enddo
102 format(10f10.4)
close(11)
c-----output time coefficients matrix of EOF
open(14,file='ecof.dat',form='binary')
do i=1,n
write(14)(ecof(j,i),j=1,10)
end do
close(14)
open(14,file='ecof-l.txt')
do i=1,n
write(14,102)(ecof(j,i),j=1,10)
end do
close(14)
c-----output loading vectors of REOF
OPEN(21,FILE='regvt.dat',form='binary')
do i=1,m
write(21)(regvt(i,j),j=1,20)
end do
close(21)
OPEN(21,FILE='regvt.txt')
do i=1,m
write(21,103)(regvt(i,j),j=1,20)
end do
103 format(20f10.4)
close(21)
c-----output time coefficients matrix of REOF
OPEN(22,FILE='recof.dat',form='binary')
do j=1,n
write(22)(recof(i,j),i=1,np)
end do
close(22)
OPEN(22,FILE='recof.txt')
do j=1,n
write(22,104)(recof(i,j),i=1,np)
end do
104 format(33f10.4)
close(22)
stop
end
c-----
subroutine reof(m,n,mnl,np,f,ks,er,egvt,ecof,rer,regvt,recof)
c-----input array
dimension f(m,n)
c-----work arrays
dimension cov(mnl,mnl),s(mnl,mnl),d(mnl)
c-----output arrays
dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
dimension rer(np,2),regvt(m,np),recof(np,n)
c---- Preprocessing data
call transf(m,n,f,ks)
c---- Crossed product matrix of the data f
call crossproduct(m,n,mnl,f,cov,sum)
c---- Eigenvalues and eigenvectors by the Jacobi method
eps=1e-7
call jcb(mnl,cov,s,d,eps)
c---- Specified eigenvalues
call arrang(m,n,mnl,d,s,er,tr)
write(*,*)'error=',sum-tr
c---- Normalized eignvectors and their time coefficients
call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
write(*,*)'EOF accomplished'
c-----linear transposed (rotate)
call rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
write(*,*)'REOF accomplished'
c---- Specified eignvectors and time coefficients with
c explained variances of REOF (if you need)
call arrange2(m,n,np,regvt,recof,rer)
return
end
subroutine transf(m,n,f,ks)
!-----input array
dimension f(m,n)
!-----work array
dimension fw(n)
if(ks.eq.-1)return
if(ks.eq.0)then !anomaly of f
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
do j=1,n
f(i,j)=f(i,j)-af
enddo
enddo
return
endif
if(ks.eq.1)then !normalizing f
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
do j=1,n
f(i,j)=(f(i,j)-af)/sf
enddo
enddo
endif
return
end
c-----
subroutine crossproduct(m,n,mnl,f,cov,sum)
!-----input array
dimension f(m,n)
!-----work array
dimension cov(mnl,mnl)
if(n-m) 10,50,50
10 do 20 i=1,mnl
do 20 j=i,mnl
cov(i,j)=0.0
do is=1,m
cov(i,j)=cov(i,j)+f(is,i)*f(is,j)
enddo
cov(j,i)=cov(i,j)
20 continue
sum=0.
do i=1,mnl
sum=sum+cov(i,i)
end do
sum=sum/(mnl*1.)
return
50 do 60 i=1,mnl
do 60 j=i,mnl
cov(i,j)=0.0
do js=1,n
cov(i,j)=cov(i,j)+f(i,js)*f(j,js)
enddo
cov(j,i)=cov(i,j)
60 continue
sum=0.
do i=1,mnl
sum=sum+cov(i,i)
end do
sum=sum/(mnl*1.)
return
end
c-----
subroutine jcb(m,a,s,d,eps)
!-----input array
dimension a(m,m)
!-----output arrays
dimension s(m,m),d(m)
do 30 i=1,m
do 30 j=1,i
if(i-j) 20,10,20
10 s(i,j)=1.
go to 30
20 s(i,j)=0.
s(j,i)=0.
30 continue
g=0.
do 40 i=2,m
i1=i-1
do 40 j=1,i1
40 g=g+2.*a(i,j)*a(i,j)
s1=sqrt(g)
s2=eps/float(m)*s1
s3=s1
l=0
50 s3=s3/float(m)
60 do 130 iq=2,m
iq1=iq-1
do 130 ip=1,iq1
if(abs(a(ip,iq)).lt.s3) goto 130
l=1
v1=a(ip,ip)
v2=a(ip,iq)
v3=a(iq,iq)
u=0.5*(v1-v3)
if(u.eq.0.0) g=1.
if(abs(u).ge.1e-10) g=-sign(1.,u)*v2/sqrt(v2*v2+u*u)
st=g/sqrt(2.*(1.+sqrt(1.-g*g)))
ct=sqrt(1.-st*st)
do 110 i=1,m
g=a(i,ip)*ct-a(i,iq)*st
a(i,iq)=a(i,ip)*st+a(i,iq)*ct
a(i,ip)=g
g=s(i,ip)*ct-s(i,iq)*st
s(i,iq)=s(i,ip)*st+s(i,iq)*ct
110 s(i,ip)=g
do 120 i=1,m
a(ip,i)=a(i,ip)
120 a(iq,i)=a(i,iq)
g=2.*v2*st*ct
a(ip,ip)=v1*ct*ct+v3*st*st-g
a(iq,iq)=v1*st*st+v3*ct*ct+g
a(ip,iq)=(v1-v3)*st*ct+v2*(ct*ct-st*st)
a(iq,ip)=a(ip,iq)
130 continue
if(l-1) 150,140,150
140 l=0
go to 60
150 if(s3.gt.s2) goto 50
do 160 i=1,m
d(i)=a(i,i)
160 continue
return
end
c-----
subroutine arrang(m,n,mnl,d,s,er,tr)
!-----input arrays
dimension d(mnl),s(mnl,mnl)
!-----output array
dimension er(mnl,4)
tr=0.0
do 10 i=1,mnl
tr=tr+d(i)/(mnl*1.)
er(i,1)=d(i)/(mnl*1.)
10 continue
mnl1=mnl-1
do 20 k1=mnl1,1,-1
do 20 k2=k1,mnl1
if(er(k2,1).lt.er(k2+1,1)) then
c=er(k2+1,1)
er(k2+1,1)=er(k2,1)
er(k2,1)=c
do 15 i=1,mnl
c=s(i,k2+1)
s(i,k2+1)=s(i,k2)
s(i,k2)=c
15 continue
endif
20 continue
er(1,2)=er(1,1)
do 30 i=2,mnl
er(i,2)=er(i-1,2)+er(i,1)
30 continue
do 40 i=1,mnl
er(i,3)=er(i,1)*100./tr
er(i,4)=er(i,2)*100./tr
40 continue
return
end
c-----
subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
!-----input arrays
dimension f(m,n),s(mnl,mnl),er(mnl,4)
!-----work arrays
dimension v(m),v1(n)
!-----output arrays
dimension egvt(m,mnl),ecof(mnl,n)
do j=1,mnl
do i=1,m
egvt(i,j)=0.
enddo
do i=1,n
ecof(j,i)=0.
enddo
enddo
c-----Normalizing the input eignvectors (c=sqrt(lamda))
do 10 j=1,mnl
c=0.
do i=1,mnl
c=c+s(i,j)*s(i,j)
enddo
c=sqrt(c)
do i=1,mnl
s(i,j)=s(i,j)/c
enddo
10 continue
c-----(m<n)
if(m.le.n) then
do js=1,mnl
do i=1,m
egvt(i,js)=s(i,js)
enddo
enddo
do 30 j=1,n
do i=1,m
v(i)=f(i,j)
enddo
do is=1,mnl
do i=1,m
ecof(is,j)=ecof(is,j)+v(i)*egvt(i,is)
enddo
enddo
30 continue
else
c-----(m>n)
do 40 i=1,m
do j=1,n
v1(j)=f(i,j)
enddo
do js=1,mnl
do j=1,n
egvt(i,js)=egvt(i,js)+v1(j)*s(j,js)
enddo
enddo
40 continue
do 50 js=1,mnl
do j=1,n
ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
enddo
do i=1,m
egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
enddo
50 continue
do js=1,mnl
do i=1,m
egvt(i,js)=(sqrt(abs(er(js,1)/(mnl*1.0))))*egvt(i,js)
end do
do i=1,n
ecof(js,i)=(ecof(js,i))/(sqrt(abs(er(js,1)/(mnl*1.0))))
end do
end do
end if
return
end
c-----
subroutine rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
c-----input arrays
dimension egvt(m,mnl),ecof(mnl,n)
c-----work arrays
dimension h(m),st(np),vrlv(50),er(mnl,4),s(m,np)
c-----output arrays
dimension regvt(m,np),recof(np,n),rer(np,2)
integer t,k
write(*,*)'rotate beginning'
c-----Take out forward NP eigenvectors to rotate
do j=1,np
do i=1,m
regvt(i,j)=egvt(i,j)
end do
end do
c-----Take out forward NP time coefficients to rotate
do i=1,np
do j=1,n
recof(i,j)=ecof(i,j)
end do
end do
c-----compute the common degree
do i=1,m
h(i)=0.0
do j=1,np
h(i)=h(i)+regvt(i,j)**2
end do
end do
do i=1,m
h(i)=sqrt(h(i))
end do
c-----compute the variance of the variance contribution by per common degree
c-----(equal to the ratated one)
vrlv0=0.0
do k=1,np
g1=0.0
g2=0.0
do i=1,m
g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
end do
g2=g2**2
vrlv0=vrlv0+g1-g2
end do
c-----rotated
lll=0
write(*,*)'rotate times'
800 continue
do 505 t=1,np
do 505 k=1,np
if(t.ge.k) goto 505
call rot(regvt,recof,h,t,k,m,n,np)
505 continue
lll=lll+1
write(*,*)'LLL=',LLL
vrlv(lll)=0.0
do k=1,np
g1=0.0
g2=0.0
do i=1,m
g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
end do
g2=g2**2
vrlv(lll)=vrlv(lll)+g1-g2
end do
if(lll.lt.2) goto 800
ci=(vrlv(lll)-vrlv(lll-1))/vrlv0
if(ci.lt.0.001) goto 600
if(lll.ge.40) goto 600
goto 800
600 continue
c-----compute the explained variances
do j=1,np
st(j)=0.0
do i=1,m
st(j)=st(j)+regvt(i,j)**2
end do
end do
do j=1,np
rer(j,1)=st(j)*100/tr
end do
c-----compute the accumulated explained variances
ddd=0.0
do k=1,np
ddd=ddd+st(k)/tr
rer(k,2)=ddd*100.
enddo
return
end
c-----
subroutine rot(a,b,h,t,k,m,n,np)
c-----input arrays
dimension a(m,np),b(np,n),h(m)
c-----work arrays
dimension u(m),vr(m),wt(m),wk(m),wbt(n),wbk(n)
c-----output arrays
c dimension a(m,np),b(np,n)
integer t,k
c-----compute fi
do 25 i=1,m
u(i)=(a(i,t)/h(i))**2-(a(i,k)/h(i))**2
25 vr(i)=2.0*(a(i,t)/h(i))*(a(i,k)/h(i))
c=0.0
d=0.0
s=0.0
g=0.0
do 30 i=1,m
c=c+(u(i)**2-vr(i)**2)
d=d+2.0*u(i)*vr(i)
s=s+u(i)
30 g=g+vr(i)
tg1=d-2.0*s*g/(m*1.)
tg2=c-(s**2-g**2)/(m*1.)
fi=atan2(tg1,tg2)/4.0
c-----compute new a and b with fi
do 75 i=1,m
wt(i)=a(i,t)
75 wk(i)=a(i,k)
do 84 kk=1,n
wbt(kk)=b(t,kk)
84 wbk(kk)=b(k,kk)
do 78 i=1,m
a(i,t)=wt(i)*cos(fi)+wk(i)*sin(fi)
78 a(i,k)=wt(i)*(-sin(fi))+wk(i)*cos(fi)
do 89 it=1,n
b(t,it)=wbt(it)*cos(fi)+wbk(it)*sin(fi)
89 b(k,it)=wbt(it)*(-sin(fi))+wbk(it)*cos(fi)
return
end
c-----
subroutine arrange2(m,n,np,regvt,recof,rer)
c-----input and output arrays
dimension regvt(m,np),recof(np,n),rer(np,2)
mnl1=np-1
do 20 k1=mnl1,1,-1
do 20 k2=k1,mnl1
if(rer(k2,1).lt.rer(k2+1,1)) then
c=rer(k2+1,1)
rer(k2+1,1)=rer(k2,1)
rer(k2,1)=c
do i=1,m
d=regvt(i,k2+1)
regvt(i,k2+1)=regvt(i,k2)
regvt(i,k2)=d
end do
do i=1,n
e=recof(k2+1,i)
recof(k2+1,i)=recof(k2,i)
recof(k2,i)=e
end do
endif
20 continue
do i=1,np
rer(i,2)=0.0
end do
rer(1,2)=rer(1,1)
do 30 i=2,np
rer(i,2)=rer(i-1,2)+rer(i,1)
30 continue
return
end
subroutine meanvar(n,x,ax,sx,vx)
!-----input array
dimension x(n)
ax=0.
vx=0.
sx=0.
do 10 i=1,n
ax=ax+x(i)
10 continue
ax=ax/float(n)
do 20 i=1,n
vx=vx+(x(i)-ax)**2
20 continue
vx=vx/float(n)
sx=sqrt(vx)
return
end
|
|