diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 651376e2..5618a6c0 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -13,6 +13,7 @@ qp_basis_clean.native qp_create_ezfio_from_xyz qp_create_ezfio_from_xyz.native qp_edit +qp_edit.ml qp_edit.native qp_print qp_print.native diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index 737f03f7..9b3fa7e9 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -59,6 +59,8 @@ enddo enddo !$OMP END PARALLEL DO + + ! Check for linear dependencies in the basis set END_PROVIDER diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index a7462f94..c5e66006 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -1,3 +1,48 @@ +subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) + implicit none + BEGIN_DOC + ! Compute A = U.D.Vt + ! + ! LDx : leftmost dimension of x + ! + ! Dimsneion of A is m x n + ! + END_DOC + + integer, intent(in) :: LDA, LDU, LDVt, m, n + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU,n) + double precision,intent(out) :: Vt(LDVt,n) + double precision,intent(out) :: D(n) + double precision,allocatable :: work(:) + integer :: info, lwork, i, j, k + + double precision,allocatable :: A_tmp(:,:) + allocate (A_tmp(LDA,n)) + A_tmp = A + + ! Find optimal size for temp arrays + allocate(work(1)) + lwork = -1 + call dgesvd('A','A', n, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + lwork = work(1) + deallocate(work) + + allocate(work(lwork)) + call dgesvd('A','A', n, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + deallocate(work,A_tmp) + + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + +end + + + subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) implicit none BEGIN_DOC @@ -29,25 +74,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D, work integer :: info, lwork, i, j, k - double precision,allocatable :: overlap_tmp(:,:) - allocate (overlap_tmp(lda,n)) - overlap_tmp = overlap - - allocate(work(1)) - lwork = -1 - call dgesvd('A','A', n, n, overlap_tmp, lda, & - D, U, ldc, Vt, lda, work, lwork, info) - lwork = work(1) - deallocate(work) - allocate(work(lwork)) - call dgesvd('A','A', n, n, overlap_tmp, lda, & - D, U, ldc, Vt, lda, work, lwork, info) - deallocate(work,overlap_tmp) - if (info /= 0) then - print *, info, ': SVD failed' - stop - endif - + call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S_half,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j,k)