! Subroutine Sort_svd (U,V,W,Nx,Ny) ! Sort svd values from hgihgest to lowest ! Integer :: q, i,j,n,m,k , Nx, Ny ! Real :: temp1,temp2(Nx),temp3(Ny) ! Real :: U(Nx,Ny), V(Ny,Ny), W(Ny) ! m=Nx; n=Ny ! q=min(m,n) ! do k=1,q ! do j=1,(q-k) ! 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 ! u(i,j)=u(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 ! u(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 ! End Subroutine Sort_svd !--------------------------------------------------------------------------------- SUBROUTINE svdcmp (A,U,V,W,Nx,Ny) !(a,m,n,mp,np,w,v) INTEGER :: m,mp,n,np , Nx,Ny INTEGER :: i,its,j,jj,k,l,nm REAL :: anorm,c,f,g,h,s,scale,x,y,z !, pythag !rv1(NMAX),pythag Real :: A(Nx,Ny), U(Nx,Ny), V(Ny,Ny), W(Ny), rv1(Ny) ! print*,"A-SVD=",A !M=Nt; N=Nx M=Nx; N=Ny g=0.0 scale=0.0 anorm=0.0 do 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 k=i,m scale=scale+abs(a(k,i)) enddo if(scale.ne.0.0)then do k=i,m a(k,i)=a(k,i)/scale s=s+a(k,i)*a(k,i) enddo f=a(i,i) g=-sign(sqrt(s),f) h=f*g-s a(i,i)=f-g do j=l,n s=0.0 do k=i,m s=s+a(k,i)*a(k,j) enddo f=s/h do k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo enddo do k=i,m a(k,i)=scale*a(k,i) enddo 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 k=l,n scale=scale+abs(a(i,k)) enddo if(scale.ne.0.0)then do k=l,n a(i,k)=a(i,k)/scale s=s+a(i,k)*a(i,k) enddo f=a(i,l) g=-sign(sqrt(s),f) h=f*g-s a(i,l)=f-g do k=l,n rv1(k)=a(i,k)/h enddo do j=l,m s=0.0 do k=l,n s=s+a(j,k)*a(i,k) enddo do k=l,n a(j,k)=a(j,k)+s*rv1(k) enddo enddo do k=l,n a(i,k)=scale*a(i,k) enddo endif endif anorm=max(anorm,(abs(w(i))+abs(rv1(i)))) enddo do i=n,1,-1 if(i.lt.n)then if(g.ne.0.0)then do j=l,n v(j,i)=(a(i,j)/a(i,l))/g enddo do j=l,n s=0.0 do k=l,n s=s+a(i,k)*v(k,j) enddo do k=l,n v(k,j)=v(k,j)+s*v(k,i) enddo enddo endif do j=l,n v(i,j)=0.0 v(j,i)=0.0 enddo endif v(i,i)=1.0 g=rv1(i) l=i enddo do i=min(m,n),1,-1 l=i+1 g=w(i) do j=l,n a(i,j)=0.0 enddo if(g.ne.0.0)then g=1.0/g do j=l,n s=0.0 do k=l,m s=s+a(k,i)*a(k,j) enddo f=(s/a(i,i))*g do k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo enddo do j=i,m a(j,i)=a(j,i)*g enddo else do j= i,m a(j,i)=0.0 enddo endif a(i,i)=a(i,i)+1.0 enddo do k=n,1,-1 do its=1,30 do 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 enddo 1 c=0.0 s=1.0 do 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 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) enddo enddo 2 z=w(k) if(l.eq.k)then if(z.lt.0.0)then w(k)=-z do j=1,n v(j,k)=-v(j,k) enddo 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 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 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) enddo 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 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) enddo enddo rv1(l)=0.0 rv1(k)=f w(k)=x enddo 3 continue enddo U(:,:)=A(:,:) !print*,"U-SVD=",U END SUBROUTINE svdcmp !--------------------------------------------------------------------------------- Real 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 END FUNCTION pythag !---------------------------------------------------------------------------------