From 4a6b95aac7b423c78ccec5695810c480acee9030 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 4 Sep 2024 08:22:03 +0200 Subject: [PATCH] added local svd --- src/utils/non_sym_diag.f90 | 61 +++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 index 0f44099..1659533 100644 --- a/src/utils/non_sym_diag.f90 +++ b/src/utils/non_sym_diag.f90 @@ -512,7 +512,7 @@ subroutine impose_biorthog_svd(n, m, L, R) allocate(U(m,m), Vt(m,m), D(m)) - call svd(S, m, U, m, D, Vt, m, m, m) + call svd_local(S, m, U, m, D, Vt, m, m, m) deallocate(S) @@ -569,3 +569,62 @@ end ! --- +subroutine svd_local(A, LDA, U, LDU, D, Vt, LDVt, m, n) + + implicit none + + ! Compute A = U D Vt + + integer, intent(in) :: LDA, LDU, LDVt, m, n + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU, min(m, n)) + double precision, intent(out) :: Vt(LDVt,n) + double precision, intent(out) :: D(min(m, n)) + + integer :: info, lwork, i, j, k + double precision, allocatable :: work(:) + + double precision,allocatable :: A_tmp(:,:) + + allocate(A_tmp(LDA,n)) + do k = 1, n + do i = 1, m + A_tmp(i,k) = A(i,k) + enddo + enddo + + ! Find optimal size for temp arrays + allocate(work(1)) + lwork = -1 + call dgesvd('S','S', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 + lwork = max(int(work(1)), 10*MIN(M,N)) + deallocate(work) + + allocate(work(lwork)) + call dgesvd('S','S', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + + deallocate(A_tmp, work) + + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + + do j = 1, min(m, n) + do i = 1, m + if (dabs(U(i,j)) < 1.d-14) U(i,j) = 0.d0 + enddo + enddo + + do j = 1, n + do i = 1, n + if(dabs(Vt(i,j)) < 1.d-14) Vt(i,j) = 0.d0 + enddo + enddo + +end + +