subroutine svd_s(A, LDA, U, LDU, D, Vt, LDVt, m, n)
  implicit none
  BEGIN_DOC
  ! !!!
  ! DGESVD computes the singular value decomposition (SVD) of a real
  ! M-by-N matrix A, optionally computing the left and/or right singular
  ! vectors. The SVD is written:
  !      A = U * SIGMA * transpose(V)
  ! where SIGMA is an M-by-N matrix which is zero except for its
  ! min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  ! V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  ! are the singular values of A; they are real and non-negative, and
  ! are returned in descending order.  The first min(m,n) columns of
  ! U and V are the left and right singular vectors of A.
  !
  ! Note that the routine returns V**T, not V.
  ! !!!
  END_DOC

  integer, intent(in)             :: LDA, LDU, LDVt, m, n
  double precision, intent(in)    :: A(LDA,n)
  double precision, intent(out)   :: U(LDU,min(m,n)), Vt(LDVt,n), D(min(m,n))
  double precision,allocatable    :: work(:), A_tmp(:,:)
  integer                         :: info, lwork, i, j, k


   allocate (A_tmp(LDA,n))
   do k=1,n
     do i=1,m
       A_tmp(i,k) = A(i,k) + 1d-16
     enddo
   enddo
 
   ! Find optimal size for temp arrays
   allocate(work(1))
   lwork = -1
   call dgesvd('A', '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)), 5*MIN(M,N))
   deallocate(work)

   allocate(work(lwork))
   call dgesvd('A','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 in svds 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, LDVt
       if (dabs(Vt(i,j)) < 1.d-14) Vt(i,j) = 0.d0
     enddo
   enddo
 
end






subroutine my_ortho_qr(A,LDA,m,n)
  implicit none
  BEGIN_DOC
  ! Orthogonalization using Q.R factorization
  !
  ! A : matrix to orthogonalize
  !
  ! LDA : leftmost dimension of A
  !
  ! m : Number of rows of A
  !
  ! n : Number of columns of A
  !
  END_DOC
  integer, intent(in)             :: m, n, LDA
  double precision, intent(inout) :: A(LDA,n)
  integer                         :: LWORK, INFO, nTAU, ii, jj
  double precision, allocatable   :: TAU(:), WORK(:)
  double precision                :: Adorgqr(LDA,n)

  allocate (TAU(min(m,n)), WORK(1))

  LWORK=-1
  call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
  LWORK=max(n,int(WORK(1)))

  deallocate(WORK)
  allocate(WORK(LWORK))
  call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
  if(INFO.ne.0 ) then
   print*, 'problem in DGEQRF'
  endif

  nTAU = size(TAU)
  do jj = 1, n
    do ii = 1, LDA
      Adorgqr(ii,jj) = A(ii,jj)
    enddo
  enddo

  LWORK=-1
  call dorgqr(m, n, nTAU, Adorgqr, LDA, TAU, WORK, LWORK, INFO)
  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
  LWORK=max(n,int(WORK(1)))

  deallocate(WORK)
  allocate(WORK(LWORK))
  call dorgqr(m, n, nTAU, A, LDA, TAU, WORK, LWORK, INFO)
  if(INFO.ne.0 ) then
   print*, 'problem in DORGQR'
  endif


  deallocate(WORK,TAU)
end