From 70ad9f31b390b5931392a52ae375fab7c1e9f241 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Feb 2025 11:58:46 +0100 Subject: [PATCH] Accelerated PT2-cholesky with dgemm --- external/ezfio | 2 +- src/cipsi/selection.irp.f | 28 +++++++++++++------ src/mo_two_e_ints/map_integrals.irp.f | 40 ++++++++++++++++++++++----- 3 files changed, 54 insertions(+), 16 deletions(-) diff --git a/external/ezfio b/external/ezfio index d02132ea..dba01c4f 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d02132ea79217c16fd24242e8f8b8a6c3ff68091 +Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index af7420c8..8c22ec85 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1478,17 +1478,21 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer, parameter :: bant=1 - double precision, allocatable :: hij_cache1(:), hij_cache2(:) - allocate (hij_cache1(mo_num),hij_cache2(mo_num)) +! double precision, allocatable :: hij_cache1(:), hij_cache2(:) + double precision, allocatable :: hij_cache1(:,:), hij_cache2(:,:) +! allocate (hij_cache1(mo_num),hij_cache2(mo_num)) PROVIDE mo_integrals_threshold if(sp == 3) then ! AB + + allocate(hij_cache1(mo_num,mo_num)) + h1 = p(1,1) h2 = p(1,2) + call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache1,mo_integrals_map) do p1=1, mo_num if(bannedOrb(p1, 1)) cycle - call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) do p2=1, mo_num if(bannedOrb(p2,2)) cycle if(banned(p1, p2, bant)) cycle ! rentable? @@ -1497,7 +1501,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) - hij = hij_cache1(p2) * phase + hij = hij_cache1(p2,p1) * phase end if if (dabs(hij) < mo_integrals_threshold) cycle !DIR$ LOOP COUNT AVG(4) @@ -1507,13 +1511,18 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do end do + deallocate(hij_cache1) + else ! AA BB + + allocate(hij_cache1(mo_num,mo_num),hij_cache2(mo_num,mo_num)) + p1 = p(1,sp) p2 = p(2,sp) + call get_mo_two_e_integrals_ij(p2,p1,mo_num,hij_cache1,mo_integrals_map) + call get_mo_two_e_integrals_ij(p1,p2,mo_num,hij_cache2,mo_integrals_map) do puti=1, mo_num if (bannedOrb(puti, sp)) cycle - call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map) - call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map) do putj=puti+1, mo_num if(bannedOrb(putj, sp)) cycle if(banned(puti, putj, bant)) cycle ! rentable? @@ -1522,7 +1531,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) if (dabs(hij) < mo_integrals_threshold) cycle else - hij = hij_cache1(putj) - hij_cache2(putj) + hij = hij_cache1(putj,puti) - hij_cache2(putj,puti) if (dabs(hij) < mo_integrals_threshold) cycle hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if @@ -1532,9 +1541,12 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) enddo end do end do + + deallocate(hij_cache1,hij_cache2) + end if - deallocate(hij_cache1,hij_cache2) +! deallocate(hij_cache1,hij_cache2) end diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 5a005d7b..8f3c72d0 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -262,8 +262,14 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) out_val, 1) else integer :: isplit - out_val = 0.d0 - do isplit=1,4 + call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), & + mo_integrals_cache_min-1, 1., & + cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, & + cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., & + out_val_sp, 1) + out_val(1:mo_integrals_cache_min-1) += out_val_sp(1:mo_integrals_cache_min-1) + out_val(1:mo_integrals_cache_min-1) = out_val_sp(1:mo_integrals_cache_min-1) + do isplit=2,4 call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), & mo_integrals_cache_min-1, 1., & cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, & @@ -442,16 +448,36 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map) double precision, intent(out) :: out_array(sze,sze) type(map_type), intent(inout) :: map integer :: j - real(integral_kind), allocatable :: tmp_val(:) if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max