From 9d3d843bc7f50f74a0ca479061a9c1cc7180b293 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 16 Apr 2017 01:28:35 +0200 Subject: [PATCH] Debugged Davidson for large Ndet --- src/Davidson/u0Hu0.irp.f | 102 ++++++++++++++------- src/Determinants/spindeterminants.irp.f | 114 +++++++++++++++++------- 2 files changed, 152 insertions(+), 64 deletions(-) diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 50546a33..30ff0fa8 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -665,7 +665,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) integer, allocatable :: singles_b(:,:) integer, allocatable :: idx(:), idx0(:) logical, allocatable :: is_single_a(:) - integer :: maxab, n_singles_max, kcol_prev + integer :: maxab, n_singles_max, kcol_prev, nmax double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: v_t, s_t, u_t @@ -702,7 +702,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) !$OMP singles_alpha_size, sze_8, & !$OMP idx0, u_t, maxab, v_0, s_0) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & - !$OMP lcol, lrow, is_single_a,l_a, l_b, & + !$OMP lcol, lrow, is_single_a,l_a, l_b, nmax, & !$OMP buffer, singles, doubles, n_singles, n_doubles, & !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, v_t, s_t) @@ -756,11 +756,12 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol) l_a = psi_bilinear_matrix_columns_loc(lcol) - do while (l_a <= k_a) - l_a += 1 - enddo n_doubles=1 + + ! Loop over alpha singles + ! ----------------------- + do while ( l_a < psi_bilinear_matrix_columns_loc(lcol+1) ) lrow = psi_bilinear_matrix_rows(l_a) if (.not.is_single_a(lrow)) then @@ -778,8 +779,8 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) call get_s2(tmp_det,tmp_det2,N_int,sij) do l=1,N_st v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a) - s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a) v_t(l,l_a) = v_t(l,l_a) + hij * u_t(l,k_a) + s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a) s_t(l,l_a) = s_t(l,l_a) + sij * u_t(l,k_a) enddo enddo @@ -811,6 +812,10 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) enddo !$OMP END DO NOWAIT + + ! Single and double alpha excitations + ! =================================== + !$OMP DO SCHEDULE(static,1) do k_a=1,N_det @@ -828,23 +833,18 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) k_b = psi_bilinear_matrix_order_transp_reverse(k_a) - ! Get all single and double alpha excitations - ! =========================================== - spindet(1:N_int) = tmp_det(1:N_int,1) ! Loop inside the beta column to gather all the connected alphas - i=1 - l_a = k_a - lcol = psi_bilinear_matrix_columns(l_a) - do while (lcol == kcol) + l_a = k_a+1 + nmax = min(N_det_alpha_unique, N_det - l_a) + do i=1,nmax + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit lrow = psi_bilinear_matrix_rows(l_a) buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow) idx(i) = l_a - i = i +1 l_a = l_a+1 - if (l_a > N_det) exit - lcol = psi_bilinear_matrix_columns(l_a) enddo i = i-1 @@ -852,7 +852,6 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) buffer, idx, spindet, N_int, i, & singles, doubles, n_singles, n_doubles ) - ! Compute Hij for all alpha singles ! ---------------------------------- @@ -884,26 +883,38 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) enddo enddo + end do + !$OMP END DO NOWAIT + + + ! Single and double beta excitations + ! ================================== + + !$OMP DO SCHEDULE(static,1) + do k_b=1,N_det + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + krow = psi_bilinear_matrix_transp_rows(k_b) + kcol = psi_bilinear_matrix_transp_columns(k_b) - ! Get all single and double beta excitations - ! =========================================== + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol) spindet(1:N_int) = tmp_det(1:N_int,2) + k_a = psi_bilinear_matrix_transp_order(k_b) ! Loop inside the alpha row to gather all the connected betas - i=1 - l_b = k_b - - lrow = psi_bilinear_matrix_transp_rows(l_b) - do while (lrow == krow) + l_b = k_b+1 + nmax = min(N_det_beta_unique, N_det - l_b) + do i=1,nmax + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit lcol = psi_bilinear_matrix_transp_columns(l_b) buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol) idx(i) = l_b - i = i +1 l_b = l_b+1 - if (l_b > N_det) exit - lrow = psi_bilinear_matrix_transp_rows(l_b) enddo i = i-1 @@ -943,9 +954,25 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) enddo enddo - ! Diagonal contribution - ! --------------------- + end do + !$OMP END DO NOWAIT + + ! Diagonal contribution + ! ===================== + + !$OMP DO SCHEDULE(static,1) + do k_a=1,N_det + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol) + double precision, external :: diag_H_mat_elem, diag_S_mat_elem hij = diag_H_mat_elem(tmp_det,N_int) @@ -956,7 +983,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8) enddo end do - !$OMP END DO + !$OMP END DO NOWAIT !$OMP CRITICAL do l=1,N_st @@ -986,7 +1013,7 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze double precision, allocatable :: vt(:,:) integer, allocatable :: idx(:) - integer :: i,j, jj + integer :: i,j, jj, l double precision :: hij do i=1,n @@ -995,6 +1022,7 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze allocate(idx(0:n), vt(N_st,n)) Vt = 0.d0 + !$OMP PARALLEL DO DEFAULT(shared) PRIVATE(i,idx,jj,j,degree,exc,phase,hij,l) SCHEDULE(static,1) do i=2,n idx(0) = i call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) @@ -1012,11 +1040,19 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze ! if ((degree == 2).and.(exc(0,1,1)==1)) cycle ! if ((degree > 1)) cycle ! if (exc(0,1,2) /= 0) cycle +! if (exc(0,1,1) == 2) cycle +! if (exc(0,1,2) == 2) cycle +! if ((degree==1).and.(exc(0,1,2) == 1)) cycle call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) - vt (:,i) = vt (:,i) + hij*u_0(j,:) - vt (:,j) = vt (:,j) + hij*u_0(i,:) + do l=1,N_st + !$OMP ATOMIC + vt (l,i) = vt (l,i) + hij*u_0(j,l) + !$OMP ATOMIC + vt (l,j) = vt (l,j) + hij*u_0(i,l) + enddo enddo enddo + !$OMP END PARALLEL DO do i=1,n v_0(i,:) = v_0(i,:) + vt(:,i) enddo diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index afd34d13..73460d0b 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -389,8 +389,6 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows , (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order , (N_det) ] -&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ] -&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) ] use bitmasks implicit none BEGIN_DOC @@ -408,7 +406,7 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) PROVIDE psi_coef_sorted_bit - integer, allocatable :: to_sort(:) + integer*8, allocatable :: to_sort(:) integer, external :: get_index_in_psi_det_alpha_unique integer, external :: get_index_in_psi_det_beta_unique allocate(to_sort(N_det)) @@ -421,16 +419,47 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) enddo psi_bilinear_matrix_rows(k) = i psi_bilinear_matrix_columns(k) = j - to_sort(k) = N_det_alpha_unique * (j-1) + i + to_sort(k) = int(N_det_alpha_unique,8) * int(j-1,8) + int(i,8) psi_bilinear_matrix_order(k) = k enddo - call isort(to_sort, psi_bilinear_matrix_order, N_det) + call i8sort(to_sort, psi_bilinear_matrix_order, N_det) call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det) call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det) do l=1,N_states call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det) enddo - psi_bilinear_matrix_columns_loc(1) = 1 + deallocate(to_sort) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ] + use bitmasks + implicit none + BEGIN_DOC +! Order which allors to go from psi_bilinear_matrix to psi_det + END_DOC + integer :: k + do k=1,N_det + psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_order(k)) = k + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) ] + use bitmasks + implicit none + BEGIN_DOC +! Sparse coefficient matrix if the wave function is expressed in a bilinear form : +! D_a^t C D_b +! +! Rows are alpha determinants and columns are beta. +! +! Order refers to psi_det + END_DOC + integer :: i,j,k, l + + l = psi_bilinear_matrix_columns(1) + psi_bilinear_matrix_columns_loc(l) = 1 do k=2,N_det if (psi_bilinear_matrix_columns(k) == psi_bilinear_matrix_columns(k-1)) then cycle @@ -440,35 +469,27 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) endif enddo psi_bilinear_matrix_columns_loc(N_det_beta_unique+1) = N_det+1 - do k=1,N_det - psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_order(k)) = k - enddo - deallocate(to_sort) END_PROVIDER - BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows , (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order , (N_det) ] -&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_transp_reverse , (N_det) ] use bitmasks implicit none BEGIN_DOC -! Sparse coefficient matrix if the wave function is expressed in a bilinear form : -! D_a^t C D_b +! Transpose of psi_bilinear_matrix +! D_b^t C^t D_a ! ! Rows are Alpha determinants and columns are beta, but the matrix is stored in row major ! format -! -! Order refers to psi_bilinear_matrix END_DOC integer :: i,j,k,l PROVIDE psi_coef_sorted_bit - integer, allocatable :: to_sort(:) + integer*8, allocatable :: to_sort(:) allocate(to_sort(N_det)) do l=1,N_states do k=1,N_det @@ -480,19 +501,50 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_ psi_bilinear_matrix_transp_rows (k) = psi_bilinear_matrix_rows (k) i = psi_bilinear_matrix_transp_columns(k) j = psi_bilinear_matrix_transp_rows (k) - to_sort(k) = N_det_beta_unique * (j-1) + i + to_sort(k) = int(N_det_beta_unique,8) * int(j-1,8) + int(i,8) psi_bilinear_matrix_transp_order(k) = k enddo - call isort(to_sort, psi_bilinear_matrix_transp_order, N_det) + call i8sort(to_sort, psi_bilinear_matrix_transp_order, N_det) call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det) do l=1,N_states call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det) enddo + deallocate(to_sort) +END_PROVIDER + +BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows_loc, (N_det_alpha_unique+1) ] + use bitmasks + implicit none + BEGIN_DOC +! Location of the columns in the psi_bilinear_matrix + END_DOC + integer :: i,j,k, l + + l = psi_bilinear_matrix_transp_rows(1) + psi_bilinear_matrix_transp_rows_loc(l) = 1 + do k=2,N_det + if (psi_bilinear_matrix_transp_rows(k) == psi_bilinear_matrix_transp_rows(k-1)) then + cycle + else + l = psi_bilinear_matrix_transp_rows(k) + psi_bilinear_matrix_transp_rows_loc(l) = k + endif + enddo + psi_bilinear_matrix_transp_rows_loc(N_det_beta_unique+1) = N_det+1 +END_PROVIDER + +BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_transp_reverse , (N_det) ] + use bitmasks + implicit none + BEGIN_DOC +! Order which allows to go from psi_bilinear_matrix_order_transp to psi_bilinear_matrix + END_DOC + integer :: k + do k=1,N_det psi_bilinear_matrix_order_transp_reverse(psi_bilinear_matrix_transp_order(k)) = k enddo - deallocate(to_sort) END_PROVIDER @@ -641,17 +693,17 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree -! select case (Nint) -! case (1) -! call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles) -! return -! case (2) -! call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) -! return -! case (3) -! call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) -! return -! end select + select case (Nint) + case (1) + call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles) + return + case (2) + call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) + return + case (3) + call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) + return + end select size_buffer_align = align_double(size_buffer)