From b773a361b204a1f8424ae26bc2ac1a1a9d424a9f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 31 May 2024 20:07:29 +0200 Subject: [PATCH 1/5] fixed uninitialized variable in cholesky' --- src/ao_two_e_ints/cholesky.irp.f | 38 ++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 5fbd166c..09d86679 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -22,6 +22,9 @@ END_PROVIDER ! = (ik|jl) = sum_a (ik|a).(a|jl) ! ! Last dimension of cholesky_ao is cholesky_ao_num + ! + ! https://mogp-emulator.readthedocs.io/en/latest/methods/proc/ProcPivotedCholesky.html + ! https://doi.org/10.1016/j.apnum.2011.10.001 : Page 4, Algorithm 1 END_DOC integer :: rank, ndim @@ -86,20 +89,25 @@ END_PROVIDER call print_memory_usage() allocate(L(ndim,1)) +!print *, 'allocate : (L(ndim,1))', memory_of_double(ndim) print *, '' print *, 'Cholesky decomposition of AO integrals' print *, '======================================' print *, '' - print *, '============ =============' - print *, ' Rank Threshold' - print *, '============ =============' + print *, '============ ============ =============' + print *, ' Rank Block size Threshold' + print *, '============ ============ =============' rank = 0 allocate( D(ndim), Lset(ndim), Dset(ndim) ) allocate( addr(3,ndim) ) +!print *, 'allocate : (D(ndim))', memory_of_int(ndim) +!print *, 'allocate : (Lset(ndim))', memory_of_int(ndim) +!print *, 'allocate : (Dset(ndim))', memory_of_int(ndim) +!print *, 'allocate : (3,addr(ndim))', memory_of_int(3*ndim) ! 1. k=0 @@ -151,9 +159,10 @@ END_PROVIDER ! a. i = i+1 - s = 0.01d0 ! Inrease s until the arrays fit in memory + s = 0.01d0 + block_size = max(N,24) do while (.True.) ! b. @@ -168,6 +177,7 @@ END_PROVIDER endif enddo + call total_memory(mem) mem = mem & + np*memory_of_double(nq) &! Delta(np,nq) @@ -176,23 +186,28 @@ END_PROVIDER if (mem > qp_max_mem) then s = s*2.d0 + block_size = block_size / 2 else exit endif if ((s > 1.d0).or.(nq == 0)) then call print_memory_usage() - print *, 'Not enough memory. Reduce cholesky threshold' + print *, 'Required peak memory: ', mem, 'Gb' + call total_memory(mem) + print *, 'Already used memory: ', mem, 'Gb' + 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) +!print *, 'allocate : L(ndim,rank+nq)', memory_of_double(ndim*(rank+nq)) + if (ierr /= 0) then call print_memory_usage() print *, irp_here, ': allocation failed : (L(ndim,rank+nq))' @@ -210,6 +225,8 @@ END_PROVIDER deallocate(L_old) allocate(Delta(np,nq), stat=ierr) +!print *, 'allocate : Delta(np,nq)', memory_of_double(np*nq) + if (ierr /= 0) then call print_memory_usage() print *, irp_here, ': allocation failed : (Delta(np,nq))' @@ -217,6 +234,8 @@ END_PROVIDER endif allocate(Ltmp_p(np,block_size), stat=ierr) +!print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double(np*block_size) + if (ierr /= 0) then call print_memory_usage() print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))' @@ -224,6 +243,8 @@ END_PROVIDER endif allocate(Ltmp_q(nq,block_size), stat=ierr) +!print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double(nq*block_size) + if (ierr /= 0) then call print_memory_usage() print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))' @@ -232,6 +253,7 @@ END_PROVIDER allocate(computed(nq)) +!print *, 'allocate : computed(nq)', memory_of_int(nq) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j) @@ -353,7 +375,7 @@ END_PROVIDER enddo - print '(I10, 4X, ES12.3)', rank, Qmax + print '(I10, 4X, I10, 4X, ES12.3)', rank, block_size, Qmax deallocate(computed) deallocate(Delta) @@ -380,6 +402,8 @@ END_PROVIDER enddo allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) +!print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', ao_num*ao_num*(rank*1_8) * 8_8 / 1024_8**3 + if (ierr /= 0) then call print_memory_usage() print *, irp_here, ': Allocation failed' From 0dca6cfde4ba658518637c2f42ecb2e45c04de6a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 31 May 2024 20:20:04 +0200 Subject: [PATCH 2/5] block size in cholesky --- src/ao_two_e_ints/cholesky.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 09d86679..5a44571c 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -186,7 +186,7 @@ END_PROVIDER if (mem > qp_max_mem) then s = s*2.d0 - block_size = block_size / 2 + block_size = max(block_size / 2, 1) else exit endif From 0a3d462510a6146da21770ff6fff2f7a1794a0fb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 31 May 2024 20:30:48 +0200 Subject: [PATCH 3/5] Clean up openmp in cholesky --- src/ao_two_e_ints/cholesky.irp.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 5a44571c..cfd57050 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -255,7 +255,7 @@ END_PROVIDER allocate(computed(nq)) !print *, 'allocate : computed(nq)', memory_of_int(nq) - !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j) + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q,j) !$OMP DO do q=1,nq @@ -351,7 +351,7 @@ END_PROVIDER ! iii. f = 1.d0/dsqrt(Qmax) - !$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared) + !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared) !$OMP DO do p=1,np Ltmp_p(p,iblock) = Ltmp_p(p,iblock) * f From b743201efe692294db887f175dceb02a81f73422 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 31 May 2024 20:50:30 +0200 Subject: [PATCH 4/5] Use integer*8 in cholesky --- src/ao_two_e_ints/cholesky.irp.f | 152 ++++++++++++++++++------------- src/utils/memory.irp.f | 20 ++++ 2 files changed, 109 insertions(+), 63 deletions(-) diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index cfd57050..daa29079 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -27,19 +27,22 @@ END_PROVIDER ! https://doi.org/10.1016/j.apnum.2011.10.001 : Page 4, Algorithm 1 END_DOC - integer :: rank, ndim + integer*8 :: ndim8 + integer :: rank double precision :: tau double precision, pointer :: L(:,:), L_old(:,:) double precision :: s - double precision, parameter :: dscale = 1.d0 + double precision :: dscale double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:) - integer, allocatable :: Lset(:), Dset(:), addr(:,:) + integer, allocatable :: addr1(:,:), addr2(:,:) + integer**, allocatable :: Lset(:), Dset(:), addr3(:,:) logical, allocatable :: computed(:) - integer :: i,j,k,m,p,q, qj, dj, p2, q2 + integer :: i,j,k,m,p,q, dj, p2, q2 + integer*8 :: i8, j8, p8, qj8 integer :: N, np, nq double precision :: Dmax, Dmin, Qmax, f @@ -47,15 +50,15 @@ END_PROVIDER logical, external :: ao_two_e_integral_zero double precision, external :: ao_two_e_integral - integer :: block_size, iblock, ierr + integer :: block_size, iblock double precision :: mem - double precision, external :: memory_of_double, memory_of_int + double precision, external :: memory_of_double8, memory_of_int8 integer, external :: getUnitAndOpen - integer :: iunit + integer :: iunit, ierr - ndim = ao_num*ao_num + ndim8 = ao_num*ao_num*1_8 deallocate(cholesky_ao) if (read_ao_cholesky) then @@ -83,13 +86,13 @@ END_PROVIDER tau = ao_cholesky_threshold - mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim) + mem = 6.d0 * memory_of_double8(ndim8) + 6.d0 * memory_of_int8(ndim8) call check_mem(mem, irp_here) call print_memory_usage() - allocate(L(ndim,1)) -!print *, 'allocate : (L(ndim,1))', memory_of_double(ndim) + allocate(L(ndim8,1)) +print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8) print *, '' print *, 'Cholesky decomposition of AO integrals' @@ -102,36 +105,36 @@ END_PROVIDER rank = 0 - allocate( D(ndim), Lset(ndim), Dset(ndim) ) - allocate( addr(3,ndim) ) -!print *, 'allocate : (D(ndim))', memory_of_int(ndim) -!print *, 'allocate : (Lset(ndim))', memory_of_int(ndim) -!print *, 'allocate : (Dset(ndim))', memory_of_int(ndim) -!print *, 'allocate : (3,addr(ndim))', memory_of_int(3*ndim) + allocate( D(ndim8), Lset(ndim8), Dset(ndim8) ) + allocate( addr1(ndim8), addr2(ndim8), addr3(ndim8), ) +print *, 'allocate : (D(ndim8))', memory_of_int8(ndim8) +print *, 'allocate : (Lset(ndim8))', memory_of_int8(ndim8) +print *, 'allocate : (Dset(ndim8))', memory_of_int8(ndim8) +print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) ! 1. k=0 do j=1,ao_num do i=1,ao_num k = k+1 - addr(1,k) = i - addr(2,k) = j - addr(3,k) = (i-1)*ao_num + j + addr1(k) = i + addr2(k) = j + addr3(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 - D(i) = ao_two_e_integral(addr(1,i), addr(2,i), & - addr(1,i), addr(2,i)) + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(guided) + do i8=1,ndim8 + D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), & + addr1(i8), addr2(i8)) enddo !$OMP END PARALLEL DO else - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided) - do i=1,ndim - D(i) = get_ao_two_e_integral(addr(1,i), addr(1,i), & - addr(2,i), addr(2,i), & + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(guided) + do i8=1,ndim8 + D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), & + addr2(i8), addr2(i8), & ao_integrals_map) enddo !$OMP END PARALLEL DO @@ -140,12 +143,21 @@ END_PROVIDER Dmax = maxval(D) ! 2. - np=0 - do p=1,ndim - if ( dscale*dscale*Dmax*D(p) > tau*tau ) then - np = np+1 - Lset(np) = p - endif + np = huge(1_4) + dscale = 1.d0 + do while (np == huge(1_4)) + np=0 + do p8=1,ndim8 + if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then + np = np+1 + Lset(np) = p8 + if (np == huge(1_4)) then + ! Overflow detected + dscale = dscale*0.5d0 + exit + endif + endif + enddo enddo ! 3. @@ -155,7 +167,7 @@ END_PROVIDER i = 0 ! 5. - do while ( (Dmax > tau).and.(rank < ndim) ) + do while ( (Dmax > tau).and.(rank < min(ndim8,huge(1_4)) ) ! a. i = i+1 @@ -181,7 +193,8 @@ END_PROVIDER call total_memory(mem) mem = mem & + np*memory_of_double(nq) &! Delta(np,nq) - + (rank+nq)* memory_of_double(ndim) &! L(ndim,rank+nq) + + (rank+nq)* memory_of_double8(ndim8) +&! L(ndim8,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 @@ -205,19 +218,19 @@ END_PROVIDER ! d., e. L_old => L - allocate(L(ndim,rank+nq), stat=ierr) -!print *, 'allocate : L(ndim,rank+nq)', memory_of_double(ndim*(rank+nq)) + allocate(L(ndim8,rank+nq), stat=ierr) +print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq)) if (ierr /= 0) then call print_memory_usage() - print *, irp_here, ': allocation failed : (L(ndim,rank+nq))' + print *, irp_here, ': allocation failed : (L(ndim8,rank+nq))' stop -1 endif !$OMP PARALLEL DO PRIVATE(k,j) do k=1,rank - do j=1,ndim - L(j,k) = L_old(j,k) + do j8=1,ndim8 + L(j8,k) = L_old(j8,k) enddo enddo !$OMP END PARALLEL DO @@ -225,7 +238,7 @@ END_PROVIDER deallocate(L_old) allocate(Delta(np,nq), stat=ierr) -!print *, 'allocate : Delta(np,nq)', memory_of_double(np*nq) +print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8) if (ierr /= 0) then call print_memory_usage() @@ -234,7 +247,7 @@ END_PROVIDER endif allocate(Ltmp_p(np,block_size), stat=ierr) -!print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double(np*block_size) +print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double8(np*block_size*1_8) if (ierr /= 0) then call print_memory_usage() @@ -243,7 +256,7 @@ END_PROVIDER endif allocate(Ltmp_q(nq,block_size), stat=ierr) -!print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double(nq*block_size) +print *, 'allocate : Ltmp_q(nq,block_size)', memory_of_double8(nq*block_size*1_8) if (ierr /= 0) then call print_memory_usage() @@ -253,8 +266,9 @@ END_PROVIDER allocate(computed(nq)) -!print *, 'allocate : computed(nq)', memory_of_int(nq) +print *, 'allocate : computed(nq)', memory_of_int(nq) +print *, 'p1' !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q,j) !$OMP DO @@ -296,7 +310,7 @@ END_PROVIDER iblock = 0 do j=1,nq - if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit + if ( (Qmax <= Dmin).or.(N+j*1_8 > ndim8) ) exit ! i. rank = N+j @@ -308,28 +322,28 @@ END_PROVIDER ! ii. do dj=1,nq - qj = Dset(dj) - if (D(qj) == Qmax) then + qj8 = Dset(dj) + if (D(qj8) == Qmax) then exit endif enddo - L(1:ndim, rank) = 0.d0 + L(1:ndim8, rank) = 0.d0 if (.not.computed(dj)) then m = dj !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided) do k=np,1,-1 - if (.not.ao_two_e_integral_zero( addr(1,Lset(k)), addr(1,Dset(m)),& - addr(2,Lset(k)), addr(2,Dset(m)) ) ) then + if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& + addr2(Lset(k)), addr2(Dset(m)) ) ) then if (do_direct_integrals) then Delta(k,m) = Delta(k,m) + & - ao_two_e_integral(addr(1,Lset(k)), addr(2,Lset(k)),& - addr(1,Dset(m)), addr(2,Dset(m))) + ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& + addr1(Dset(m)), addr2(Dset(m))) else Delta(k,m) = Delta(k,m) + & - get_ao_two_e_integral( addr(1,Lset(k)), addr(1,Dset(m)),& - addr(2,Lset(k)), addr(2,Dset(m)), ao_integrals_map) + get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& + addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map) endif endif enddo @@ -391,18 +405,28 @@ END_PROVIDER Dmax = max(Dmax, D(Lset(p))) enddo - np=0 - do p=1,ndim - if ( dscale*dscale*Dmax*D(p) > tau*tau ) then - np = np+1 - Lset(np) = p - endif + np = huge(1_4) + dscale = 1.d0 + do while (np == huge(1_4)) + np=0 + do p8=1,ndim8 + if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then + np = np+1 + Lset(np) = p8 + if (np == huge(1_4)) then + ! Overflow detected + dscale = dscale*0.5d0 + exit + endif + endif + enddo enddo + enddo allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) -!print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', ao_num*ao_num*(rank*1_8) * 8_8 / 1024_8**3 +print *, 'allocate : cholesky_ao(ao_num,ao_num,rank)', memory_of_double8(ao_num*ao_num*rank*1_8) if (ierr /= 0) then call print_memory_usage() @@ -411,7 +435,9 @@ END_PROVIDER endif !$OMP PARALLEL DO PRIVATE(k) do k=1,rank - call dcopy(ndim, L(1,k), 1, cholesky_ao(1,1,k), 1) + do j=1,ao_num + call dcopy(ao_num, L((j-1)*ao_num+1,k), 1, cholesky_ao(1,j,k), 1) + enddo enddo !$OMP END PARALLEL DO deallocate(L) diff --git a/src/utils/memory.irp.f b/src/utils/memory.irp.f index e69bf71e..043562db 100644 --- a/src/utils/memory.irp.f +++ b/src/utils/memory.irp.f @@ -79,6 +79,26 @@ IRP_ENDIF call unlock_io() end function +double precision function memory_of_double8(n) + implicit none + BEGIN_DOC +! Computes the memory required for n double precision elements in gigabytes. + END_DOC + integer*8, intent(in) :: n + double precision, parameter :: f = 8.d0 / (1024.d0*1024.d0*1024.d0) + memory_of_double8 = dble(n) * f +end function + +double precision function memory_of_int8(n) + implicit none + BEGIN_DOC +! Computes the memory required for n double precision elements in gigabytes. + END_DOC + integer*8, intent(in) :: n + double precision, parameter :: f = 4.d0 / (1024.d0*1024.d0*1024.d0) + memory_of_int8 = dble(n) * f +end function + double precision function memory_of_double(n) implicit none BEGIN_DOC From 38d386d36c78ff87a6c1062d065ab2b12fa4dcc9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 2 Jun 2024 19:03:05 +0200 Subject: [PATCH 5/5] Debug cholesky --- config/gfortran_debug_mkl.cfg | 63 ++++++++++++ src/ao_two_e_ints/cholesky.irp.f | 169 +++++++++++++++++++------------ 2 files changed, 166 insertions(+), 66 deletions(-) create mode 100644 config/gfortran_debug_mkl.cfg diff --git a/config/gfortran_debug_mkl.cfg b/config/gfortran_debug_mkl.cfg new file mode 100644 index 00000000..1dc3f2b2 --- /dev/null +++ b/config/gfortran_debug_mkl.cfg @@ -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 + diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index daa29079..3d0baa48 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -29,21 +29,20 @@ END_PROVIDER integer*8 :: ndim8 integer :: rank - double precision :: tau + double precision :: tau, tau2 double precision, pointer :: L(:,:), L_old(:,:) - double precision :: s - double precision :: dscale + double precision :: dscale, dscale_tmp double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:) - integer, allocatable :: addr1(:,:), addr2(:,:) - integer**, allocatable :: Lset(:), Dset(:), addr3(:,:) + integer, allocatable :: addr1(:), addr2(:) + integer*8, allocatable :: Lset(:), Dset(:), addr3(:) logical, allocatable :: computed(:) integer :: i,j,k,m,p,q, dj, p2, q2 integer*8 :: i8, j8, p8, qj8 - integer :: N, np, nq + integer :: N, np, nq, npmax double precision :: Dmax, Dmin, Qmax, f double precision, external :: get_ao_two_e_integral @@ -53,6 +52,7 @@ END_PROVIDER integer :: block_size, iblock double precision :: mem + double precision, external :: memory_of_double, memory_of_int double precision, external :: memory_of_double8, memory_of_int8 integer, external :: getUnitAndOpen @@ -61,6 +61,9 @@ END_PROVIDER ndim8 = ao_num*ao_num*1_8 deallocate(cholesky_ao) + +! TODO : Save L() to disk + if (read_ao_cholesky) then print *, 'Reading Cholesky vectors from disk...' iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') @@ -85,6 +88,7 @@ END_PROVIDER endif tau = ao_cholesky_threshold + tau2 = tau*tau mem = 6.d0 * memory_of_double8(ndim8) + 6.d0 * memory_of_int8(ndim8) call check_mem(mem, irp_here) @@ -92,25 +96,25 @@ END_PROVIDER call print_memory_usage() allocate(L(ndim8,1)) -print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8) +!print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8) print *, '' print *, 'Cholesky decomposition of AO integrals' print *, '======================================' print *, '' - print *, '============ ============ =============' - print *, ' Rank Block size Threshold' - print *, '============ ============ =============' + print *, '============ =============' + print *, ' Rank Threshold' + print *, '============ =============' rank = 0 allocate( D(ndim8), Lset(ndim8), Dset(ndim8) ) - allocate( addr1(ndim8), addr2(ndim8), addr3(ndim8), ) -print *, 'allocate : (D(ndim8))', memory_of_int8(ndim8) -print *, 'allocate : (Lset(ndim8))', memory_of_int8(ndim8) -print *, 'allocate : (Dset(ndim8))', memory_of_int8(ndim8) -print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) + allocate( addr1(ndim8), addr2(ndim8), addr3(ndim8) ) +!print *, 'allocate : (D(ndim8))', memory_of_int8(ndim8) +!print *, 'allocate : (Lset(ndim8))', memory_of_int8(ndim8) +!print *, 'allocate : (Dset(ndim8))', memory_of_int8(ndim8) +!print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) ! 1. k=0 @@ -124,14 +128,14 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) enddo 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 D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), & addr1(i8), addr2(i8)) enddo !$OMP END PARALLEL DO else - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(guided) + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,16) do i8=1,ndim8 D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), & addr2(i8), addr2(i8), & @@ -143,17 +147,22 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) Dmax = maxval(D) ! 2. - np = huge(1_4) + npmax = huge(1_4)*1_8 + np = npmax dscale = 1.d0 - do while (np == huge(1_4)) + dscale_tmp = Dmax + do while (np == npmax) np=0 do p8=1,ndim8 - if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then + if ( dscale_tmp*D(p8) > tau2 ) then np = np+1 Lset(np) = p8 - if (np == huge(1_4)) then + if (np == npmax) then ! Overflow detected - dscale = dscale*0.5d0 + dscale = dscale*0.1d0 + dscale_tmp = dscale*dscale*Dmax +!print *, 'Overflow detected ' +!print *, 'dscale = ', dscale exit endif endif @@ -167,7 +176,7 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) i = 0 ! 5. - do while ( (Dmax > tau).and.(rank < min(ndim8,huge(1_4)) ) + do while ( (Dmax > tau).and.(rank < min(ndim8,huge(1_4))) ) ! a. i = i+1 @@ -191,15 +200,13 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) call total_memory(mem) - mem = mem & - + np*memory_of_double(nq) &! Delta(np,nq) - + (rank+nq)* memory_of_double8(ndim8) -&! L(ndim8,rank+nq) - + (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size) + mem = mem & + + np*memory_of_double(nq) &! Delta(np,nq) + + (rank+nq)* memory_of_double8(ndim8) &! L(ndim8,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 - block_size = max(block_size / 2, 1) else exit endif @@ -219,7 +226,7 @@ print *, 'allocate : (4,addr(ndim8))', memory_of_int8(4_8*ndim8) L_old => L 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 call print_memory_usage() @@ -227,7 +234,7 @@ print *, 'allocate : L(ndim8,rank+nq)', memory_of_double8(ndim8*(rank+nq)) stop -1 endif - !$OMP PARALLEL DO PRIVATE(k,j) + !$OMP PARALLEL DO PRIVATE(k,j8) do k=1,rank do j8=1,ndim8 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) 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 call print_memory_usage() @@ -247,7 +254,7 @@ print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8) endif 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 call print_memory_usage() @@ -256,7 +263,7 @@ print *, 'allocate : Ltmp_p(np,block_size)', memory_of_double8(np*block_size*1_8 endif 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 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)) -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) +!print *, 'computed' + !$OMP DO + do q=1,nq + computed(q) = .False. + enddo + !$OMP ENDDO NOWAIT + +!print *, 'Delta' !$OMP DO do q=1,nq do j=1,np Delta(j,q) = 0.d0 enddo - computed(q) = .False. enddo !$OMP ENDDO NOWAIT - !$OMP DO +!print *, 'Ltmp_p' do k=1,N + !$OMP DO do p=1,np Ltmp_p(p,k) = L(Lset(p),k) enddo + !$OMP END DO NOWAIT + + !$OMP DO do q=1,nq Ltmp_q(q,k) = L(Dset(q),k) enddo + !$OMP END DO NOWAIT enddo - !$OMP END DO NOWAIT !$OMP BARRIER !$OMP END PARALLEL +!print *, 'p2', np, nq, N if (N>0) then call dgemm('N','T', np, nq, N, -1.d0, & Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) @@ -315,8 +335,10 @@ print *, 'p1' rank = N+j if (iblock == block_size) then +!print *, 'dgemm' call dgemm('N','T',np,nq,block_size,-1.d0, & Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) + iblock = 0 endif @@ -328,36 +350,47 @@ print *, 'p1' endif enddo - L(1:ndim8, rank) = 0.d0 + do i8=1,ndim8 + L(i8, rank) = 0.d0 + enddo if (.not.computed(dj)) then m = dj - !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided) - 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 - if (do_direct_integrals) then - Delta(k,m) = Delta(k,m) + & - ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& - addr1(Dset(m)), addr2(Dset(m))) - else - Delta(k,m) = Delta(k,m) + & - get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map) - endif - endif - enddo - !$OMP END PARALLEL DO + if (do_direct_integrals) then + !$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) + & + ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& + addr1(Dset(m)), addr2(Dset(m))) + endif + enddo + !$OMP END PARALLEL DO + 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) + & + get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& + addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map) + endif + enddo + !$OMP END PARALLEL DO + endif computed(dj) = .True. endif iblock = iblock+1 +!print *, iblock do p=1,np Ltmp_p(p,iblock) = Delta(p,dj) enddo ! iv. 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,& Ltmp_p(1,iblock), 1) endif @@ -365,6 +398,7 @@ print *, 'p1' ! iii. f = 1.d0/dsqrt(Qmax) +!print *, 'p4' !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared) !$OMP DO do p=1,np @@ -379,7 +413,6 @@ print *, 'p1' Ltmp_q(q,iblock) = L(Dset(q), rank) enddo !$OMP END DO - !$OMP END PARALLEL Qmax = D(Dset(1)) @@ -389,12 +422,12 @@ print *, 'p1' 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_q) + deallocate(computed) + deallocate(Delta) ! i. N = rank @@ -405,17 +438,21 @@ print *, 'p1' Dmax = max(Dmax, D(Lset(p))) enddo - np = huge(1_4) + np = npmax dscale = 1.d0 - do while (np == huge(1_4)) + dscale_tmp = Dmax + do while (np == npmax) np=0 do p8=1,ndim8 - if ( dscale*dscale*Dmax*D(p8) > tau*tau ) then + if ( dscale_tmp*D(p8) > tau2 ) then np = np+1 Lset(np) = p8 - if (np == huge(1_4)) then + if (np == npmax) then ! Overflow detected dscale = dscale*0.5d0 + dscale_tmp = dscale*dscale*Dmax +!print *, 'Overflow detected ' +!print *, 'dscale = ', dscale exit endif endif @@ -426,7 +463,7 @@ print *, 'p1' enddo 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 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) do k=1,rank 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 !$OMP END PARALLEL DO