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 ! ! Dimension 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,min(m,n)) double precision,intent(out) :: Vt(LDVt,n) double precision,intent(out) :: D(min(m,n)) double precision,allocatable :: work(:) integer :: info, lwork, i, j, k double precision,allocatable :: A_tmp(:,:) allocate (A_tmp(LDA,n)) do k=1,n do i=1,m A_tmp(i,k) = A(i,k) enddo enddo ! Find optimal size for temp arrays allocate(work(1)) lwork = -1 call dgesvd('S','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)), 10*MIN(M,N)) deallocate(work) allocate(work(lwork)) call dgesvd('S','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 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,n if (dabs(Vt(i,j)) < 1.d-14) Vt(i,j) = 0.d0 enddo enddo end subroutine svd_symm(A,LDA,U,LDU,D,Vt,LDVt,m,n) implicit none BEGIN_DOC ! Compute A = U.D.Vt ! ! LDx : leftmost dimension of x ! ! Dimension 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,min(m,n)) double precision,intent(out) :: Vt(LDVt,n) double precision,intent(out) :: D(min(m,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', 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','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) deallocate(A_tmp,work) if (info /= 0) then print *, info, ': SVD failed' stop endif ! Iterative refinement ! -------------------- ! https://doi.org/10.1016/j.cam.2019.112512 integer :: iter double precision,allocatable :: R(:,:), S(:,:), T(:,:) double precision,allocatable :: sigma(:), F(:,:), G(:,:) double precision :: alpha, beta, x, thresh allocate (R(m,m), S(n,n), T(m,n), sigma(m), F(m,m), G(n,n), A_tmp(m,n)) sigma = 0.d0 R = 0.d0 S = 0.d0 T = 0.d0 F = 0.d0 G = 0.d0 thresh = 1.d-8 call restore_symmetry(m,m,U,size(U,1),thresh) call restore_symmetry(n,n,Vt,size(Vt,1),thresh) do iter=1,4 do k=1,n A_tmp(1:m,k) = D(k) * U(1:m,k) enddo call dgemm('N','N',m,n,n,1.d0,A_tmp,size(A_tmp,1), & Vt,size(Vt,1),0.d0,R,size(R,1)) print *, maxval(dabs(R(1:m,1:n) - A(1:m,1:n))) call dgemm('T','N',m,m,m,-1.d0,U,size(U,1), & U,size(U,1),0.d0,R,size(R,1)) do i=1,m R(i,i) = R(i,i) + 1.d0 enddo call dgemm('N','T',n,n,n,-1.d0,Vt,size(Vt,1), & Vt,size(Vt,1),0.d0,S,size(S,1)) do i=1,n S(i,i) = S(i,i) + 1.d0 enddo call dgemm('T','N',m,n,m,1.d0,U,size(U,1), & A,size(A,1),0.d0,A_tmp,size(A_tmp,1)) call dgemm('N','T',m,n,n,1.d0,A_tmp,size(A_tmp,1), & Vt,size(Vt,1),0.d0,T,size(T,1)) do i=1,n sigma(i) = T(i,i)/(1.d0 - (R(i,i)+S(i,i))*0.5d0) F(i,i) = 0.5d0*R(i,i) G(i,i) = 0.5d0*S(i,i) enddo do j=1,n do i=1,n if (i == j) cycle alpha = T(i,j) + sigma(j) * R(i,j) beta = T(i,j) + sigma(j) * S(i,j) x = 1.d0 / (sigma(j)*sigma(j) - sigma(i)*sigma(i)) F(i,j) = (alpha * sigma(j) + beta * sigma(i)) * x G(i,j) = (alpha * sigma(i) + beta * sigma(j)) * x enddo enddo do i=1,n x = 1.d0/sigma(i) do j=n+1,m F(i,j) = -T(j,i) * x enddo enddo do i=n+1,m do j=1,n F(i,j) = R(i,j) - F(j,i) enddo enddo do i=n+1,m do j=n+1,m F(i,j) = R(i,j)*0.5d0 enddo enddo D(1:min(n,m)) = sigma(1:min(n,m)) call dgemm('N','N',m,m,m,1.d0,U,size(U,1),F,size(F,1), & 0.d0, A_tmp, size(A_tmp,1)) do j=1,m do i=1,m U(i,j) = U(i,j) + A_tmp(i,j) enddo enddo call dgemm('T','N',n,n,n,1.d0,G,size(G,1),Vt,size(Vt,1), & 0.d0, A_tmp, size(A_tmp,1)) do j=1,n do i=1,n Vt(i,j) = Vt(i,j) + A_tmp(i,j) enddo enddo thresh = 0.01d0 * thresh call restore_symmetry(m,m,U,size(U,1),thresh) call restore_symmetry(n,n,Vt,size(Vt,1),thresh) enddo deallocate(A_tmp,R,S,F,G,sigma) end subroutine eigSVD(A,LDA,U,LDU,D,Vt,LDVt,m,n) implicit none BEGIN_DOC ! Algorithm 3 of https://arxiv.org/pdf/1810.06860.pdf ! ! A(m,n) = U(m,n) D(n) Vt(n,n) with m>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) integer :: i,j,k if (m= 0.d0) then exit endif D(k) = dsqrt(-D(k)) call dscal(n, 1.d0/D(k), V(1,k), 1) enddo D(k:n) = 0.d0 k=k-1 call dgemm('N','N',m,n,k,1.d0,A,size(A,1),V,size(V,1),0.d0,U,size(U,1)) end subroutine randomized_svd(A,LDA,U,LDU,D,Vt,LDVt,m,n,q,r) implicit none include 'constants.include.F' BEGIN_DOC ! Randomized SVD: rank r, q power iterations ! ! 1. Sample column space of A with P: Z = A.P where P is random with r+p columns. ! ! 2. Power iterations : Z <- X . (Xt.Z) ! ! 3. Z = Q.R ! ! 4. Compute SVD on projected Qt.X = U' . S. Vt ! ! 5. U = Q U' END_DOC integer, intent(in) :: LDA, LDU, LDVt, m, n, q, r double precision, intent(in) :: A(LDA,n) double precision, intent(out) :: U(LDU,r) double precision,intent(out) :: Vt(LDVt,r) double precision,intent(out) :: D(r) integer :: i, j, k double precision,allocatable :: Z(:,:), P(:,:), Y(:,:), UY(:,:) double precision :: r1,r2 allocate(P(n,r), Z(m,r)) ! P is a normal random matrix (n,r) do k=1,r do i=1,n call random_number(r1) call random_number(r2) r1 = dsqrt(-2.d0*dlog(r1)) r2 = dtwo_pi*r2 P(i,k) = r1*dcos(r2) enddo enddo ! Z(m,r) = A(m,n).P(n,r) call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1)) ! Power iterations do k=1,q ! P(n,r) = At(n,m).Z(m,r) call dgemm('T','N',n,r,m,1.d0,A,size(A,1),Z,size(Z,1),0.d0,P,size(P,1)) ! Z(m,r) = A(m,n).P(n,r) call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1)) enddo deallocate(P) ! QR factorization of Z call ortho_svd(Z,size(Z,1),m,r) allocate(Y(r,n), UY(r,r)) ! Y(r,n) = Zt(r,m).A(m,n) call dgemm('T','N',r,n,m,1.d0,Z,size(Z,1),A,size(A,1),0.d0,Y,size(Y,1)) ! SVD of Y call svd(Y,size(Y,1),UY,size(UY,1),D,Vt,size(Vt,1),r,n) deallocate(Y) ! U(m,r) = Z(m,r).UY(r,r) call dgemm('N','N',m,r,r,1.d0,Z,size(Z,1),UY,size(UY,1),0.d0,U,size(U,1)) deallocate(UY,Z) end subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) implicit none BEGIN_DOC ! Compute A = U.D.Vt ! ! LDx : leftmost dimension of x ! ! Dimension of A is m x n ! A,U,Vt are complex*16 ! D is double precision END_DOC integer, intent(in) :: LDA, LDU, LDVt, m, n complex*16, intent(in) :: A(LDA,n) complex*16, intent(out) :: U(LDU,m) complex*16, intent(out) :: Vt(LDVt,n) double precision,intent(out) :: D(min(m,n)) complex*16,allocatable :: work(:) double precision,allocatable :: rwork(:) integer :: info, lwork, i, j, k, lrwork complex*16,allocatable :: A_tmp(:,:) allocate (A_tmp(LDA,n)) A_tmp = A lrwork = 5*min(m,n) ! Find optimal size for temp arrays allocate(work(1),rwork(lrwork)) lwork = -1 call zgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, rwork, info) lwork = int(work(1)) deallocate(work) allocate(work(lwork)) call zgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, rwork, info) deallocate(work,rwork,A_tmp) if (info /= 0) then print *, info, ': SVD failed' stop endif end subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. ! ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! ! N : Overlap matrix is NxN (array is (LDA,N) ) ! ! C : Coefficients of the vectors to orthogonalize. On exit, ! orthogonal vectors ! ! LDC : leftmost dimension of C ! ! m : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC integer, intent(in) :: lda, ldc, n integer, intent(out) :: m complex*16, intent(in) :: overlap(lda,n) double precision, intent(in) :: cutoff complex*16, intent(inout) :: C(ldc,n) complex*16, allocatable :: U(:,:) complex*16, allocatable :: Vt(:,:) double precision, allocatable :: D(:) complex*16, allocatable :: S(:,:) !DIR$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j double precision :: local_cutoff if (n < 2) then return endif allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) D(:) = dsqrt(D(:)) local_cutoff = dsqrt(cutoff)*D(1) ! such that D(i)/D(1) > dsqrt(cutoff) is kept m=n do i=1,n if ( D(i) >= local_cutoff ) then D(i) = 1.d0/D(i) else m = i-1 print *, 'Removed Linear dependencies below:', local_cutoff exit endif enddo do i=m+1,n D(i) = 0.d0 enddo do j=1,n do i=1,n S(i,j) = U(i,j)*D(j) enddo enddo do j=1,n do i=1,n U(i,j) = C(i,j) enddo enddo call zgemm('N','N',n,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) deallocate (U, Vt, D, S) end subroutine ortho_qr_complex(A,LDA,m,n) implicit none BEGIN_DOC ! Orthogonalization using Q.R factorization ! ! A : matrix to orthogonalize ! ! LDA : leftmost dimension of A ! ! n : Number of rows of A ! ! m : Number of columns of A ! END_DOC integer, intent(in) :: m,n, LDA complex*16, intent(inout) :: A(LDA,n) integer :: lwork, info integer, allocatable :: jpvt(:) complex*16, allocatable :: tau(:), work(:) allocate (jpvt(n), tau(n), work(1)) LWORK=-1 call zgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO ) LWORK=2*int(WORK(1)) deallocate(WORK) allocate(WORK(LWORK)) call zgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO ) call zungqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO) deallocate(WORK,jpvt,tau) end subroutine ortho_qr_unblocked_complex(A,LDA,m,n) implicit none BEGIN_DOC ! Orthogonalization using Q.R factorization ! ! A : matrix to orthogonalize ! ! LDA : leftmost dimension of A ! ! n : Number of rows of A ! ! m : Number of columns of A ! END_DOC integer, intent(in) :: m,n, LDA double precision, intent(inout) :: A(LDA,n) integer :: info integer, allocatable :: jpvt(:) double precision, allocatable :: tau(:), work(:) print *, irp_here, ': TO DO' stop -1 ! allocate (jpvt(n), tau(n), work(n)) ! call dgeqr2( m, n, A, LDA, TAU, WORK, INFO ) ! call dorg2r(m, n, n, A, LDA, tau, WORK, INFO) ! deallocate(WORK,jpvt,tau) end subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.S^-1/2 orthogonalization. ! ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! ! N : Overlap matrix is NxN (array is (LDA,N) ) ! ! C : Coefficients of the vectors to orthogonalize. On exit, ! orthogonal vectors ! ! LDC : leftmost dimension of C ! ! M : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC integer, intent(in) :: LDA, ldc, n, m complex*16, intent(in) :: overlap(lda,n) complex*16, intent(inout) :: C(ldc,n) complex*16, allocatable :: U(:,:) complex*16, allocatable :: Vt(:,:) double precision, allocatable :: D(:) complex*16, allocatable :: S(:,:) double precision, intent(in) :: cutoff double precision :: local_cutoff integer :: info, i, j, k, mm if (n < 2) then return endif allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n)) call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) D(:) = dsqrt(D(:)) local_cutoff = dsqrt(cutoff)*D(1) ! such that D(i)/D(1) > dsqrt(cutoff) is kept mm=n do i=1,n if ( D(i) >= local_cutoff) then D(i) = 1.d0/D(i) else mm = mm-1 D(i) = 0.d0 endif do j=1,n S(j,i) = (0.d0,0.d0) enddo enddo if (mm < n) then print *, 'Removed Linear dependencies below ', local_cutoff endif !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S,U,D,Vt,n,C,m,local_cutoff) & !$OMP PRIVATE(i,j,k) do k=1,n if (D(k) /= 0.d0) then !$OMP DO SCHEDULE(STATIC) do j=1,n do i=1,n S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j) enddo enddo !$OMP END DO NOWAIT endif enddo !$OMP BARRIER !$OMP DO do j=1,n do i=1,m U(i,j) = C(i,j) enddo enddo !$OMP END DO !$OMP END PARALLEL call zgemm('N','N',m,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) deallocate(U,Vt,S,D) end subroutine get_inverse_complex(A,LDA,m,C,LDC) implicit none BEGIN_DOC ! Returns the inverse of the square matrix A END_DOC integer, intent(in) :: m, LDA, LDC complex*16, intent(in) :: A(LDA,m) complex*16, intent(out) :: C(LDC,m) integer :: info,lwork integer, allocatable :: ipiv(:) complex*16,allocatable :: work(:) allocate (ipiv(m), work(m*m)) lwork = size(work) C(1:m,1:m) = A(1:m,1:m) call zgetrf(m,m,C,size(C,1),ipiv,info) if (info /= 0) then print *, info stop 'error in inverse (zgetrf)' endif call zgetri(m,C,size(C,1),ipiv,work,lwork,info) if (info /= 0) then print *, info stop 'error in inverse (zgetri)' endif deallocate(ipiv,work) end subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff) implicit none BEGIN_DOC ! Find C = A^-1 END_DOC integer, intent(in) :: m,n, LDA, LDC complex*16, intent(in) :: A(LDA,n) double precision, intent(in) :: cutoff complex*16, intent(out) :: C(LDC,m) double precision, allocatable :: D(:), rwork(:) complex*16, allocatable :: U(:,:), Vt(:,:), work(:), A_tmp(:,:) integer :: info, lwork integer :: i,j,k double precision :: d1 allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n),rwork(5*n)) do j=1,n do i=1,m A_tmp(i,j) = A(i,j) enddo enddo lwork = -1 call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info) if (info /= 0) then print *, info, ': SVD failed' stop endif lwork = int(real(work(1))) deallocate(work) allocate(work(lwork)) call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info) if (info /= 0) then print *, info, ':: SVD failed' stop 1 endif d1 = D(1) do i=1,n if (D(i) > cutoff*d1) then D(i) = 1.d0/D(i) else D(i) = 0.d0 endif enddo C = (0.d0,0.d0) do i=1,m do j=1,n do k=1,n C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) enddo enddo enddo deallocate(U,D,Vt,work,A_tmp,rwork) end subroutine lapack_diagd_diag_in_place_complex(eigvalues,eigvectors,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H(complex) ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax ! double precision, intent(out) :: eigvectors(nmax,n) complex*16, intent(inout) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) ! double precision, intent(in) :: H(nmax,n) complex*16,allocatable :: work(:) integer ,allocatable :: iwork(:) ! complex*16,allocatable :: A(:,:) double precision, allocatable :: rwork(:) integer :: lrwork, lwork, info, i,j,l,k, liwork ! print*,'Diagonalization by jacobi' ! print*,'n = ',n lwork = 2*n*n + 2*n lrwork = 2*n*n + 5*n+ 1 liwork = 5*n + 3 allocate (work(lwork),iwork(liwork),rwork(lrwork)) lwork = -1 liwork = -1 lrwork = -1 ! get optimal work size call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = int( real(work(1))) liwork = iwork(1) lrwork = int(rwork(1)) deallocate (work,iwork,rwork) allocate (work(lwork),iwork(liwork),rwork(lrwork)) call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) deallocate(work,iwork,rwork) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'ZHEEVD Failed; calling ZHEEV' lwork = 2*n - 1 lrwork = 3*n - 2 allocate(work(lwork),rwork(lrwork)) lwork = -1 call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info) if (info < 0) then print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = int(work(1)) deallocate(work) allocate(work(lwork)) call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info) if (info /= 0 ) then write(*,*)'ZHEEV Failed' stop 1 endif deallocate(work,rwork) end if end subroutine lapack_diagd_diag_complex(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H(complex) ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax ! double precision, intent(out) :: eigvectors(nmax,n) complex*16, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) ! double precision, intent(in) :: H(nmax,n) complex*16, intent(in) :: H(nmax,n) double precision, allocatable :: eigenvalues(:) complex*16,allocatable :: work(:) integer ,allocatable :: iwork(:) complex*16,allocatable :: A(:,:) double precision, allocatable :: rwork(:) integer :: lrwork, lwork, info, i,j,l,k, liwork allocate(A(nmax,n),eigenvalues(n)) ! print*,'Diagonalization by jacobi' ! print*,'n = ',n A=H lwork = 2*n*n + 2*n lrwork = 2*n*n + 5*n+ 1 liwork = 5*n + 3 allocate (work(lwork),iwork(liwork),rwork(lrwork)) lwork = -1 liwork = -1 lrwork = -1 ! get optimal work size call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = int( real(work(1))) liwork = iwork(1) lrwork = int(rwork(1)) deallocate (work,iwork,rwork) allocate (work(lwork),iwork(liwork),rwork(lrwork)) call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) deallocate(work,iwork,rwork) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'ZHEEVD Failed; calling ZHEEV' lwork = 2*n - 1 lrwork = 3*n - 2 allocate(work(lwork),rwork(lrwork)) lwork = -1 call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info) if (info < 0) then print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = int(work(1)) deallocate(work) allocate(work(lwork)) call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info) if (info /= 0 ) then write(*,*)'ZHEEV Failed' stop 1 endif deallocate(work,rwork) end if eigvectors = (0.d0,0.d0) eigvalues = 0.d0 do j = 1, n eigvalues(j) = eigenvalues(j) do i = 1, n eigvectors(i,j) = A(i,j) enddo enddo deallocate(A,eigenvalues) end subroutine lapack_diagd_complex(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H(complex) ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax ! double precision, intent(out) :: eigvectors(nmax,n) complex*16, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) ! double precision, intent(in) :: H(nmax,n) complex*16, intent(in) :: H(nmax,n) double precision, allocatable :: eigenvalues(:) complex*16,allocatable :: work(:) integer ,allocatable :: iwork(:) complex*16,allocatable :: A(:,:) double precision, allocatable :: rwork(:) integer :: lrwork, lwork, info, i,j,l,k, liwork allocate(A(nmax,n),eigenvalues(n)) ! print*,'Diagonalization by jacobi' ! print*,'n = ',n A=H lwork = 2*n*n + 2*n lrwork = 2*n*n + 5*n+ 1 liwork = 5*n + 3 allocate (work(lwork),iwork(liwork),rwork(lrwork)) lwork = -1 liwork = -1 lrwork = -1 call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = max(int( work( 1 ) ),lwork) liwork = iwork(1) lrwork = max(int(rwork(1),4),lrwork) deallocate (work,iwork,rwork) allocate (work(lwork),iwork(liwork),rwork(lrwork)) call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) deallocate(work,iwork,rwork) if (info < 0) then print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'ZHEEVD Failed' stop 1 end if eigvectors = (0.d0,0.d0) eigvalues = 0.d0 do j = 1, n eigvalues(j) = eigenvalues(j) do i = 1, n eigvectors(i,j) = A(i,j) enddo enddo deallocate(A,eigenvalues) end subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H (complex) ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax complex*16, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) complex*16, intent(in) :: H(nmax,n) double precision,allocatable :: eigenvalues(:) complex*16,allocatable :: work(:) complex*16,allocatable :: A(:,:) double precision,allocatable :: rwork(:) integer :: lwork, info, i,j,l,k,lrwork allocate(A(nmax,n),eigenvalues(n)) ! print*,'Diagonalization by jacobi' ! print*,'n = ',n A=H !lwork = 2*n*n + 6*n+ 1 lwork = 2*n - 1 lrwork = 3*n - 2 allocate (work(lwork),rwork(lrwork)) lwork = -1 call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, info ) if (info < 0) then print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' stop 2 endif lwork = int( work( 1 ) ) deallocate (work) allocate (work(lwork)) call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, info ) deallocate(work,rwork) if (info < 0) then print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'ZHEEV Failed : ', info do i=1,n do j=1,n print *, H(i,j) enddo enddo stop 1 end if eigvectors = (0.d0,0.d0) eigvalues = 0.d0 do j = 1, n eigvalues(j) = eigenvalues(j) do i = 1, n eigvectors(i,j) = A(i,j) enddo enddo deallocate(A,eigenvalues) end subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda) implicit none BEGIN_DOC ! performs u1 += u0 * matrix END_DOC integer, intent(in) :: sze,lda complex*16, intent(in) :: u0(sze) complex*16, intent(inout) :: u1(sze) complex*16, intent(in) :: matrix(lda,sze) integer :: i,j integer :: incx,incy incx = 1 incy = 1 !call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy) call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy) end subroutine ortho_canonical(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. ! ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! ! N : Overlap matrix is NxN (array is (LDA,N) ) ! ! C : Coefficients of the vectors to orthogonalize. On exit, ! orthogonal vectors ! ! LDC : leftmost dimension of C ! ! m : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC integer, intent(in) :: lda, ldc, n integer, intent(out) :: m double precision, intent(in) :: overlap(lda,n) double precision, intent(in) :: cutoff double precision, intent(inout) :: C(ldc,n) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) double precision, allocatable :: S(:,:) !DIR$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j double precision :: local_cutoff if (n < 2) then return endif allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) D(:) = dsqrt(D(:)) local_cutoff = dsqrt(cutoff)*D(1) ! such that D(i)/D(1) > dsqrt(cutoff) is kept m=n do i=1,n if ( D(i) >= local_cutoff ) then D(i) = 1.d0/D(i) else m = i-1 print *, 'Removed Linear dependencies below:', local_cutoff exit endif enddo do i=m+1,n D(i) = 0.d0 enddo do j=1,n do i=1,n S(i,j) = U(i,j)*D(j) enddo enddo do j=1,n do i=1,n U(i,j) = C(i,j) enddo enddo call dgemm('N','N',n,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) deallocate (U, Vt, D, S) end subroutine ortho_svd(A,LDA,m,n) implicit none BEGIN_DOC ! Orthogonalization via fast SVD ! ! 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) if (m < n) then call ortho_qr(A,LDA,m,n) endif double precision, allocatable :: U(:,:), D(:), Vt(:,:) allocate(U(m,n), D(n), Vt(n,n)) call SVD(A,LDA,U,size(U,1),D,Vt,size(Vt,1),m,n) integer :: i,j do j=1,n do i=1,m A(i,j) = U(i,j) enddo enddo deallocate(U,D, Vt) end subroutine 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 double precision, allocatable :: TAU(:), WORK(:) 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 ) LWORK=-1 call dorgqr(m, n, 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 dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO) deallocate(WORK,TAU) end subroutine ortho_qr_unblocked(A,LDA,m,n) implicit none BEGIN_DOC ! Orthogonalization using Q.R factorization ! ! A : matrix to orthogonalize ! ! LDA : leftmost dimension of A ! ! n : Number of rows of A ! ! m : Number of columns of A ! END_DOC integer, intent(in) :: m,n, LDA double precision, intent(inout) :: A(LDA,n) integer :: info double precision, allocatable :: TAU(:), WORK(:) allocate (TAU(n), WORK(n)) call dgeqr2( m, n, A, LDA, TAU, WORK, INFO ) call dorg2r(m, n, n, A, LDA, TAU, WORK, INFO) deallocate(WORK,TAU) end subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.S^-1/2 orthogonalization. ! ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! ! N : Overlap matrix is NxN (array is (LDA,N) ) ! ! C : Coefficients of the vectors to orthogonalize. On exit, ! orthogonal vectors ! ! LDC : leftmost dimension of C ! ! M : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC integer, intent(in) :: LDA, ldc, n, m double precision, intent(in) :: overlap(lda,n) double precision, intent(in) :: cutoff double precision, intent(inout) :: C(ldc,n) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) double precision, allocatable :: S(:,:) integer :: info, i, j, k, mm double precision :: local_cutoff if (n < 2) then return endif allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) D(:) = dsqrt(D(:)) local_cutoff = dsqrt(cutoff)*D(1) ! such that D(i)/D(1) > dsqrt(cutoff) is kept mm=n do i=1,n if ( D(i) >= local_cutoff) then D(i) = 1.d0/D(i) else mm = mm-1 D(i) = 0.d0 endif do j=1,n S(j,i) = 0.d0 enddo enddo if (mm < n) then print *, 'Removed Linear dependencies below ', local_cutoff endif !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S,U,D,Vt,n,C,m,cutoff) & !$OMP PRIVATE(i,j,k) do k=1,n if (D(k) /= 0.d0) then !$OMP DO do j=1,n do i=1,n S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j) enddo enddo !$OMP END DO NOWAIT endif enddo !$OMP BARRIER !$OMP DO do j=1,n do i=1,m U(i,j) = C(i,j) enddo enddo !$OMP END DO !$OMP END PARALLEL call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) deallocate(U,Vt,S,D) end subroutine get_inverse(A,LDA,m,C,LDC) implicit none BEGIN_DOC ! Returns the inverse of the square matrix A END_DOC integer, intent(in) :: m, LDA, LDC double precision, intent(in) :: A(LDA,m) double precision, intent(out) :: C(LDC,m) integer :: info,lwork integer, allocatable :: ipiv(:) double precision,allocatable :: work(:) allocate (ipiv(m), work(m*m)) lwork = size(work) C(1:m,1:m) = A(1:m,1:m) call dgetrf(m,m,C,size(C,1),ipiv,info) if (info /= 0) then print *, info stop 'error in inverse (dgetrf)' endif call dgetri(m,C,size(C,1),ipiv,work,lwork,info) if (info /= 0) then print *, info stop 'error in inverse (dgetri)' endif deallocate(ipiv,work) end subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff) BEGIN_DOC ! Find C = A^-1 END_DOC implicit none integer, intent(in) :: m, n, LDA, LDC double precision, intent(in) :: A(LDA,n) double precision, intent(in) :: cutoff double precision, intent(out) :: C(LDC,m) integer :: info, lwork integer :: i, j, k, n_svd double precision :: D1_inv double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:) allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n)) do j=1,n do i=1,m A_tmp(i,j) = A(i,j) enddo enddo lwork = -1 call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info) if (info /= 0) then print *, info, ': SVD failed' stop endif LWORK=max(5*min(m,n),int(WORK(1))) deallocate(work) allocate(work(lwork)) call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info) if (info /= 0) then print *, info, ':: SVD failed' stop 1 endif if(D(1) .lt. 1d-14) then print*, ' largest singular value is very small:', D(1) n_svd = 1 else n_svd = 0 D1_inv = 1.d0 / D(1) do i = 1, n if(D(i)*D1_inv > cutoff) then D(i) = 1.d0 / D(i) n_svd = n_svd + 1 else D(i) = 0.d0 endif enddo endif !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, j) & !$OMP SHARED (n, n_svd, D, Vt) !$OMP DO do j = 1, n do i = 1, n_svd Vt(i,j) = D(i) * Vt(i,j) enddo enddo !$OMP END DO !$OMP END PARALLEL call dgemm('T', 'T', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1)) ! C = 0.d0 ! do i=1,m ! do j=1,n ! do k=1,n_svd ! C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) ! enddo ! enddo ! enddo deallocate(U,D,Vt,work,A_tmp) end subroutine find_rotation(A,LDA,B,m,C,n) implicit none BEGIN_DOC ! Find A.C = B END_DOC integer, intent(in) :: m,n, LDA double precision, intent(in) :: A(LDA,n), B(LDA,n) double precision, intent(out) :: C(n,n) double precision, allocatable :: A_inv(:,:) allocate(A_inv(LDA,n)) call get_pseudo_inverse(A,LDA,m,n,A_inv,LDA,1.d-10) integer :: i,j,k call dgemm('N','N',n,n,m,1.d0,A_inv,n,B,LDA,0.d0,C,n) deallocate(A_inv) end subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n) implicit none BEGIN_DOC ! Apply the rotation found by find_rotation END_DOC integer, intent(in) :: m,n, LDA, LDB, LDR double precision, intent(in) :: R(LDR,n) double precision, intent(in) :: A(LDA,n) double precision, intent(out) :: B(LDB,n) call dgemm('N','N',m,n,n,1.d0,A,LDA,R,LDR,0.d0,B,LDB) end subroutine lapack_diagd(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax double precision, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) double precision, intent(in) :: H(nmax,n) double precision,allocatable :: eigenvalues(:) double precision,allocatable :: work(:) integer ,allocatable :: iwork(:) double precision,allocatable :: A(:,:) integer :: lwork, info, i,j,l,k, liwork allocate(A(nmax,n),eigenvalues(n)) ! print*,'Diagonalization by jacobi' ! print*,'n = ',n A=H lwork = 1 liwork = 1 allocate (work(lwork),iwork(liwork)) lwork = -1 liwork = -1 call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & iwork, liwork, info ) if (info < 0) then print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value' stop 2 endif ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 LWORK = max(int(work(1)), 2*n*n + 6*n+ 1) liwork = max(iwork(1), 5*n + 3) deallocate (work,iwork) allocate (work(lwork),iwork(liwork)) call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & iwork, liwork, info ) deallocate(work,iwork) if (info < 0) then print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'DSYEVD Failed' stop 1 end if eigvectors = 0.d0 eigvalues = 0.d0 do j = 1, n eigvalues(j) = eigenvalues(j) do i = 1, n eigvectors(i,j) = A(i,j) enddo enddo deallocate(A,eigenvalues) end subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H ! ! H is untouched between input and ouptut ! ! eigevalues(i) = ith lowest eigenvalue of the H matrix ! ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector ! END_DOC integer, intent(in) :: n,nmax double precision, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) double precision, intent(in) :: H(nmax,n) double precision,allocatable :: eigenvalues(:) double precision,allocatable :: work(:) double precision,allocatable :: A(:,:) integer :: lwork, info, i,j,l,k, liwork allocate(A(nmax,n),eigenvalues(n)) A=H lwork = 1 allocate (work(lwork)) lwork = -1 call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & info ) if (info < 0) then print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' stop 2 endif ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 LWORK = max(int(work(1)), 2*n*n + 6*n+ 1) deallocate (work) allocate (work(lwork)) call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & info ) deallocate(work) if (info < 0) then print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' stop 2 else if( info > 0 ) then write(*,*)'DSYEV Failed : ', info do i=1,n do j=1,n print *, H(i,j) enddo enddo stop 1 end if eigvectors = 0.d0 eigvalues = 0.d0 do j = 1, n eigvalues(j) = eigenvalues(j) do i = 1, n eigvectors(i,j) = A(i,j) enddo enddo deallocate(A,eigenvalues) end subroutine nullify_small_elements(m,n,A,LDA,thresh) implicit none integer, intent(in) :: m,n,LDA double precision, intent(inout) :: A(LDA,n) double precision, intent(in) :: thresh integer :: i,j double precision :: amax ! Find max value amax = 0.d0 do j=1,n do i=1,m amax = max(dabs(A(i,j)), amax) enddo enddo if (amax == 0.d0) return amax = 1.d0/amax ! Remove tiny elements do j=1,n do i=1,m if ( (dabs(A(i,j) * amax) < thresh).or.(dabs(A(i,j)) < 1.d-99) ) then A(i,j) = 0.d0 endif enddo enddo end subroutine restore_symmetry(m,n,A,LDA,thresh) implicit none BEGIN_DOC ! Tries to find the matrix elements that are the same, and sets them ! to the average value. ! If restore_symm is False, only nullify small elements END_DOC integer, intent(in) :: m,n,LDA double precision, intent(inout) :: A(LDA,n) double precision, intent(in) :: thresh double precision, allocatable :: copy(:), copy_sign(:) integer, allocatable :: key(:), ii(:), jj(:) integer :: sze, pi, pf, idx, i,j,k double precision :: average, val, thresh2 thresh2 = dsqrt(thresh) sze = m * n allocate(copy(sze),copy_sign(sze),key(sze),ii(sze),jj(sze)) ! Copy to 1D !$OMP PARALLEL if (m>100) & !$OMP SHARED(A,m,n,sze,copy_sign,copy,key,ii,jj) & !$OMP PRIVATE(i,j,k) & !$OMP DEFAULT(NONE) !$OMP DO do j = 1, n do i = 1, m k = i+(j-1)*m copy(k) = A(i,j) copy_sign(k) = sign(1.d0,copy(k)) copy(k) = -dabs(copy(k)) key(k) = k ii(k) = i jj(k) = j enddo enddo !$OMP END DO !$OMP END PARALLEL ! Sort call dsort(copy,key,sze) call iset_order(ii,key,sze) call iset_order(jj,key,sze) call dset_order(copy_sign,key,sze) !TODO ! Parallelization with OMP ! Symmetrize i = 1 do while (i < sze) pi = i pf = i ! Exit if the remaining elements are below thresh if (-copy(i) <= thresh) exit val = 1d0/copy(i) do while (dabs(val * copy(pf+1) - 1d0) < thresh2) pf = pf + 1 ! if pf == sze, copy(pf+1) will not be valid if (pf == sze) then exit endif enddo ! if pi and pf are different do the average from pi to pf if (pf - pi > 0) then average = 0d0 do j = pi, pf average = average + copy(j) enddo average = average / (pf-pi+1.d0) do j = pi, pf copy(j) = average enddo ! Update i i = pf endif ! Update i i = i + 1 enddo ! To nullify the remaining elements that are below the threshold if (i == sze) then if (-copy(i) <= thresh) then copy(i) = 0d0 endif else copy(i:) = 0.d0 endif !$OMP PARALLEL if (sze>10000) & !$OMP SHARED(m,sze,copy_sign,copy,key,A,ii,jj) & !$OMP PRIVATE(i,j,k,idx) & !$OMP DEFAULT(NONE) ! copy -> A !$OMP DO do k = 1, sze i = ii(k) j = jj(k) A(i,j) = sign(copy(k),copy_sign(k)) enddo !$OMP END DO !$OMP END PARALLEL deallocate(copy,copy_sign,key,ii,jj) end subroutine diag_nonsym_right(n, A, A_ldim, V, V_ldim, energy, E_ldim) implicit none integer, intent(in) :: n, A_ldim, V_ldim, E_ldim double precision, intent(in) :: A(A_ldim,n) double precision, intent(out) :: energy(E_ldim), V(V_ldim,n) character*1 :: JOBVL, JOBVR, BALANC, SENSE integer :: i, j integer :: ILO, IHI, lda, ldvl, ldvr, LWORK, INFO double precision :: ABNRM integer, allocatable :: iorder(:), IWORK(:) double precision, allocatable :: WORK(:), SCALE_array(:), RCONDE(:), RCONDV(:) double precision, allocatable :: Atmp(:,:), WR(:), WI(:), VL(:,:), VR(:,:), Vtmp(:) double precision, allocatable :: energy_loc(:), V_loc(:,:) allocate( Atmp(n,n), WR(n), WI(n), VL(1,1), VR(n,n) ) do i = 1, n do j = 1, n Atmp(j,i) = A(j,i) enddo enddo JOBVL = "N" ! computes the left eigenvectors JOBVR = "V" ! computes the right eigenvectors BALANC = "B" ! Diagonal scaling and Permutation for optimization SENSE = "V" ! Determines which reciprocal condition numbers are computed lda = n ldvr = n ldvl = 1 allocate( WORK(1), SCALE_array(n), RCONDE(n), RCONDV(n), IWORK(2*n-2) ) LWORK = -1 ! to ask for the optimal size of WORK call dgeevx( BALANC, JOBVL, JOBVR, SENSE & ! CHARACTERS , n, Atmp, lda & ! MATRIX TO DIAGONALIZE , WR, WI & ! REAL AND IMAGINARY PART OF EIGENVALUES , VL, ldvl, VR, ldvr & ! LEFT AND RIGHT EIGENVECTORS , ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV & ! OUTPUTS OF OPTIMIZATION , WORK, LWORK, IWORK, INFO ) if(INFO .ne. 0) then print*, 'dgeevx failed !!', INFO stop endif LWORK = max(int(work(1)), 1) ! this is the optimal size of WORK deallocate(WORK) allocate(WORK(LWORK)) call dgeevx( BALANC, JOBVL, JOBVR, SENSE & , n, Atmp, lda & , WR, WI & , VL, ldvl, VR, ldvr & , ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV & , WORK, LWORK, IWORK, INFO ) if(INFO .ne. 0) then print*, 'dgeevx failed !!', INFO stop endif deallocate( WORK, SCALE_array, RCONDE, RCONDV, IWORK ) deallocate( VL, Atmp ) allocate( energy_loc(n), V_loc(n,n) ) energy_loc = 0.d0 V_loc = 0.d0 i = 1 do while(i .le. n) ! print*, i, WR(i), WI(i) if( dabs(WI(i)) .gt. 1e-7 ) then print*, ' Found an imaginary component to eigenvalue' print*, ' Re(i) + Im(i)', i, WR(i), WI(i) energy_loc(i) = WR(i) do j = 1, n V_loc(j,i) = WR(i) * VR(j,i) - WI(i) * VR(j,i+1) enddo energy_loc(i+1) = WI(i) do j = 1, n V_loc(j,i+1) = WR(i) * VR(j,i+1) + WI(i) * VR(j,i) enddo i = i + 2 else energy_loc(i) = WR(i) do j = 1, n V_loc(j,i) = VR(j,i) enddo i = i + 1 endif enddo deallocate(WR, WI, VR) ! ordering ! do j = 1, n ! write(444, '(100(1X, F16.10))') (V_loc(j,i), i=1,5) ! enddo allocate( iorder(n) ) do i = 1, n iorder(i) = i enddo call dsort(energy_loc, iorder, n) do i = 1, n energy(i) = energy_loc(i) do j = 1, n V(j,i) = V_loc(j,iorder(i)) enddo enddo deallocate(iorder) ! do j = 1, n ! write(445, '(100(1X, F16.10))') (V_loc(j,i), i=1,5) ! enddo deallocate(V_loc, energy_loc) end subroutine diag_nonsym_right ! --- ! Taken from GammCor thanks to Michal Hapka :-) subroutine pivoted_cholesky( A, rank, tol, ndim, U) ! ! A = U**T * U ! ! matrix A is destroyed inside this subroutine ! Cholesky vectors are stored in U ! dimension of U: U(1:rank, 1:n) ! U is allocated inside this subroutine ! rank is the number of Cholesky vectors depending on tol ! integer :: ndim integer, intent(inout) :: rank double precision, intent(inout) :: A(ndim, ndim) double precision, intent(out) :: U(ndim, rank) double precision, intent(in) :: tol integer, dimension(:), allocatable :: piv double precision, dimension(:), allocatable :: work character, parameter :: uplo = 'L' integer :: LDA integer :: info integer :: k, l, rank0 rank0 = rank LDA = ndim allocate(piv(ndim)) allocate(work(2*ndim)) call dpstrf(uplo, ndim, A, LDA, piv, rank, tol, work, info) if (rank > rank0) then print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling' stop end if do k = 1, ndim A(k,k+1:ndim) = 0.00D+0 end do ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer ! to do the swapping in-place U(:,:) = 0.00D+0 do k = 1, ndim l = piv(k) U(l, 1:rank) = A(k,1:rank) end do end subroutine pivoted_cholesky subroutine exp_matrix(X,n,exp_X) implicit none double precision, intent(in) :: X(n,n) integer, intent(in):: n double precision, intent(out):: exp_X(n,n) BEGIN_DOC ! exponential of the matrix X: X has to be ANTI HERMITIAN !! ! ! taken from Hellgaker, jorgensen, Olsen book ! ! section evaluation of matrix exponential (Eqs. 3.1.29 to 3.1.31) END_DOC integer :: i double precision, allocatable :: r2_mat(:,:),eigvalues(:),eigvectors(:,:) double precision, allocatable :: matrix_tmp1(:,:),eigvalues_mat(:,:),matrix_tmp2(:,:) include 'constants.include.F' allocate(r2_mat(n,n),eigvalues(n),eigvectors(n,n)) allocate(eigvalues_mat(n,n),matrix_tmp1(n,n),matrix_tmp2(n,n)) ! r2_mat = X^2 in the 3.1.30 call get_A_squared(X,n,r2_mat) call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n) eigvalues=-eigvalues do i = 1,n ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 eigvalues(i) = dsqrt(eigvalues(i)) enddo if(.false.)then !!! For debugging and following the book intermediate ! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30 ! matrix_tmp1 = W t^2 print*,'eigvalues = ' do i = 1, n print*,i,eigvalues(i) write(*,'(100(F16.10,X))')eigvectors(:,i) enddo eigvalues_mat=0.d0 do i = 1,n eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i) enddo call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) print*,'r2_mat = ' do i = 1, n write(*,'(100(F16.10,X))')r2_mat(:,i) enddo print*,'r2_mat new = ' do i = 1, n write(*,'(100(F16.10,X))')matrix_tmp2(:,i) enddo endif ! building the exponential ! exp(X) = W cos(t) W^T + W t^-1 sin(t) W^T X as in Eq. 3.1.31 ! matrix_tmp1 = W cos(t) do i = 1,n eigvalues_mat(i,i) = dcos(eigvalues(i)) enddo ! matrix_tmp2 = W cos(t) call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) ! matrix_tmp2 = W cos(t) W^T call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) exp_X = matrix_tmp2 ! matrix_tmp2 = W t^-1 sin(t) W^T X do i = 1,n if(dabs(eigvalues(i)).gt.1.d-4)then eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i) else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6 eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 & + eigvalues(i)*eigvalues(i)*eigvalues(i)*eigvalues(i)*c_1_3*0.025d0 endif enddo ! matrix_tmp1 = W t^-1 sin(t) call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) ! matrix_tmp2 = W t^-1 sin(t) W^T call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) ! exp_X += matrix_tmp2 X call dgemm('N','N',n,n,n,1.d0,matrix_tmp2,size(matrix_tmp2,1), & X,size(X,1),1.d0,exp_X,size(exp_X,1)) end subroutine exp_matrix_taylor(X,n,exp_X,converged) implicit none BEGIN_DOC ! exponential of a general real matrix X using the Taylor expansion of exp(X) ! ! returns the logical converged which checks the convergence END_DOC double precision, intent(in) :: X(n,n) integer, intent(in):: n double precision, intent(out):: exp_X(n,n) logical :: converged double precision :: f integer :: i,iter double precision, allocatable :: Tpotmat(:,:),Tpotmat2(:,:) allocate(Tpotmat(n,n),Tpotmat2(n,n)) BEGIN_DOC ! exponential of X using Taylor expansion END_DOC Tpotmat(:,:)=0.D0 exp_X(:,:) =0.D0 do i=1,n Tpotmat(i,i)=1.D0 exp_X(i,i) =1.d0 end do iter=0 converged=.false. do while (.not.converged) iter+=1 f = 1.d0 / dble(iter) Tpotmat2(:,:) = Tpotmat(:,:) * f call dgemm('N','N', n,n,n,1.d0, & Tpotmat2, size(Tpotmat2,1), & X, size(X,1), 0.d0, & Tpotmat, size(Tpotmat,1)) exp_X(:,:) = exp_X(:,:) + Tpotmat(:,:) converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) end do if(.not.converged)then print*,'Warning !! exp_matrix_taylor did not converge !' endif end subroutine get_A_squared(A,n,A2) implicit none BEGIN_DOC ! A2 = A A where A is n x n matrix. Use the dgemm routine END_DOC double precision, intent(in) :: A(n,n) integer, intent(in) :: n double precision, intent(out):: A2(n,n) call dgemm('N','N',n,n,n,1.d0,A,size(A,1),A,size(A,1),0.d0,A2,size(A2,1)) end subroutine get_AB_prod(A,n,m,B,l,AB) implicit none BEGIN_DOC ! AB = A B where A is n x m, B is m x l. Use the dgemm routine END_DOC double precision, intent(in) :: A(n,m),B(m,l) integer, intent(in) :: n,m,l double precision, intent(out):: AB(n,l) if(size(A,2).ne.m.or.size(B,1).ne.m)then print*,'error in get_AB_prod ! ' print*,'matrices do not have the good dimension ' print*,'size(A,2) = ',size(A,2) print*,'size(B,1) = ',size(B,1) print*,'m = ',m stop endif call dgemm('N','N',n,l,m,1.d0,A,size(A,1),B,size(B,1),0.d0,AB,size(AB,1)) end