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
+
+