mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-31 08:35:53 +01:00
Fix segfault in cholesky due to array syntax
This commit is contained in:
parent
b102cc816a
commit
a4d7648bb0
@ -23,38 +23,38 @@ END_PROVIDER
|
||||
!
|
||||
! Last dimension of cholesky_ao is cholesky_ao_num
|
||||
END_DOC
|
||||
|
||||
|
||||
integer :: rank, ndim
|
||||
double precision :: tau
|
||||
double precision, pointer :: L(:,:), L_old(:,:)
|
||||
|
||||
|
||||
|
||||
|
||||
double precision :: s
|
||||
double precision, parameter :: dscale = 1.d0
|
||||
|
||||
|
||||
double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
|
||||
integer, allocatable :: Lset(:), Dset(:), addr(:,:)
|
||||
logical, allocatable :: computed(:)
|
||||
|
||||
|
||||
integer :: i,j,k,m,p,q, qj, dj, p2, q2
|
||||
integer :: N, np, nq
|
||||
|
||||
|
||||
double precision :: Dmax, Dmin, Qmax, f
|
||||
double precision, external :: get_ao_two_e_integral
|
||||
logical, external :: ao_two_e_integral_zero
|
||||
|
||||
|
||||
double precision, external :: ao_two_e_integral
|
||||
integer :: block_size, iblock, ierr
|
||||
|
||||
|
||||
double precision :: mem
|
||||
double precision, external :: memory_of_double, memory_of_int
|
||||
|
||||
|
||||
integer, external :: getUnitAndOpen
|
||||
integer :: iunit
|
||||
|
||||
|
||||
ndim = ao_num*ao_num
|
||||
deallocate(cholesky_ao)
|
||||
|
||||
|
||||
if (read_ao_cholesky) then
|
||||
print *, 'Reading Cholesky vectors from disk...'
|
||||
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
|
||||
@ -63,11 +63,11 @@ END_PROVIDER
|
||||
read(iunit) cholesky_ao
|
||||
close(iunit)
|
||||
cholesky_ao_num = rank
|
||||
|
||||
|
||||
else
|
||||
|
||||
|
||||
PROVIDE nucl_coord
|
||||
|
||||
|
||||
if (do_direct_integrals) then
|
||||
if (ao_two_e_integral(1,1,1,1) < huge(1.d0)) then
|
||||
! Trigger providers inside ao_two_e_integral
|
||||
@ -76,16 +76,16 @@ END_PROVIDER
|
||||
else
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
endif
|
||||
|
||||
|
||||
tau = ao_cholesky_threshold
|
||||
|
||||
|
||||
mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim)
|
||||
call check_mem(mem, irp_here)
|
||||
|
||||
|
||||
call print_memory_usage()
|
||||
|
||||
allocate(L(ndim,1))
|
||||
|
||||
|
||||
print *, ''
|
||||
print *, 'Cholesky decomposition of AO integrals'
|
||||
print *, '======================================'
|
||||
@ -93,13 +93,13 @@ END_PROVIDER
|
||||
print *, '============ ============='
|
||||
print *, ' Rank Threshold'
|
||||
print *, '============ ============='
|
||||
|
||||
|
||||
|
||||
|
||||
rank = 0
|
||||
|
||||
|
||||
allocate( D(ndim), Lset(ndim), Dset(ndim) )
|
||||
allocate( addr(3,ndim) )
|
||||
|
||||
|
||||
! 1.
|
||||
k=0
|
||||
do j=1,ao_num
|
||||
@ -110,7 +110,7 @@ END_PROVIDER
|
||||
addr(3,k) = (i-1)*ao_num + j
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
if (do_direct_integrals) then
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided)
|
||||
do i=1,ndim
|
||||
@ -127,9 +127,9 @@ END_PROVIDER
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
endif
|
||||
|
||||
|
||||
Dmax = maxval(D)
|
||||
|
||||
|
||||
! 2.
|
||||
np=0
|
||||
do p=1,ndim
|
||||
@ -138,26 +138,26 @@ END_PROVIDER
|
||||
Lset(np) = p
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
! 3.
|
||||
N = 0
|
||||
|
||||
|
||||
! 4.
|
||||
i = 0
|
||||
|
||||
|
||||
! 5.
|
||||
do while ( (Dmax > tau).and.(rank < ndim) )
|
||||
! a.
|
||||
i = i+1
|
||||
|
||||
|
||||
s = 0.01d0
|
||||
|
||||
|
||||
! Inrease s until the arrays fit in memory
|
||||
do while (.True.)
|
||||
|
||||
|
||||
! b.
|
||||
Dmin = max(s*Dmax,tau)
|
||||
|
||||
|
||||
! c.
|
||||
nq=0
|
||||
do p=1,np
|
||||
@ -166,30 +166,30 @@ END_PROVIDER
|
||||
Dset(nq) = Lset(p)
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
call total_memory(mem)
|
||||
mem = mem &
|
||||
+ np*memory_of_double(nq) &! Delta(np,nq)
|
||||
+ (rank+nq)* memory_of_double(ndim) &! L(ndim,rank+nq)
|
||||
+ (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
|
||||
|
||||
|
||||
if (mem > qp_max_mem) then
|
||||
s = s*2.d0
|
||||
else
|
||||
exit
|
||||
endif
|
||||
|
||||
|
||||
if ((s > 1.d0).or.(nq == 0)) then
|
||||
call print_memory_usage()
|
||||
print *, 'Not enough memory. Reduce cholesky threshold'
|
||||
stop -1
|
||||
endif
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
! d., e.
|
||||
block_size = max(N,24)
|
||||
|
||||
|
||||
L_old => L
|
||||
allocate(L(ndim,rank+nq), stat=ierr)
|
||||
if (ierr /= 0) then
|
||||
@ -197,48 +197,52 @@ END_PROVIDER
|
||||
print *, irp_here, ': allocation failed : (L(ndim,rank+nq))'
|
||||
stop -1
|
||||
endif
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(k)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(k,j)
|
||||
do k=1,rank
|
||||
L(:,k) = L_old(:,k)
|
||||
do j=1,ndim
|
||||
L(j,k) = L_old(j,k)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
|
||||
deallocate(L_old)
|
||||
|
||||
|
||||
allocate(Delta(np,nq), stat=ierr)
|
||||
if (ierr /= 0) then
|
||||
call print_memory_usage()
|
||||
print *, irp_here, ': allocation failed : (Delta(np,nq))'
|
||||
stop -1
|
||||
endif
|
||||
|
||||
|
||||
allocate(Ltmp_p(np,block_size), stat=ierr)
|
||||
if (ierr /= 0) then
|
||||
call print_memory_usage()
|
||||
print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))'
|
||||
stop -1
|
||||
endif
|
||||
|
||||
|
||||
allocate(Ltmp_q(nq,block_size), stat=ierr)
|
||||
if (ierr /= 0) then
|
||||
call print_memory_usage()
|
||||
print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))'
|
||||
stop -1
|
||||
endif
|
||||
|
||||
|
||||
|
||||
|
||||
allocate(computed(nq))
|
||||
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j)
|
||||
|
||||
|
||||
!$OMP DO
|
||||
do q=1,nq
|
||||
Delta(:,q) = 0.d0
|
||||
do j=1,np
|
||||
Delta(j,q) = 0.d0
|
||||
enddo
|
||||
computed(q) = .False.
|
||||
enddo
|
||||
!$OMP ENDDO NOWAIT
|
||||
|
||||
|
||||
!$OMP DO
|
||||
do k=1,N
|
||||
do p=1,np
|
||||
@ -249,36 +253,36 @@ END_PROVIDER
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
|
||||
|
||||
!$OMP BARRIER
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
if (N>0) then
|
||||
call dgemm('N','T', np, nq, N, -1.d0, &
|
||||
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
|
||||
endif
|
||||
|
||||
|
||||
! f.
|
||||
Qmax = D(Dset(1))
|
||||
do q=1,nq
|
||||
Qmax = max(Qmax, D(Dset(q)))
|
||||
enddo
|
||||
|
||||
|
||||
! g.
|
||||
|
||||
|
||||
iblock = 0
|
||||
do j=1,nq
|
||||
|
||||
|
||||
if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit
|
||||
! i.
|
||||
rank = N+j
|
||||
|
||||
|
||||
if (iblock == block_size) then
|
||||
call dgemm('N','T',np,nq,block_size,-1.d0, &
|
||||
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
|
||||
iblock = 0
|
||||
endif
|
||||
|
||||
|
||||
! ii.
|
||||
do dj=1,nq
|
||||
qj = Dset(dj)
|
||||
@ -286,9 +290,9 @@ END_PROVIDER
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
L(1:ndim, rank) = 0.d0
|
||||
|
||||
|
||||
if (.not.computed(dj)) then
|
||||
m = dj
|
||||
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided)
|
||||
@ -314,16 +318,16 @@ END_PROVIDER
|
||||
do p=1,np
|
||||
Ltmp_p(p,iblock) = Delta(p,dj)
|
||||
enddo
|
||||
|
||||
|
||||
! iv.
|
||||
if (iblock > 1) then
|
||||
call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,&
|
||||
Ltmp_p(1,iblock), 1)
|
||||
endif
|
||||
|
||||
|
||||
! iii.
|
||||
f = 1.d0/dsqrt(Qmax)
|
||||
|
||||
|
||||
!$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared)
|
||||
!$OMP DO
|
||||
do p=1,np
|
||||
@ -332,38 +336,38 @@ END_PROVIDER
|
||||
D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,iblock) * Ltmp_p(p,iblock)
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
|
||||
!$OMP DO
|
||||
do q=1,nq
|
||||
Ltmp_q(q,iblock) = L(Dset(q), rank)
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
Qmax = D(Dset(1))
|
||||
do q=1,nq
|
||||
Qmax = max(Qmax, D(Dset(q)))
|
||||
enddo
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
print '(I10, 4X, ES12.3)', rank, Qmax
|
||||
|
||||
|
||||
deallocate(computed)
|
||||
deallocate(Delta)
|
||||
deallocate(Ltmp_p)
|
||||
deallocate(Ltmp_q)
|
||||
|
||||
|
||||
! i.
|
||||
N = rank
|
||||
|
||||
|
||||
! j.
|
||||
Dmax = D(Lset(1))
|
||||
do p=1,np
|
||||
Dmax = max(Dmax, D(Lset(p)))
|
||||
enddo
|
||||
|
||||
|
||||
np=0
|
||||
do p=1,ndim
|
||||
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
|
||||
@ -371,9 +375,9 @@ END_PROVIDER
|
||||
Lset(np) = p
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
|
||||
if (ierr /= 0) then
|
||||
call print_memory_usage()
|
||||
@ -387,10 +391,10 @@ END_PROVIDER
|
||||
!$OMP END PARALLEL DO
|
||||
deallocate(L)
|
||||
cholesky_ao_num = rank
|
||||
|
||||
|
||||
print *, '============ ============='
|
||||
print *, ''
|
||||
|
||||
|
||||
if (write_ao_cholesky) then
|
||||
print *, 'Writing Cholesky vectors to disk...'
|
||||
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W')
|
||||
@ -401,9 +405,9 @@ END_PROVIDER
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
|
||||
print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
|
||||
print *, ''
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
@ -252,7 +252,8 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
|
||||
do a = 1, nV
|
||||
do j = 1, nO
|
||||
do i = 1, nO
|
||||
W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
|
||||
! W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
|
||||
W_oovo(i,j,a,u) = 2d0 * cc_space_v_oovo(i,j,a,u) - cc_space_v_oovo(j,i,a,u)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -514,10 +515,10 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
|
||||
cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, &
|
||||
0.d0, tmpB1, nV*block_size)
|
||||
|
||||
call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, 1.d0, &
|
||||
cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
|
||||
tmp_cc2, cholesky_mo_num, &
|
||||
1.d0, tmpB1, nV*block_size)
|
||||
call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, &
|
||||
1.d0, cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
|
||||
tmp_cc2, cholesky_mo_num, &
|
||||
1.d0, tmpB1, nV*block_size)
|
||||
|
||||
do beta = iblock, min(nV, iblock+block_size-1)
|
||||
do b = 1, nV
|
||||
@ -1107,37 +1108,27 @@ subroutine compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
|
||||
double precision, intent(in) :: t2(nO, nO, nV, nV)
|
||||
double precision, intent(out) :: g_occ(nO, nO)
|
||||
|
||||
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
|
||||
g_occ = H_oo
|
||||
|
||||
call dgemm('N','N',nO,nO,nV, &
|
||||
1d0, t1, size(t1,1), &
|
||||
cc_space_f_vo, size(cc_space_f_vo,1), &
|
||||
0d0, g_occ, size(g_occ,1))
|
||||
1d0, g_occ, size(g_occ,1))
|
||||
|
||||
!$omp parallel &
|
||||
!$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) &
|
||||
!$omp private(i,j,a,u) &
|
||||
!$omp default(none)
|
||||
!$omp do
|
||||
do i = 1, nO
|
||||
do u = 1, nO
|
||||
g_occ(u,i) = g_occ(u,i) + H_oo(u,i)
|
||||
enddo
|
||||
enddo
|
||||
!$omp end do
|
||||
double precision, allocatable :: X(:)
|
||||
allocate(X(cholesky_mo_num))
|
||||
call dgemv('N',cholesky_mo_num,nO*nV,2.d0, &
|
||||
cc_space_v_ov_chol, cholesky_mo_num, &
|
||||
t1, 1, 0.d0, X, 1)
|
||||
|
||||
!$omp do
|
||||
do i = 1, nO
|
||||
do j = 1, nO
|
||||
do a = 1, nV
|
||||
do u = 1, nO
|
||||
g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$omp end do
|
||||
!$omp end parallel
|
||||
call dgemv('T',cholesky_mo_num,nO*nO,1.d0, &
|
||||
cc_space_v_oo_chol, cholesky_mo_num, &
|
||||
X, 1, 1.d0, g_occ, 1)
|
||||
deallocate(X)
|
||||
|
||||
call dgemv('T',nO*nV,nO*nO,-1.d0, &
|
||||
cc_space_v_ovoo, nO*nV, &
|
||||
t1, 1, 1.d0, g_occ, 1)
|
||||
|
||||
end
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user