c Fortran program to perform singular value decomposition on c a space-time array c c Input file is the space-time array in binary format (AT array) c c Output files are: c V array (spacedim x spacedim): Normalized vectors (var=1) c U array (timedim x timedim): Normalized vectors (var=1) c Singular values (sigma) (maxdim) c c c Specify space and time dimensions c c maxdim is the larger of timedim or spacedim c mindim is the smaller of timedim or spacedim c parameter(timedim=10) parameter(spacedim=20) parameter(mindim=10) parameter(maxdim=20) c c Specify variables and arrays c c Integer variables c m = timedim c n = spacedim c q = smaller of timedim or spacedim (mindim) c i,j,pass = counters c c Arrays c c at = original data in format space-time d a = original data in format time-space c v = V array (**Note it's not VT) c w = Sigma array: Singular values c c temp arrays used for data sorting integer m,n,q,i,j,pass real at(spacedim,timedim) real a(timedim,spacedim),v(spacedim,spacedim) real w(maxdim) real temp1,temp2(timedim),temp3(spacedim) c Input and output files open (10,file='Aarray.bin', + form='unformatted', access='direct', + recl=4*spacedim*timedim) open (11,file='Uarray.bin', + form='unformatted', access='direct', + recl=4*timedim*timedim) open (12,file='Varray.bin', + form='unformatted', access='direct', + recl=4*spacedim*spacedim) open (13, file='Sigmavalues.asci', + form='formatted') c Read in data (AT array) read(10,rec=1) ((at(j,i),j=1,spacedim),i=1,timedim) n=spacedim m=timedim c Transpose data so it is in the form time-space do j=1,n do i=1,m a(i,j)=at(j,i) end do end do c Call SVDcmp subroutine print*, 'To svdcmp' call svdcmp(a,m,n,m,n,w,v) print*, 'Exit svdcmp' c Sort the data from highest to lowest eigenvalue c Null space excluded q=mindim do pass=1,q do j=1,(q-pass) if(w(j).le.w(j+1)) then temp1=w(j) do i=1,m temp2(i)=a(i,j) end do do i=1,n temp3(i)=v(i,j) end do w(j)=w(j+1) do i=1,m a(i,j)=a(i,j+1) end do do i=1,n v(i,j)=v(i,j+1) end do w(j+1)=temp1 do i=1,m a(i,j+1)=temp2(i) end do do i=1,n v(i,j+1)=temp3(i) end do end if end do end do c Output U and V arrays, Singular values write(11,rec=1) ((a(i,j),i=1,timedim),j=1,timedim) write(12,rec=1) ((v(i,j),i=1,spacedim),j=1,spacedim) do i=1,q write(13,*) w(i) end do print*, 'SVD Analysis complete' end c SVDcmp subroutine c Given a matrix (1:m,1:n) with physical dimensions mp by np, c this routine computes its singular value decomposition, c A = U W VT. The matrix U replaces a on output. The diagonal c matrix of singular values W is output as a vector w(1:n) c The matrix V (not the transpose VT) is the output as v(1:n,1:n) c SUBROUTINE svdcmp(a,m,n,mp,np,w,v) INTEGER m,mp,n,np,NMAX REAL a(mp,np),v(np,np),w(np) PARAMETER (NMAX=10000) CU USES pythag INTEGER i,its,j,jj,k,l,nm REAL anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythag 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 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,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=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 C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3. C C FUNCTION pythag(a,b) REAL a,b,pythag REAL absa,absb absa=abs(a) absb=abs(b) if(absa.gt.absb)then pythag=absa*sqrt(1.+(absb/absa)**2) else if(absb.eq.0.)then pythag=0. else pythag=absb*sqrt(1.+(absa/absb)**2) endif endif return END C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.