Merge pull request #96 from nekStab/stdlib_svd

Switch to `stdlib_svd` + minor typo corrections.
nekStab · Jun 15, 2024 · beffa01 · beffa01
2 parents 86ec306 + fd546e9
commit beffa01
Show file tree

Hide file tree

Showing 6 changed files with 58 additions and 262 deletions.
diff --git a/src/Constants.f90 b/src/Constants.f90
@@ -11,7 +11,7 @@ module lightkrylov_constants
 
     integer , parameter, public :: dp = selected_real_kind(15, 307)
     !! Definition of the double precision data type.
-    real(dp), parameter, public :: atol_dp = 10.0_sp ** -precision(1.0_dp)
+    real(dp), parameter, public :: atol_dp = 10.0_dp ** -precision(1.0_dp)
     !! Definition of the absolute tolerance for double precision computations.
     real(dp), parameter, public :: rtol_dp = sqrt(atol_dp)
     !! Definition of the relative tolerance for double precision computations.

diff --git a/src/IterativeSolvers.f90 b/src/IterativeSolvers.f90
@@ -5,7 +5,7 @@ module lightkrylov_IterativeSolvers
     use stdlib_sorting, only: sort_index
     use stdlib_optval, only: optval
     use stdlib_io_npy, only: save_npy
-    use stdlib_linalg, only: lstsq
+    use stdlib_linalg, only: lstsq, svd
     use stdlib_stats, only: median
 
     use lightkrylov_constants
@@ -1230,17 +1230,19 @@ subroutine svds_rsp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_rsp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_sp ; umat = 0.0_sp ; vmat = 0.0_sp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_rsp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = transpose(vmat(:k, :k))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_rsp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1325,17 +1327,19 @@ subroutine svds_rdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_rdp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_dp ; umat = 0.0_dp ; vmat = 0.0_dp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_rdp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = transpose(vmat(:k, :k))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_rdp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1420,17 +1424,19 @@ subroutine svds_csp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_csp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_sp ; umat = 0.0_sp ; vmat = 0.0_sp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_csp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_csp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1515,17 +1521,19 @@ subroutine svds_cdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_cdp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_dp ; umat = 0.0_dp ; vmat = 0.0_dp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_cdp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_cdp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------

diff --git a/src/IterativeSolvers.fypp b/src/IterativeSolvers.fypp
@@ -7,7 +7,7 @@ module lightkrylov_IterativeSolvers
     use stdlib_sorting, only: sort_index
     use stdlib_optval, only: optval
     use stdlib_io_npy, only: save_npy
-    use stdlib_linalg, only: lstsq
+    use stdlib_linalg, only: lstsq, svd
     use stdlib_stats, only: median
 
     use lightkrylov_constants
