c SVD analysis (mx >= nx) c============================================================= subroutine svdsub(mx,nx,nt,md,f1,f2,v1,v2,a,b,s) real v1(mx,md),v2(nx,md),a(md,nt),b(md,nt),s(nx) parameter (m=794,n=318,ntt=70,mode=10) c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ real left,right,stdl,stdr,f1(mx,nt),f2(nx,nt) real svlue,vecl,vecr real ak,bk,cf common/data1/left(ntt,m),right(ntt,n) common/data2/stdl(m),stdr(n) common/data3/svlue(n),vecl(m,n),vecr(n,n) common/data4/ak(ntt,mode),bk(ntt,mode),cf(mode) do it=1,nt do i=1,mx left(it,i)=f1(i,it) if (f1(i,it).eq.-8888.8) left(it,i)=0. end do do j=1,nx right(it,j)=f2(j,it) if (f2(j,it).eq.-8888.8) right(it,j)=0. end do end do c left field: left (anomaly) ----> stdl (standard deviation) do i=1,mx stdl(i)=0. do it=1,nt stdl(i)=stdl(i) + left(it,i)**2 end do stdl(i)=sqrt(stdl(i)/real(nt)) end do c print *,stdl c right field: right -----> stdr do i=1,nx stdr(i)=0. do it=1,nt stdr(i)=stdr(i) + right(it,i)**2 end do stdr(i)=sqrt(stdr(i)/real(nt)) end do c without the following do-loop c svd is based on covariance matrix c************************************* c do it=1,nt c do i=1,mx c left(it,i)=left(it,i)/stdl(i) c end do c do j=1,nx c right(it,j)=right(it,j)/stdr(j) c end do c end do c compute covariance matrix c print*,'computing cov_matrix...' call cov_mat(mx,nx,nt) c SVD analysis of the matrix c print*,'doing SVD...' call solv_svd(mx,nx,md,nt) do k=1,md do i=1,mx c v1(i,k)=vecl(i,k) *sqrt(svlue(k)) v1(i,k)=vecl(i,k)*sqrt(svlue(k)) end do do i=1,nx c v2(i,k)=vecr(i,k) *sqrt(svlue(k)) v2(i,k)=vecr(i,k)*sqrt(svlue(k)) end do end do do k=1,md do i=1,nt a(k,i)=ak(i,k) b(k,i)=bk(i,k) end do end do do i=1,nx s(i)=svlue(i) end do return end c========================================================== subroutine cov_mat(mx,nx,nt) parameter (m=794,n=318,ntt=70,mode=10) c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ real left,right,stdl,stdr real svlue,vecl,vecr common/data1/left(ntt,m),right(ntt,n) common/data2/stdl(m),stdr(n) common/data3/svlue(n),vecl(m,n),vecr(n,n) rnt=1./float(nt-1) do j=1,nx do i=1,mx cij=0. do k=1,nt cij=cij+left(k,i)*right(k,j) enddo vecl(i,j)=cij*rnt enddo enddo return end c========================================================= c SVD analysis of mp by np matrix vecl, return matrice c vecl(mp by np),svlue(np), and vecr(np by np), where c vecl and vecr are column-orthogonal matrice. c========================================================= subroutine solv_svd(mp,np,md,nt) parameter (m=794,n=318,ntt=70,mode=10) c ^^^^^^^^^^^^^^^^^^^^^^^^^^^ real left,right,stdl,stdr real svlue,vecl,vecr real ak,bk,cf common/data1/left(ntt,m),right(ntt,n) common/data2/stdl(m),stdr(n) common/data3/svlue(n),vecl(m,n),vecr(n,n) common/data4/ak(ntt,mode),bk(ntt,mode),cf(mode) dimension rakbk(mode) c do svd call svdcmp(vecl,mp,np,m,n,svlue,vecr) c sort out U & V according to sigular value W call svdsrt(svlue,vecl,vecr,mp,np,m,n) c get total CF cfsum=0. do i=1,np c cfsum=cfsum+svlue(i) cfsum=cfsum+svlue(i)**2. c ^^^^^^^^^^ huiwang 2/12/95 enddo c print*,' ' c print*,'Total cov=',cfsum c get CF do i=1,md c cf(i)=svlue(i)/cfsum * 100 cf(i)=svlue(i)**2./cfsum * 100 c ^^^^^ <-------------huiwang 2/12/95 enddo c print*,' CF(%):',(cf(i),i=1,md) c check orthogonality do i=1,md do j=1,md oa=0. do k=1,mp oa=oa + vecl(k,i)*vecl(k,j) end do c print *, 'left orthogonal?',oa end do end do do i=1,md do j=1,md oa=0. do k=1,np oa=oa + vecr(k,i)*vecr(k,j) end do c print *, 'right orthogonal?',oa end do end do c get time coefficients c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do k=1,nt do j=1,md tmp=0. do i=1,mp tmp=tmp+vecl(i,j)*left(k,i) enddo ak(k,j)=tmp tmp=0. do i=1,np tmp=tmp+vecr(i,j)*right(k,i) enddo bk(k,j)=tmp enddo enddo c c output for projection c cc open (66,file='prjt_JJA_us-atl_md1.gr', cc $ access='direct',form='unformatted',recl=np*4) cc write(66,rec=1) (vecr(i,1),i=1,np) cc open (67,file='prjt_JJA_us-atl_md2.gr', cc $ access='direct',form='unformatted',recl=np*4) cc write(67,rec=1) (vecr(i,2),i=1,np) c-------------------------------------------------------------- c total variance for left and right fields <-- huiwang 2.13.95 tvar1=0. do i=1,mp tv1=0. do k=1,nt tv1=tv1 + left(k,i)**2 end do tvar1=tvar1 + tv1/real(nt) end do tvar2=0. do i=1,np tv2=0. do k=1,nt tv2=tv2 + right(k,i)**2 end do tvar2=tvar2 + tv2/real(nt) end do print *, 'total variance LEFT=',tvar1,' RIGHT=',tvar2 c variance for single svd mode ssv=0. do j=1,md svar1=0. do i=1,mp sv1=0. do k=1,nt sv1=sv1 + (ak(k,j)*vecl(i,j))**2 end do svar1=svar1 + sv1/real(nt) end do print *,'LEFT MODE',j,' vari=',svar1,' ',svar1/tvar1*100,'%' ssv=ssv + svar1 end do print *,'LEFT vari=',ssv,' ',ssv/tvar1*100,'%' ssv=0. do j=1,md svar2=0. do i=1,np sv2=0. do k=1,nt sv2=sv2 + (bk(k,j)*vecr(i,j))**2 end do svar2=svar2 + sv2/real(nt) end do print *,'RIGHT MODE',j,' vari=',svar2,' ',svar2/tvar2*100,'%' ssv= ssv + svar2 end do print *,'RIGHT vari=',ssv,' ',ssv/tvar2*100,'%' print*,' ' print*,'Total cov=',cfsum c check reconstructed data at one grid point c print *,'check reconstructed data' i=50 do k=1,nt sv1=0. do j=1,md sv1=sv1 + bk(k,j)*vecr(i,j) end do c print *, sv1,' ',right(k,i) end do c correlation b/t ak & bk do i=1,md aa=0. bb=0. cc=0. do it=1,nt aa=aa + ak(it,i)**2 bb=bb + bk(it,i)**2 cc=cc + ak(it,i)*bk(it,i) end do rakbk(i)=cc/sqrt(aa*bb) print *,i,' CF(%):',cf(i),' r(ak,bk):',rakbk(i) enddo c print*,' r(ak,bk):',(rakbk(i),i=1,md) return end c======================================================= subroutine svdcmp(a,m,n,mp,np,w,v) parameter (nmax=794) ! changed from 100 to 500 dimension a(mp,np),w(np),v(np,np),rv1(nmax) if(n.gt.nmax) then print*,'nmax,n::',nmax,n stop 'svdcmp' endif g=0.0 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 if (i.ne.n) then 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 endif 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 if (i.ne.m) then 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 endif 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=n,1,-1 l=i+1 g=w(i) if (i.lt.n) then do 33 j=l,n a(i,j)=0.0 33 continue endif if (g.ne.0.0) then g=1.0/g if (i.ne.n) then 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 endif 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) go to 2 if ((abs(w(nm))+anorm).eq.anorm) go to 1 41 continue 1 c=0.0 s=1.0 do 43 i=l,k f=s*rv1(i) if ((abs(f)+anorm).ne.anorm) then g=w(i) h=sqrt(f*f+g*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 endif 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 go to 3 endif ccc if (its.eq.30) pause 'no convergence in 30 iterations' if (its.eq.30) print*,'no convergence in 30 iterations' 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=sqrt(f*f+1.0) 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=sqrt(f*f+h*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 nm=1,n x=v(nm,j) z=v(nm,i) v(nm,j)= (x*c)+(z*s) v(nm,i)=-(x*s)+(z*c) 45 continue z=sqrt(f*f+h*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 nm=1,m y=a(nm,j) z=a(nm,i) a(nm,j)= (y*c)+(z*s) a(nm,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 c=================================================== c given singualr values w, singular vectors u and v c this subroutine sorts w into descending order,and c rearranges the columns of u and v correspondingly c=================================================== subroutine svdsrt(w,u,v,m,n,mp,np) dimension w(np),u(mp,np),v(np,np) do i=1,n-1 k=i p=w(i) do j=i+1,n if(w(j).ge.p)then k=j p=w(j) endif enddo if(k.ne.i)then w(k)=w(i) w(i)=p do j=1,n p=v(j,i) v(j,i)=v(j,k) v(j,k)=p enddo do j=1,m p=u(j,i) u(j,i)=u(j,k) u(j,k)=p enddo endif enddo return end