diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 0edd3f63..8b21798f 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -737,12 +737,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi = mat(istate, p1, p2) - do jstate=1,N_states - pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) - enddo - - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - pt2_data % pt2(istate) += e_pert(istate) + pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate) + pt2_data % variance(istate) = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi + pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate) !!!DEBUG ! delta_E = E0(istate) - Hii + E_shift @@ -1573,7 +1570,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * phase end if - mat(:, p1, p2) += coefs(:) * hij + mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij end do end do else ! AA BB @@ -1590,7 +1587,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, else hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do end if @@ -1649,18 +1646,18 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) end do do putj=hfix+1, mo_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) end do if(ma == 1) then - mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num) + mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num) else - mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num) + mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num) end if end if @@ -1674,22 +1671,22 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, putj = p1 if(.not. banned(putj,puti,bant)) then hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) - tmp_row(:,puti) += hij * coefs(:) + tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) end if putj = p2 if(.not. banned(putj,puti,bant)) then hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) - tmp_row2(:,puti) += hij * coefs(:) + tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) end if end do if(mi == 1) then - mat(:,:,p1) += tmp_row(:,:) - mat(:,:,p2) += tmp_row2(:,:) + mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:) + mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:) else - mat(:,p1,:) += tmp_row(:,:) - mat(:,p2,:) += tmp_row2(:,:) + mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:) + mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:) end if else if(p(0,ma) == 3) then @@ -1702,16 +1699,16 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(:,putj) += hij * coefs(:) + tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) end do do putj=hfix+1,mo_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(:,putj) += hij * coefs(:) + tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) end do - mat(:, :puti-1, puti) += tmp_row(:,:puti-1) - mat(:, puti, puti:) += tmp_row(:,puti:) + mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1) + mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:) end do else hfix = h(1,mi) @@ -1725,19 +1722,19 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, putj = p2 if(.not. banned(puti,putj,1)) then hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) - tmp_row(:,puti) += hij * coefs(:) + tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) end if putj = p1 if(.not. banned(puti,putj,1)) then hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) - tmp_row2(:,puti) += hij * coefs(:) + tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) end if end do - mat(:,:p2-1,p2) += tmp_row(:,:p2-1) - mat(:,p2,p2:) += tmp_row(:,p2:) - mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) - mat(:,p1,p1:) += tmp_row2(:,p1:) + mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1) + mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:) + mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1) + mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:) end if end if deallocate(lbanned) @@ -1760,7 +1757,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) call i_h_j(gen, det, N_int, hij) - mat(:, p1, p2) += coefs(:) * hij + mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij end do end do end @@ -1813,9 +1810,9 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) if(ma == 1) then - mat(:, putj, puti) += coefs(:) * hij + mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij else - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end if end do else @@ -1831,7 +1828,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p1 = p(turn2(i), 1) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do end if @@ -1851,7 +1848,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p1 = p(i1, ma) p2 = p(i2, ma) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do else if(tip == 3) then @@ -1865,7 +1862,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p2 = p(i, ma) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) - mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij + mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij end do else ! tip == 4 puti = p(1, sp) @@ -1876,7 +1873,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, h1 = h(1, mi) h2 = h(2, mi) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end if end if end if diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index e56fdec4..99bc3816 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -211,6 +211,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, double precision :: rss, mem, ratio double precision, allocatable :: utl(:,:) integer, parameter :: block_size=128 + logical :: u_is_sparse ! call resident_memory(rss) ! mem = dble(singles_beta_csc_size) / 1024.d0**3 @@ -222,6 +223,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, ! endif compute_singles=.True. + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 allocate(idx0(maxab)) @@ -249,7 +251,7 @@ compute_singles=.True. !$OMP singles_alpha_csc,singles_alpha_csc_idx, & !$OMP singles_beta_csc,singles_beta_csc_idx) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & - !$OMP lcol, lrow, l_a, l_b, utl, kk, & + !$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, & !$OMP buffer, doubles, n_doubles, umax, & !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, ratio, & @@ -266,6 +268,19 @@ compute_singles=.True. kcol_prev=-1 + ! Check if u has multiple zeros + kk=0 + do k=1,N_det + umax = 0.d0 + do l=1,N_st + umax = max(umax, dabs(u_t(l,k))) + enddo + if (umax < 1.d-20) then + kk = kk+1 + endif + enddo + u_is_sparse = N_det / kk < 10 + ASSERT (iend <= N_det) ASSERT (istart > 0) ASSERT (istep > 0) @@ -405,16 +420,25 @@ compute_singles=.True. do k = 1,n_singles_a,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (k+kk > n_singles_a) exit - l_a = singles_a(k+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (k+kk > n_singles_a) exit + l_a = singles_a(k+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (k+kk > n_singles_a) exit + l_a = singles_a(k+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -497,16 +521,25 @@ compute_singles=.True. do i=1,n_singles_a,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (i+kk > n_singles_a) exit - l_a = singles_a(i+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_singles_a) exit + l_a = singles_a(i+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_singles_a) exit + l_a = singles_a(i+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -534,16 +567,25 @@ compute_singles=.True. do i=1,n_doubles,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (i+kk > n_doubles) exit - l_a = doubles(i+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_a = doubles(i+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_a = doubles(i+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -611,19 +653,29 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(1000) do i=1,n_singles_b,block_size umax = 0.d0 - do kk=0,block_size-1 - if (i+kk > n_singles_b) exit - l_b = singles_b(i+kk) - ASSERT (l_b <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_singles_b) exit + l_b = singles_b(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) - l_a = psi_bilinear_matrix_transp_order(l_b) - ASSERT (l_a <= N_det) - - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_singles_b) exit + l_b = singles_b(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -649,18 +701,28 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(50000) do i=1,n_doubles,block_size umax = 0.d0 - do kk=0,block_size-1 - if (i+kk > n_doubles) exit - l_b = doubles(i+kk) - ASSERT (l_b <= N_det) - l_a = psi_bilinear_matrix_transp_order(l_b) - ASSERT (l_a <= N_det) - - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_b = doubles(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_b = doubles(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -688,10 +750,14 @@ compute_singles=.True. ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - umax = 0.d0 - do l=1,N_st - umax = max(umax, dabs(u_t(l,k_a))) - enddo + if (u_is_sparse) then + umax = 0.d0 + do l=1,N_st + umax = max(umax, dabs(u_t(l,k_a))) + enddo + else + umax = 1.d0 + endif if (umax < 1.d-20) cycle krow = psi_bilinear_matrix_rows(k_a)