mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-10-05 15:56:00 +02:00
128 lines
3.2 KiB
FortranFixed
128 lines
3.2 KiB
FortranFixed
|
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
|