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