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