fortran使用MKL函数库中的gesvd求矩阵的SVD分解

下面的代码是使用MKl函数库中的gesvd进行矩阵的svd分解
运行环境: win10 + vs2013 + ivf2013
代码如下：
program svd
        use lapack95
        implicit none
        integer           :: i, m, n
        real, allocatable :: A(:,:), U(:,:), S(:,:), V(:,:), s1(:)
        
        write(*,'(g0)') '请输入矩阵的大小m,n: '
        read(*,*) m, n
        
        allocate( A(m,n), U(m,m), S(m,n), s1(min(m,n)), V(n,n) )

        call random_seed()
        call random_number(A)

        write(*,'(g0)') 'A is...'
        do i = 1, size(A,1)
                write(*,'(*(f11.4))') A(i,:)
        end do
        

        call gesvd( A, s1, U, V )  !// 这里返回的V是V的转置Vt
        
        write(*,'(g0)') 'U is...'
        do i = 1, size(U,1)
                write(*,'(*(f11.4))') U(i,:)
        end do
        
        write(*,'(g0)') 'S is...'
        S = 0.d0
        forall( i = 1:min(m,n) ) S(i,i) = s1(i)
        
        do i = 1, size(S,1)
                write(*,'(*(f11.4))') S(i,:)
        end do
        
        write(*,'(g0)') 'V is...'
        do i = 1, size(V,1)
                write(*,'(*(f11.4))') V(i,:)
        end do
        
        write(*,'(g0)') 'cheching, A is ...'
        A = matmul( matmul(U,S), V )  !// A = U*S*V'
        do i = 1, size(A,1)
                write(*,'(*(f11.4))') A(i,:)
        end do
        
        deallocate( A, U, S, V, s1 )
        
end program svd

运行结果如下：
请输入矩阵的大小m,n:
3 3
A is...
     0.4229     0.7888     0.6318
     0.4695     0.8313     0.3344
     0.7387     0.8524     0.5978
U is...
    -0.5571    -0.7084     0.4333
    -0.5130     0.7039     0.4913
    -0.6530     0.0514    -0.7556
S is...
     1.9437     0.0000     0.0000
     0.0000     0.2064     0.0000
     0.0000     0.0000     0.1796
V is...
    -0.4933    -0.7318    -0.4702
     0.3336     0.3401    -0.8793
    -0.8034     0.5906    -0.0764
cheching, A is ...
     0.4229     0.7888     0.6318
     0.4695     0.8313     0.3344
     0.7387     0.8524     0.5978

说明计算正确