PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 04:13:55 +01:00
Code Releases Activity

Debug cholesky

Browse Source
This commit is contained in:
Anthony Scemama 2024-06-02 19:03:05 +02:00
parent b743201efe
commit 38d386d36c
2 changed files with 166 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
config/gfortran_debug_mkl.cfg Normal file
View File

@ -0,0 +1,63 @@
# Common flags
##############
#
# -ffree-line-length-none : Needed for IRPF90 which produces long lines
# -lblas -llapack : Link with libblas and liblapack libraries provided by the system
# -I . : Include the curent directory (Mandatory)
#
# --ninja : Allow the utilisation of ninja. (Mandatory)
# --align=32 : Align all provided arrays on a 32-byte boundary
#
#
[COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC -std=legacy
LAPACK_LIB : -I${MKLROOT}/include -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_core -lpthread -lm -ldl -lmkl_gnu_thread -lgomp -fopenmp
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
# Global options
################
#
# 1 : Activate
# 0 : Deactivate
#
[OPTION]
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags
####################
#
# -Ofast : Disregard strict standards compliance. Enables all -O3 optimizations.
# It also enables optimizations that are not valid
# for all standard-compliant programs. It turns on
# -ffast-math and the Fortran-specific
# -fno-protect-parens and -fstack-arrays.
[OPT]
FCFLAGS : -Ofast
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
# Debugging flags
#################
#
# -fcheck=all : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
#
[DEBUG]
#FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
FCFLAGS : -g -mavx -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan
# OpenMP flags
#################
#
[OPENMP]
FC : -fopenmp
IRPF90_FLAGS : --openmp

143
src/ao_two_e_ints/cholesky.irp.f
View File

@ -29,21 +29,20 @@ END_PROVIDER
integer*8 :: ndim8 integer*8 :: ndim8
integer :: rank integer :: rank
double precision :: tau double precision :: tau, tau2
double precision, pointer :: L(:,:), L_old(:,:) double precision, pointer :: L(:,:), L_old(:,:)
double precision :: s double precision :: s
double precision :: dscale double precision :: dscale, dscale_tmp
double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:) double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
integer, allocatable :: addr1(:,:), addr2(:,:) integer, allocatable :: addr1(:), addr2(:)
integer**, allocatable :: Lset(:), Dset(:), addr3(:,:) integer*8, allocatable :: Lset(:), Dset(:), addr3(:)
logical, allocatable :: computed(:) logical, allocatable :: computed(:)
integer :: i,j,k,m,p,q, dj, p2, q2 integer :: i,j,k,m,p,q, dj, p2, q2
integer*8 :: i8, j8, p8, qj8 integer*8 :: i8, j8, p8, qj8
integer :: N, np, nq integer :: N, np, nq, npmax
double precision :: Dmax, Dmin, Qmax, f double precision :: Dmax, Dmin, Qmax, f
double precision, external :: get_ao_two_e_integral double precision, external :: get_ao_two_e_integral
@ -53,6 +52,7 @@ END_PROVIDER
integer :: block_size, iblock integer :: block_size, iblock
double precision :: mem double precision :: mem
double precision, external :: memory_of_double, memory_of_int
double precision, external :: memory_of_double8, memory_of_int8 double precision, external :: memory_of_double8, memory_of_int8
integer, external :: getUnitAndOpen integer, external :: getUnitAndOpen
@ -61,6 +61,9 @@ END_PROVIDER
ndim8 = ao_num*ao_num*1_8 ndim8 = ao_num*ao_num*1_8
deallocate(cholesky_ao) deallocate(cholesky_ao)
! TODO : Save L() to disk
if (read_ao_cholesky) then if (read_ao_cholesky) then
print *, 'Reading Cholesky vectors from disk...' print *, 'Reading Cholesky vectors from disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
@ -85,6 +88,7 @@ END_PROVIDER
endif endif
tau = ao_cholesky_threshold tau = ao_cholesky_threshold
tau2 = tau*tau
mem = 6.d0 * memory_of_double8(ndim8) + 6.d0 * memory_of_int8(ndim8) mem = 6.d0 * memory_of_double8(ndim8) + 6.d0 * memory_of_int8(ndim8)
call check_mem(mem, irp_here) call check_mem(mem, irp_here)
@ -92,25 +96,25 @@ END_PROVIDER
call print_memory_usage() call print_memory_usage()
allocate(L(ndim8,1)) allocate(L(ndim8,1))
print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8) !print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8)
print *, '' print *, ''
print *, 'Cholesky decomposition of AO integrals' print *, 'Cholesky decomposition of AO integrals'
print *, '======================================' print *, '======================================'
print *, '' print *, ''
print *, '============ ============ =============' print *, '============ ============='
print *, ' Rank Block size Threshold' print *, ' Rank Threshold'
print *, '============ ============ =============' print *, '============ ============='
rank = 0 rank = 0
allocate( D(ndim8), Lset(ndim8), Dset(ndim8) ) allocate( D(ndim8), Lset(ndim8), Dset(ndim8) )
allocate( addr1(ndim8), addr2(ndim8), addr3(ndim8), ) allocate( addr1(ndim8), addr2(ndim8), addr3(ndim8) )
print *, 'allocate : (D(ndim8))', memory_of_int8(ndim8) !print *, 'allocate : (D(ndim8))', memory_of_int8(ndim8)
print *, 'allocate : (Lset(ndim8))', memory_of_int8(ndim8) !print *, 'allocate : (Lset(ndim8))', memory_of_int8(ndim8)
print *, 'allocate : (Dset(ndim8))', memory_of_int8(ndim8) !print *, 'allocate : (Dset(ndim8))', memory_of_int8(ndim8)
print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) !print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
! 1. ! 1.
k=0 k=0
@ -124,14 +128,14 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
enddo enddo
if (do_direct_integrals) then if (do_direct_integrals) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(guided) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,16)
do i8=1,ndim8 do i8=1,ndim8
D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), & D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), &
addr1(i8), addr2(i8)) addr1(i8), addr2(i8))
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else else
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(guided) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,16)
do i8=1,ndim8 do i8=1,ndim8
D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), & D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), &
addr2(i8), addr2(i8), & addr2(i8), addr2(i8), &
@ -143,17 +147,22 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
Dmax = maxval(D) Dmax = maxval(D)
! 2. ! 2.
np = huge(1_4) npmax = huge(1_4)*1_8
np = npmax
dscale = 1.d0 dscale = 1.d0
do while (np == huge(1_4)) dscale_tmp = Dmax
do while (np == npmax)
np=0 np=0
do p8=1,ndim8 do p8=1,ndim8
if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then if ( dscale_tmp*D(p8) > tau2 ) then
np = np+1 np = np+1
Lset(np) = p8 Lset(np) = p8
if (np == huge(1_4)) then if (np == npmax) then
! Overflow detected ! Overflow detected
dscale = dscale*0.5d0 dscale = dscale*0.1d0
dscale_tmp = dscale*dscale*Dmax
!print *, 'Overflow detected '
!print *, 'dscale = ', dscale
exit exit
endif endif
endif endif
@ -167,7 +176,7 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
i = 0 i = 0
! 5. ! 5.
do while ( (Dmax > tau).and.(rank < min(ndim8,huge(1_4)) ) do while ( (Dmax > tau).and.(rank < min(ndim8,huge(1_4))) )
! a. ! a.
i = i+1 i = i+1
@ -193,13 +202,11 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
call total_memory(mem) call total_memory(mem)
mem = mem & mem = mem &
+ np*memory_of_double(nq) &! Delta(np,nq) + np*memory_of_double(nq) &! Delta(np,nq)
+ (rank+nq)* memory_of_double8(ndim8) + (rank+nq)* memory_of_double8(ndim8) &! L(ndim8,rank+nq)
&! L(ndim8,rank+nq)
+ (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size) + (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
if (mem > qp_max_mem) then if (mem > qp_max_mem) then
s = s*2.d0 s = s*2.d0
block_size = max(block_size / 2, 1)
else else
exit exit
endif endif
@ -219,7 +226,7 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8)
L_old => L L_old => L
allocate(L(ndim8,rank+nq), stat=ierr) allocate(L(ndim8,rank+nq), stat=ierr)
print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq)) !print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq))
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -227,7 +234,7 @@ print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq))
stop -1 stop -1
endif endif
!$OMP PARALLEL DO PRIVATE(k,j) !$OMP PARALLEL DO PRIVATE(k,j8)
do k=1,rank do k=1,rank
do j8=1,ndim8 do j8=1,ndim8
L(j8,k) = L_old(j8,k) L(j8,k) = L_old(j8,k)
@ -238,7 +245,7 @@ print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq))
deallocate(L_old) deallocate(L_old)
allocate(Delta(np,nq), stat=ierr) allocate(Delta(np,nq), stat=ierr)
print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8) !print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -247,7 +254,7 @@ print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8)
endif endif
allocate(Ltmp_p(np,block_size), stat=ierr) allocate(Ltmp_p(np,block_size), stat=ierr)
print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double8(np*block_size*1_8) !print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double8(np*block_size*1_8), np, block_size
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -256,7 +263,7 @@ print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double8(np*block_size*1_8
endif endif
allocate(Ltmp_q(nq,block_size), stat=ierr) allocate(Ltmp_q(nq,block_size), stat=ierr)
print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double8(nq*block_size*1_8) !print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double8(nq*block_size*1_8), nq, block_size
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -266,34 +273,47 @@ print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double8(nq*block_size*1_8
allocate(computed(nq)) allocate(computed(nq))
print *, 'allocate : computed(nq)', memory_of_int(nq) !print *, 'allocate : computed(nq)', memory_of_int(nq)
print *, 'p1' !print *, 'N, rank, block_size', N, rank, block_size
!print *, 'p1'
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q,j) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q,j)
!print *, 'computed'
!$OMP DO
do q=1,nq
computed(q) = .False.
enddo
!$OMP ENDDO NOWAIT
!print *, 'Delta'
!$OMP DO !$OMP DO
do q=1,nq do q=1,nq
do j=1,np do j=1,np
Delta(j,q) = 0.d0 Delta(j,q) = 0.d0
enddo enddo
computed(q) = .False. enddo
!$OMP ENDDO NOWAIT
!print *, 'Ltmp_p'
do k=1,N
!$OMP DO
do p=1,np
Ltmp_p(p,k) = L(Lset(p),k)
enddo enddo
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
!$OMP DO !$OMP DO
do k=1,N
do p=1,np
Ltmp_p(p,k) = L(Lset(p),k)
enddo
do q=1,nq do q=1,nq
Ltmp_q(q,k) = L(Dset(q),k) Ltmp_q(q,k) = L(Dset(q),k)
enddo enddo
enddo
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
enddo
!$OMP BARRIER !$OMP BARRIER
!$OMP END PARALLEL !$OMP END PARALLEL
!print *, 'p2', np, nq, N
if (N>0) then if (N>0) then
call dgemm('N','T', np, nq, N, -1.d0, & call dgemm('N','T', np, nq, N, -1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
@ -315,8 +335,10 @@ print *, 'p1'
rank = N+j rank = N+j
if (iblock == block_size) then if (iblock == block_size) then
!print *, 'dgemm'
call dgemm('N','T',np,nq,block_size,-1.d0, & call dgemm('N','T',np,nq,block_size,-1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
iblock = 0 iblock = 0
endif endif
@ -328,36 +350,47 @@ print *, 'p1'
endif endif
enddo enddo
L(1:ndim8, rank) = 0.d0 do i8=1,ndim8
L(i8, rank) = 0.d0
enddo
if (.not.computed(dj)) then if (.not.computed(dj)) then
m = dj m = dj
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided) if (do_direct_integrals) then
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,16)
do k=np,1,-1 do k=np,1,-1
if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)) ) ) then addr2(Lset(k)), addr2(Dset(m)) ) ) then
if (do_direct_integrals) then
Delta(k,m) = Delta(k,m) + & Delta(k,m) = Delta(k,m) + &
ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),&
addr1(Dset(m)), addr2(Dset(m))) addr1(Dset(m)), addr2(Dset(m)))
endif
enddo
!$OMP END PARALLEL DO
else else
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,16)
do k=np,1,-1
if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)) ) ) then
Delta(k,m) = Delta(k,m) + & Delta(k,m) = Delta(k,m) + &
get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map) addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map)
endif endif
endif
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif
computed(dj) = .True. computed(dj) = .True.
endif endif
iblock = iblock+1 iblock = iblock+1
!print *, iblock
do p=1,np do p=1,np
Ltmp_p(p,iblock) = Delta(p,dj) Ltmp_p(p,iblock) = Delta(p,dj)
enddo enddo
! iv. ! iv.
if (iblock > 1) then if (iblock > 1) then
!print *, 'dgemv', iblock
call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,& call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,&
Ltmp_p(1,iblock), 1) Ltmp_p(1,iblock), 1)
endif endif
@ -365,6 +398,7 @@ print *, 'p1'
! iii. ! iii.
f = 1.d0/dsqrt(Qmax) f = 1.d0/dsqrt(Qmax)
!print *, 'p4'
!$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared) !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared)
!$OMP DO !$OMP DO
do p=1,np do p=1,np
@ -379,7 +413,6 @@ print *, 'p1'
Ltmp_q(q,iblock) = L(Dset(q), rank) Ltmp_q(q,iblock) = L(Dset(q), rank)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
Qmax = D(Dset(1)) Qmax = D(Dset(1))
@ -389,12 +422,12 @@ print *, 'p1'
enddo enddo
print '(I10, 4X, I10, 4X, ES12.3)', rank, block_size, Qmax print '(I10, 4X, ES12.3)', rank, Qmax
deallocate(computed)
deallocate(Delta)
deallocate(Ltmp_p) deallocate(Ltmp_p)
deallocate(Ltmp_q) deallocate(Ltmp_q)
deallocate(computed)
deallocate(Delta)
! i. ! i.
N = rank N = rank
@ -405,17 +438,21 @@ print *, 'p1'
Dmax = max(Dmax, D(Lset(p))) Dmax = max(Dmax, D(Lset(p)))
enddo enddo
np = huge(1_4) np = npmax
dscale = 1.d0 dscale = 1.d0
do while (np == huge(1_4)) dscale_tmp = Dmax
do while (np == npmax)
np=0 np=0
do p8=1,ndim8 do p8=1,ndim8
if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then if ( dscale_tmp*D(p8) > tau2 ) then
np = np+1 np = np+1
Lset(np) = p8 Lset(np) = p8
if (np == huge(1_4)) then if (np == npmax) then
! Overflow detected ! Overflow detected
dscale = dscale*0.5d0 dscale = dscale*0.5d0
dscale_tmp = dscale*dscale*Dmax
!print *, 'Overflow detected '
!print *, 'dscale = ', dscale
exit exit
endif endif
endif endif
@ -426,7 +463,7 @@ print *, 'p1'
enddo enddo
allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', memory_of_double8(ao_num*ao_num*rank*1_8) !print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', memory_of_double8(ao_num*ao_num*rank*1_8)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -436,7 +473,7 @@ print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', memory_of_double8(ao_num*
!$OMP PARALLEL DO PRIVATE(k) !$OMP PARALLEL DO PRIVATE(k)
do k=1,rank do k=1,rank
do j=1,ao_num do j=1,ao_num
call dcopy(ao_num, L((j-1)*ao_num+1,k), 1, cholesky_ao(1,j,k), 1) cholesky_ao(1:ao_num,j,k) = L((j-1)*ao_num+1:j*ao_num,k)
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO