From 07bfa1cf70733d7906084c3e0bc54e6c936b7b63 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Feb 2021 15:35:10 +0100 Subject: [PATCH] Removed orthonormalization in Davidson --- src/davidson/diagonalization_h_dressed.irp.f | 38 +++++++++-------- .../diagonalization_hcsf_dressed.irp.f | 42 +++++++++---------- .../diagonalization_hs2_dressed.irp.f | 38 +++++++++-------- 3 files changed, 62 insertions(+), 56 deletions(-) diff --git a/src/davidson/diagonalization_h_dressed.irp.f b/src/davidson/diagonalization_h_dressed.irp.f index ebecf623..870c58c7 100644 --- a/src/davidson/diagonalization_h_dressed.irp.f +++ b/src/davidson/diagonalization_h_dressed.irp.f @@ -298,14 +298,6 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - if (disk_based) then - call ortho_qr_unblocked(U,size(U,1),sze,shift2) - call ortho_qr_unblocked(U,size(U,1),sze,shift2) - else - call ortho_qr(U,size(U,1),sze,shift2) - call ortho_qr(U,size(U,1),sze,shift2) - endif - if ((sze > 100000).and.distributed_davidson) then call H_u_0_nstates_zmq (W(1,shift+1),U(1,shift+1),N_st_diag,sze) else @@ -359,11 +351,30 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia call dgemm('T','N', shift2, shift2, sze, & 1.d0, U, size(U,1), W, size(W,1), & 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) ! Diagonalize h ! --------------- - call lapack_diag(lambda,y,h,size(h,1),shift2) + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif ! Compute Energy for each eigenvector ! ----------------------------------- @@ -459,7 +470,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia endif do k=1,N_st - if (residual_norm(k) > 1.e8) then + if (residual_norm(k) > 1.d8) then print *, 'Davidson failed' stop -1 endif @@ -497,13 +508,6 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia enddo enddo - if (disk_based) then - call ortho_qr_unblocked(U,size(U,1),sze,N_st_diag) - call ortho_qr_unblocked(U,size(U,1),sze,N_st_diag) - else - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) - endif ! Adjust the phase do j=1,N_st_diag diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 54db55a8..d60b848b 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -15,7 +15,6 @@ subroutine davidson_diag_h_csf(dets_in,u_in,dim_in,energies,sze,sze_csf,N_st,N_s ! ! N_st : Number of eigenstates ! - ! Initial guess vectors are not necessarily orthonormal END_DOC integer, intent(in) :: dim_in, sze, sze_csf, N_st, N_st_diag, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) @@ -80,7 +79,6 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze ! - ! Initial guess vectors are not necessarily orthonormal END_DOC integer, intent(in) :: dim_in, sze, sze_csf, N_st, N_st_diag_in, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) @@ -302,14 +300,6 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - if (disk_based) then - call ortho_qr_unblocked(U_csf,size(U_csf,1),sze_csf,shift2) - call ortho_qr_unblocked(U_csf,size(U_csf,1),sze_csf,shift2) - else - call ortho_qr(U_csf,size(U_csf,1),sze_csf,shift2) - call ortho_qr(U_csf,size(U_csf,1),sze_csf,shift2) - endif - call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U) if ((sze > 100000).and.distributed_davidson) then call H_u_0_nstates_zmq (W,U,N_st_diag,sze) @@ -366,11 +356,30 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N call dgemm('T','N', shift2, shift2, sze_csf, & 1.d0, U_csf, size(U_csf,1), W_csf, size(W_csf,1), & 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze_csf, & + 1.d0, U_csf, size(U_csf,1), U_csf, size(U_csf,1), & + 0.d0, s_tmp, size(s_tmp,1)) ! Diagonalize h ! --------------- - call lapack_diag(lambda,y,h,size(h,1),shift2) + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif ! Compute Energy for each eigenvector ! ----------------------------------- @@ -438,9 +447,6 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) do k=1,N_st_diag do i=1,sze -! U_csf(i,shift2+k) = & -! (lambda(k) * U_csf(i,shift2+k) - W_csf(i,shift2+k) ) & -! /max(H_jj_csf(i) - lambda (k),1.d-2) U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) & /max(H_jj(i) - lambda (k),1.d-2) enddo @@ -509,14 +515,6 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N enddo enddo - if (disk_based) then - call ortho_qr_unblocked(U_csf,size(U_csf,1),sze_csf,N_st_diag) - call ortho_qr_unblocked(U_csf,size(U_csf,1),sze_csf,N_st_diag) - else - call ortho_qr(U_csf,size(U_csf,1),sze_csf,N_st_diag) - call ortho_qr(U_csf,size(U_csf,1),sze_csf,N_st_diag) - endif - call convertWFfromCSFtoDET(N_st_diag,U_csf,U) ! Adjust the phase diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 8c47c3e0..7461037e 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -353,14 +353,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - if (disk_based) then - call ortho_qr_unblocked(U,size(U,1),sze,shift2) - call ortho_qr_unblocked(U,size(U,1),sze,shift2) - else - call ortho_qr(U,size(U,1),sze,shift2) - call ortho_qr(U,size(U,1),sze,shift2) - endif - if ((sze > 100000).and.distributed_davidson) then call H_S2_u_0_nstates_zmq (W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze) else @@ -443,6 +435,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ call dgemm('T','N', shift2, shift2, sze, & 1.d0, U, size(U,1), W, size(W,1), & 0.d0, h, size(h_p,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) ! Penalty method ! -------------- @@ -467,7 +462,23 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! Diagonalize h_p ! --------------- - call lapack_diag(lambda,y,h_p,size(h_p,1),shift2) + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif ! Compute Energy for each eigenvector ! ----------------------------------- @@ -616,7 +627,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ endif do k=1,N_st - if (residual_norm(k) > 1.e8) then + if (residual_norm(k) > 1.d8) then print *, 'Davidson failed' stop -1 endif @@ -662,13 +673,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ enddo enddo - if (disk_based) then - call ortho_qr_unblocked(U,size(U,1),sze,N_st_diag) - call ortho_qr_unblocked(U,size(U,1),sze,N_st_diag) - else - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) - endif ! Adjust the phase do j=1,N_st_diag