@@ -466,17 +466,23 @@ contains
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_${type[0]}$${kind}$')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_${kind}$ ; umat = 0.0_${kind}$ ; vmat = 0.0_${kind}$
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_${type[0]}$${kind}$(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                #:if type[0] == "r"
+                vmat(:k, :k) = transpose(vmat(:k, :k))
+                #:else
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
+                #:endif
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_${type[0]}$${kind}$(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------

diff --git a/src/Logger.f90 b/src/Logger.f90
@@ -139,20 +139,7 @@ subroutine check_info(info, origin, module, procedure)
                call logger%log_error(origin, module=module, procedure=procedure, stat=info, errmsg=trim(msg))
                ierr = -1
             end if
-         else if (trim(to_lower(origin)) == 'gesvd') then
-            ! GESVD
-            if (info < 0) then
-               write(msg, *) "The ", -info, "-th argument has illegal value."
-               call logger%log_error(origin, module=module, procedure=procedure, stat=info, errmsg=trim(msg))
-               ierr = -1
-            else
-               write(msg, *) "The QR alg. did not converge. ", info, &
-                           & "superdiagonals of an intermediate bidiagonal form B ", &
-                           & "did not converge to zero."
-               call logger%log_error(origin, module=module, procedure=procedure, stat=info, errmsg=trim(msg))
-               ierr = -1
-            end if
-         !
+        !
          !   LightKrylov_Utils
          !
          else if (trim(to_lower(origin)) == 'sqrtm') then

diff --git a/src/Utils.f90 b/src/Utils.f90
@@ -29,8 +29,6 @@ module lightkrylov_utils
     public :: assert_shape
     ! Compute B = inv(A) in-place for dense matrices.
     public :: inv
-    ! Compute USV^T = svd(A) for dense matrices.
-    public :: svd
     ! Compute AX = XD for general dense matrices.
     public :: eig
     ! Compute AX = XD for symmetric/hermitian matrices.
@@ -73,13 +71,6 @@ module lightkrylov_utils
         module procedure inv_cdp
     end interface
 
-    interface svd
-        module procedure svd_rsp
-        module procedure svd_rdp
-        module procedure svd_csp
-        module procedure svd_cdp
-    end interface
-
     interface eig
         module procedure eig_rsp
         module procedure eig_rdp
@@ -365,41 +356,6 @@ subroutine inv_rsp(A)
         return
     end subroutine inv_rsp
 
-    subroutine svd_rsp(A, U, S, V)
-        !! Singular value decomposition of a dense matrix.
-        real(sp), intent(in) :: A(:, :)
-        !! Matrix to be factorized.
-        real(sp), intent(out) :: U(:, :)
-        !! Left singular vectors.
-        real(sp), intent(out) :: S(:)
-        !! Singular values.
-        real(sp), intent(out) :: V(:, :)
-        !! Right singular vectors.
-
-        ! Internal variables.
-        character :: jobu = "S", jobvt = "S"
-        integer :: m, n, lda, ldu, ldvt, lwork, info
-        real(sp), allocatable :: work(:)
-        real(sp) :: A_tilde(size(A, 1), size(A, 2)), Vt(min(size(A, 1), size(A, 2)), size(A, 2))
-
-        ! Setup variables.
-        m = size(A, 1) ; n = size(A, 2)
-        lda = m ; ldu = m ; ldvt = n
-        lwork = max(1, 3*min(m, n) + max(m, n), 5*min(m, n)) ; allocate(work(lwork))
-
-        ! Shape assertion.
-        call assert_shape(U, [m, m], "svd", "U")
-        call assert_shape(V, [n, n], "svd", "V")
-
-        ! SVD computation.
-        A_tilde = A
-        call gesvd(jobu, jobvt, m, n, A_tilde, lda, S, U, ldu, Vt, ldvt, work, lwork, info)
-        v = transpose(vt)
-        call check_info(info, 'GESVD', module=this_module, procedure='svd_rsp')
-
-        return
-    end subroutine svd_rsp
-
     subroutine eig_rsp(A, vecs, vals)
         !! Eigenvalue decomposition of a dense matrix using LAPACK.
         real(sp), intent(in) :: A(:, :)
@@ -592,41 +548,6 @@ subroutine inv_rdp(A)
         return
     end subroutine inv_rdp
 
-    subroutine svd_rdp(A, U, S, V)
-        !! Singular value decomposition of a dense matrix.
-        real(dp), intent(in) :: A(:, :)
-        !! Matrix to be factorized.
-        real(dp), intent(out) :: U(:, :)
-        !! Left singular vectors.
-        real(dp), intent(out) :: S(:)
-        !! Singular values.
-        real(dp), intent(out) :: V(:, :)
-        !! Right singular vectors.
-
-        ! Internal variables.
-        character :: jobu = "S", jobvt = "S"
-        integer :: m, n, lda, ldu, ldvt, lwork, info
-        real(dp), allocatable :: work(:)
-        real(dp) :: A_tilde(size(A, 1), size(A, 2)), Vt(min(size(A, 1), size(A, 2)), size(A, 2))
-
-        ! Setup variables.
-        m = size(A, 1) ; n = size(A, 2)
-        lda = m ; ldu = m ; ldvt = n
-        lwork = max(1, 3*min(m, n) + max(m, n), 5*min(m, n)) ; allocate(work(lwork))
-
-        ! Shape assertion.
-        call assert_shape(U, [m, m], "svd", "U")
-        call assert_shape(V, [n, n], "svd", "V")
-
-        ! SVD computation.
-        A_tilde = A
-        call gesvd(jobu, jobvt, m, n, A_tilde, lda, S, U, ldu, Vt, ldvt, work, lwork, info)
-        v = transpose(vt)
-        call check_info(info, 'GESVD', module=this_module, procedure='svd_rdp')
-
-        return
-    end subroutine svd_rdp
-
     subroutine eig_rdp(A, vecs, vals)
         !! Eigenvalue decomposition of a dense matrix using LAPACK.
         real(dp), intent(in) :: A(:, :)
@@ -819,43 +740,6 @@ subroutine inv_csp(A)
         return
     end subroutine inv_csp
 
-    subroutine svd_csp(A, U, S, V)
-        !! Singular value decomposition of a dense matrix.
-        complex(sp), intent(in) :: A(:, :)
-        !! Matrix to be factorized.
-        complex(sp), intent(out) :: U(:, :)
-        !! Left singular vectors.
-        real(sp), intent(out) :: S(:)
-        !! Singular values.
-        complex(sp), intent(out) :: V(:, :)
-        !! Right singular vectors.
-
-        ! Internal variables.
-        character :: jobu = "S", jobvt = "S"
-        integer :: m, n, lda, ldu, ldvt, lwork, info
-        complex(sp), allocatable :: work(:)
-        complex(sp) :: A_tilde(size(A, 1), size(A, 2)), Vt(min(size(A, 1), size(A, 2)), size(A, 2))
-        real(sp), allocatable :: rwork(:)
-
-        ! Setup variables.
-        m = size(A, 1) ; n = size(A, 2)
-        lda = m ; ldu = m ; ldvt = n
-        lwork = max(1, 3*min(m, n) + max(m, n), 5*min(m, n)) ; allocate(work(lwork))
-
-        ! Shape assertion.
-        call assert_shape(U, [m, m], "svd", "U")
-        call assert_shape(V, [n, n], "svd", "V")
-
-        ! SVD computation.
-        A_tilde = A
-        allocate(rwork(5*min(m, n)))
-        call gesvd(jobu, jobvt, m, n, A_tilde, lda, S, U, ldu, Vt, ldvt, work, lwork, rwork, info)
-        v = transpose(conjg(vt))
-        call check_info(info, 'GESVD', module=this_module, procedure='svd_csp')
-
-        return
-    end subroutine svd_csp
-
     subroutine eig_csp(A, vecs, vals)
         !! Eigenvalue decomposition of a dense matrix using LAPACK.
         complex(sp), intent(in) :: A(:, :)
@@ -1043,43 +927,6 @@ subroutine inv_cdp(A)
         return
     end subroutine inv_cdp
 
-    subroutine svd_cdp(A, U, S, V)
-        !! Singular value decomposition of a dense matrix.
-        complex(dp), intent(in) :: A(:, :)
-        !! Matrix to be factorized.
-        complex(dp), intent(out) :: U(:, :)
-        !! Left singular vectors.
-        real(dp), intent(out) :: S(:)
-        !! Singular values.
-        complex(dp), intent(out) :: V(:, :)
-        !! Right singular vectors.
-
-        ! Internal variables.
-        character :: jobu = "S", jobvt = "S"
-        integer :: m, n, lda, ldu, ldvt, lwork, info
-        complex(dp), allocatable :: work(:)
-        complex(dp) :: A_tilde(size(A, 1), size(A, 2)), Vt(min(size(A, 1), size(A, 2)), size(A, 2))
-        real(dp), allocatable :: rwork(:)
-
-        ! Setup variables.
-        m = size(A, 1) ; n = size(A, 2)
-        lda = m ; ldu = m ; ldvt = n
-        lwork = max(1, 3*min(m, n) + max(m, n), 5*min(m, n)) ; allocate(work(lwork))
-
-        ! Shape assertion.
-        call assert_shape(U, [m, m], "svd", "U")
-        call assert_shape(V, [n, n], "svd", "V")
-
-        ! SVD computation.
-        A_tilde = A
-        allocate(rwork(5*min(m, n)))
-        call gesvd(jobu, jobvt, m, n, A_tilde, lda, S, U, ldu, Vt, ldvt, work, lwork, rwork, info)
-        v = transpose(conjg(vt))
-        call check_info(info, 'GESVD', module=this_module, procedure='svd_cdp')
-
-        return
-    end subroutine svd_cdp
-
     subroutine eig_cdp(A, vecs, vals)
         !! Eigenvalue decomposition of a dense matrix using LAPACK.
         complex(dp), intent(in) :: A(:, :)