1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8module m_svd
9public :: solve_with_svd
10CONTAINS
11subroutine solve_with_svd(ain,rhs,c,info,rcond,rank_out,sigma)
12  use precision, only: dp
13
14  real(dp), intent(in)  :: ain(:,:)
15  real(dp), intent(in)  :: rhs(:)
16  real(dp), intent(out) :: c(:)
17  integer, intent(out)  :: info
18  real(dp), intent(in), optional   :: rcond
19  integer, intent(out), optional  :: rank_out
20  real(dp), intent(out), optional :: sigma(:)
21
22  real(dp), allocatable :: a(:,:), s(:), work(:), b(:)
23
24  integer, parameter :: nb = 64
25  integer :: n, lwork
26
27  REAL(DP) RCOND_local, RNORM
28  INTEGER  ::    I, J, M, RANK, LDA
29
30  REAL(DP), external ::  DNRM2
31  EXTERNAL ::      DGELSS
32
33!
34  n = size(ain,dim=1)
35  lda = n
36  m = n
37  allocate(a(n,n))
38  a = ain
39  allocate(b(n),s(n))
40
41  lwork = 3*n+nb*(n+m)
42  allocate(work(lwork))
43
44  b(:) = rhs(:)
45
46!        Choose RCOND to reflect the relative accuracy of the input data
47
48  if (present(rcond)) then
49     rcond_local = rcond
50  else
51     RCOND_local = 1.0e-6_dp
52  endif
53  ! Singular values s_i < rcond*s_1 will be neglected for
54  ! the estimation of the rank.
55
56!        Solve the least squares problem min( norm2(b - Ax) ) for the x
57!        of minimum norm.
58
59         CALL DGELSS(M,N,1,A,LDA,B,M,S,RCOND_local,RANK,WORK,LWORK,INFO)
60!
61         IF (INFO.EQ.0) THEN
62
63            c(1:n) = b(1:n)
64
65            if (present(sigma)) then
66               sigma(1:n) = s(1:n)
67            endif
68            if (present(rank_out)) then
69               rank_out = rank
70            endif
71
72         END IF
73
74         deallocate(a,s,work,b)
75
76END subroutine solve_with_svd
77end module m_svd
78