From e39a069fd552ce21b9e9f141782867f827bf44bd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 18 Sep 2019 18:40:08 +0200 Subject: [PATCH] Using dgemm in r2_cc --- devel/cc/amplitudes.irp.f | 18 ++ devel/cc/curly_W.irp.f | 15 +- devel/cc/integral_batches.irp.f | 16 ++ devel/cc/r2_cc.irp.f | 341 ++++++++++++++++++-------------- 4 files changed, 234 insertions(+), 156 deletions(-) diff --git a/devel/cc/amplitudes.irp.f b/devel/cc/amplitudes.irp.f index 34e6b99..c079bc6 100644 --- a/devel/cc/amplitudes.irp.f +++ b/devel/cc/amplitudes.irp.f @@ -92,3 +92,21 @@ BEGIN_PROVIDER [ double precision, t2_cc, (spin_occ_num,spin_occ_num,spin_vir_nu endif END_PROVIDER + +BEGIN_PROVIDER [ double precision, t2_cc2, (spin_occ_num,spin_vir_num,spin_occ_num,spin_vir_num) ] + implicit none + BEGIN_DOC + ! Amplitudes with swapped indices + END_DOC + integer :: i,j,a,b + do b=1,spin_vir_num + do a=1,spin_vir_num + do j=1,spin_occ_num + do i=1,spin_occ_num + t2_cc2(i,a,j,b) = t2_cc(i,j,a,b) + enddo + enddo + enddo + enddo +END_PROVIDER + diff --git a/devel/cc/curly_W.irp.f b/devel/cc/curly_W.irp.f index 1728cc6..89f559a 100644 --- a/devel/cc/curly_W.irp.f +++ b/devel/cc/curly_W.irp.f @@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, cWoooo, (spin_occ_num,spin_occ_num,spin_occ_n END_PROVIDER -BEGIN_PROVIDER [ double precision, cWovvo, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ] +BEGIN_PROVIDER [ double precision, cWovvo_prime, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ] implicit none BEGIN_DOC ! Curly W in occ-vir-vir-occ block @@ -50,14 +50,15 @@ BEGIN_PROVIDER [ double precision, cWovvo, (spin_occ_num,spin_vir_num,spin_vir_n integer :: a,b,e,f double precision :: x - cWovvo(:,:,:,:) = OVVO(:,:,:,:) + + cWovvo_prime(:,:,:,:) = OVVO_prime(:,:,:,:) do f=1,spin_vir_num do j=1,spin_occ_num - do e=1,spin_vir_num - do b=1,spin_vir_num + do b=1,spin_vir_num + do e=1,spin_vir_num do m=1,spin_occ_num - cWovvo(m,b,e,j) = cWovvo(m,b,e,j) + t1_cc(j,f)*OVVV(m,b,e,f) + cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) + t1_cc(j,f)*OVVV(m,b,e,f) end do end do end do @@ -69,7 +70,7 @@ BEGIN_PROVIDER [ double precision, cWovvo, (spin_occ_num,spin_vir_num,spin_vir_n do b=1,spin_vir_num do n=1,spin_occ_num do m=1,spin_occ_num - cWovvo(m,b,e,j) = cWovvo(m,b,e,j) - t1_cc(n,b)*OOVO(m,n,e,j) + cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) - t1_cc(n,b)*OOVO(m,n,e,j) end do end do end do @@ -83,7 +84,7 @@ BEGIN_PROVIDER [ double precision, cWovvo, (spin_occ_num,spin_vir_num,spin_vir_n x = 0.5d0*t2_cc(j,n,f,b) + t1_cc(j,f)*t1_cc(n,b) do e=1,spin_vir_num do m=1,spin_occ_num - cWovvo(m,b,e,j) = cWovvo(m,b,e,j) - x *OOVV(m,n,e,f) + cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) - x *OOVV(m,n,e,f) end do end do end do diff --git a/devel/cc/integral_batches.irp.f b/devel/cc/integral_batches.irp.f index 65bfeae..f3d61ed 100644 --- a/devel/cc/integral_batches.irp.f +++ b/devel/cc/integral_batches.irp.f @@ -93,6 +93,22 @@ BEGIN_PROVIDER [ double precision, OVVO, (spin_occ_num,spin_vir_num,spin_vir_num 1:spin_occ_num) END_PROVIDER +BEGIN_PROVIDER [ double precision, OVVO_prime, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + integer :: m,b,e,j + do j=1,spin_occ_num + do b=1,spin_vir_num + do e=1,spin_vir_num + do m=1,spin_occ_num + OVVO_prime(m,e,b,j) = OVVO(m,b,e,j) + enddo + enddo + enddo + enddo + FREE OVVO +END_PROVIDER BEGIN_PROVIDER [ double precision, VOVO, (spin_vir_num,spin_occ_num,spin_vir_num,spin_occ_num) ] implicit none diff --git a/devel/cc/r2_cc.irp.f b/devel/cc/r2_cc.irp.f index 8b50d4f..ff86091 100644 --- a/devel/cc/r2_cc.irp.f +++ b/devel/cc/r2_cc.irp.f @@ -1,170 +1,213 @@ BEGIN_PROVIDER [ double precision, r2_cc, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ] - implicit none - BEGIN_DOC + implicit none + BEGIN_DOC ! t2 residues in non-canonical CCSD - END_DOC + END_DOC integer :: i,j,m,n integer :: a,b,e,f double precision :: x - double precision, allocatable :: tmp_oo(:,:), tmp_oo2(:,:) - double precision, external :: u_dot_v + double precision, allocatable :: mat_A(:,:), mat_B(:,:), mat_C(:,:) - allocate ( tmp_oo(spin_occ_num,spin_occ_num) ) - r2_cc(:,:,:,:) = OOVV(:,:,:,:) - - do b=1,spin_vir_num - do a=1,spin_vir_num - if (a == b) cycle - do e=1,spin_vir_num - x = u_dot_v( t1_cc(1,b), c_spin_fock_matrix_mo_ov(1,e), spin_occ_num ) - x = c_spin_fock_matrix_mo_vv(b,e) - x*0.5d0 - if (x == 0.d0) cycle - tmp_oo(:,:) = t2_cc(:,:,a,e) * x - r2_cc(:,:,a,b) = r2_cc(:,:,a,b) + tmp_oo(:,:) - r2_cc(:,:,b,a) = r2_cc(:,:,b,a) - tmp_oo(:,:) - end do + allocate ( mat_A(spin_vir_num,spin_vir_num), mat_C(spin_vir_num,spin_vir_num), & + mat_B(spin_vir_num,spin_occ_num) ) - do i=1,spin_occ_num - do m=1,spin_occ_num - tmp_oo(m,i) = t2_cc(i,m,a,b) + do j=1,spin_occ_num + do i=1,spin_occ_num + mat_A(:,:) = t2_cc(i,j,:,:) + call dgemm('N', 'T', & + spin_vir_num, spin_vir_num, spin_vir_num, & + 1.d0, mat_A, size(mat_A,1), & + c_spin_fock_matrix_mo_vv, size(c_spin_fock_matrix_mo_vv,1),& + 0.d0, mat_C, size(mat_C,1) ) + + call dgemm('N', 'T', & + spin_vir_num, spin_occ_num, spin_vir_num, & + -0.5d0, mat_A, size(mat_A,1), & + c_spin_fock_matrix_mo_ov, size(c_spin_fock_matrix_mo_ov,1),& + 0.d0, mat_B, size(mat_B,1) ) + + call dgemm('N', 'N', & + spin_vir_num, spin_vir_num, spin_occ_num, & + 1.d0, mat_B, size(mat_B,1), & + t1_cc, size(t1_cc,1), & + 1.d0, mat_C, size(mat_C,1) ) + + do b=1,spin_vir_num + do a=1,spin_vir_num + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(a,b) - mat_C(b,a) enddo enddo - do j=1,spin_occ_num - do i=1,spin_occ_num - x = 0.d0 - do m=1,spin_occ_num - do e=1,spin_vir_num - x = x + c_spin_fock_matrix_mo_ov(m,e) * & - ( tmp_oo(m,j)*t1_cc(i,e) - & - tmp_oo(m,i)*t1_cc(j,e) ) - end do - end do - x = x + u_dot_v( tau_cc(1,1,a,b), cWoooo(1,1,i,j), spin_occ_num*spin_occ_num ) - x = x * 0.5d0 + & - u_dot_v( tmp_oo(1,j), c_spin_fock_matrix_mo_oo(1,i), spin_occ_num ) - & - u_dot_v( tmp_oo(1,i), c_spin_fock_matrix_mo_oo(1,j), spin_occ_num ) - do e=1,spin_vir_num - x = x - sum(t2_cc(j,:,a,e)*cWovvo(:,b,e,i)) - end do - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + x - end do - end do - - tmp_oo(:,:) = 0.d0 - do f=1,spin_vir_num - do e=1,spin_vir_num - x = 0.5d0*cWvvvv(a,b,e,f) - if (x == 0.d0) cycle - tmp_oo(:,:) = tmp_oo(:,:) + x*tau_cc(:,:,e,f) - end do - end do - r2_cc(:,:,a,b) = r2_cc(:,:,a,b) + tmp_oo(:,:) - - - do j=1,spin_occ_num - do e=1,spin_vir_num - do m=1,spin_occ_num - x = t1_cc(m,a)*OVVO(m,b,e,j) - do i=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t2_cc(i,m,a,e)*cWovvo(m,b,e,j) - t1_cc(i,e)*x - enddo - end do - end do - end do - + enddo + enddo + + + deallocate(mat_A,mat_B,mat_C) + + allocate ( mat_A(spin_occ_num,spin_occ_num), mat_C(spin_occ_num,spin_occ_num), & + mat_B(spin_occ_num,spin_vir_num) ) + + do b=1,spin_vir_num + do a=1,spin_vir_num + + call dgemm('N', 'N', & + spin_occ_num, spin_occ_num, spin_occ_num, & + 1.d0, t2_cc(1,1,a,b), size(t2_cc,1), & + c_spin_fock_matrix_mo_oo, size(c_spin_fock_matrix_mo_oo,1),& + 0.d0, mat_C, size(mat_C,1) ) + + + call dgemm('N', 'N', & + spin_occ_num, spin_vir_num, spin_occ_num, & + 0.5d0, t2_cc(1,1,a,b), size(t2_cc,1), & + c_spin_fock_matrix_mo_ov, size(c_spin_fock_matrix_mo_ov,1),& + 0.d0, mat_B, size(mat_B,1) ) + + call dgemm('N', 'T', & + spin_occ_num, spin_occ_num, spin_vir_num, & + 1.d0, mat_B, size(mat_B,1), & + t1_cc, size(t1_cc,1), & + 1.d0, mat_C, size(mat_C,1) ) - do e=1,spin_vir_num - do j=1,spin_occ_num - do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(j,e)*t1_cc(m,a)*OVVO(m,b,e,i) - end do - end do - end do - end do - - do e=1,spin_vir_num - do j=1,spin_occ_num - do m=1,spin_occ_num - do i=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t2_cc(i,m,b,e)*cWovvo(m,a,e,j) - end do - end do - end do - end do - - do j=1,spin_occ_num - do e=1,spin_vir_num - do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(i,e)*t1_cc(m,b)*OVVO(m,a,e,j) - end do - end do - end do - end do - - do e=1,spin_vir_num - do j=1,spin_occ_num - do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t2_cc(j,m,b,e)*cWovvo(m,a,e,i) - end do - end do - end do - end do - - do e=1,spin_vir_num - do j=1,spin_occ_num - do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t1_cc(j,e)*t1_cc(m,b)*OVVO(m,a,e,i) - end do - end do - end do - end do - - - do j=1,spin_occ_num - do e=1,spin_vir_num - do i=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(i,e)*VVVO(a,b,e,j) - end do - end do - end do - - - do e=1,spin_vir_num - do i=1,spin_occ_num - do j=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t1_cc(j,e)*VVVO(a,b,e,i) - end do - end do - end do - - do j=1,spin_occ_num do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t1_cc(m,a)*OVOO(m,b,i,j) - end do - end do - end do - - do j=1,spin_occ_num - do i=1,spin_occ_num - do m=1,spin_occ_num - r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(m,b)*OVOO(m,a,i,j) - end do - end do - end do + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(j,i) - mat_C(i,j) + enddo + enddo + end do end do + + deallocate(mat_A,mat_B,mat_C) + + call dgemm('T', 'N', & + spin_occ_num*spin_occ_num, spin_vir_num*spin_vir_num, spin_occ_num*spin_occ_num,& + 0.5d0, cWoooo, spin_occ_num*spin_occ_num, & + tau_cc, spin_occ_num*spin_occ_num, & + 1.d0, r2_cc, spin_occ_num*spin_occ_num) + + call dgemm('N', 'T', & + spin_occ_num*spin_occ_num, spin_vir_num*spin_vir_num, spin_vir_num*spin_vir_num,& + 0.5d0, tau_cc, spin_occ_num*spin_occ_num, & + cWvvvv, spin_vir_num*spin_vir_num, & + 1.d0, r2_cc, spin_occ_num*spin_occ_num) + + + double precision, allocatable :: tmp(:,:,:,:), tmp2(:,:,:,:) + allocate(tmp(spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num), & + tmp2(spin_occ_num,spin_vir_num,spin_occ_num,spin_vir_num) ) + - ! Final expression of the t2 residue + call dgemm('N', 'N', & + spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num,& + 1.0d0, t2_cc2, spin_occ_num*spin_vir_num, & + cWovvo_prime, spin_occ_num*spin_vir_num, & + 0.d0, tmp, spin_occ_num*spin_vir_num) + + + do e=1,spin_vir_num + do m=1,spin_occ_num + do a=1,spin_vir_num + do i=1,spin_occ_num + tmp2(i,a,m,e) = -t1_cc(i,e)*t1_cc(m,a) + enddo + enddo + enddo + enddo + + call dgemm('N', 'N', & + spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num,& + 1.0d0, tmp2, spin_occ_num*spin_vir_num, & + OVVO_prime, spin_occ_num*spin_vir_num, & + 1.d0, tmp, spin_occ_num*spin_vir_num) + + do b=1,spin_vir_num + do a=1,spin_vir_num + do j=1,spin_occ_num + do i=1,spin_occ_num + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + tmp(i,a,b,j) + r2_cc(j,i,a,b) = r2_cc(j,i,a,b) - tmp(i,a,b,j) + r2_cc(i,j,b,a) = r2_cc(i,j,b,a) - tmp(i,a,b,j) + r2_cc(j,i,b,a) = r2_cc(j,i,b,a) + tmp(i,a,b,j) + enddo + enddo + enddo + enddo + + + + deallocate(tmp,tmp2) + + allocate( mat_A(spin_vir_num,spin_occ_num), mat_C(spin_occ_num,spin_occ_num) ) +! +! do b=1,spin_vir_num +! do a=1,spin_vir_num +! +! do j=1,spin_occ_num +! do e=1,spin_vir_num +! mat_A(e,j) = VVVO(a,b,e,j) +! enddo +! enddo +! +! call dgemm('N', 'N', & +! spin_occ_num, spin_occ_num, spin_vir_num, & +! 1.0d0, t1_cc, spin_occ_num, & +! mat_A, size(mat_A,1), & +! 0.d0, mat_C, spin_occ_num) +! +! do j=1,spin_occ_num +! do i=1,spin_occ_num +! r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(i,j) - mat_C(j,i) +! enddo +! endd +! +! enddo +! enddo + + + do b=1,spin_vir_num + do a=1,spin_vir_num + + do j=1,spin_occ_num + do e=1,spin_vir_num + mat_A(e,j) = VVVO(a,b,e,j) + enddo + enddo + call dgemm('N', 'N', & + spin_occ_num, spin_occ_num, spin_vir_num, & + 1.0d0, t1_cc, spin_occ_num, & + mat_A, size(mat_A,1), & + 0.d0, mat_C, spin_occ_num) + do j=1,spin_occ_num + do i=1,spin_occ_num + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(i,j) - mat_C(j,i) + enddo + enddo + + end do + end do + + + do b=1,spin_vir_num + do a=1,spin_vir_num + do j=1,spin_occ_num + do i=1,spin_occ_num + do m=1,spin_occ_num + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t1_cc(m,a)*OVOO(m,b,i,j) + r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(m,b)*OVOO(m,a,i,j) + end do + end do + end do + end do + end do + + +! Final expression of the t2 residue + r2_cc(:,:,:,:) = delta_oovv(:,:,:,:)*t2_cc(:,:,:,:) - r2_cc(:,:,:,:) - + END_PROVIDER +