merged

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
+

src/ao_two_e_ints/cholesky.irp.f

@@ -22,41 +22,51 @@ END_PROVIDER
    ! <ij|kl> = (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
+   integer*8                      :: ndim8
-   double precision               :: tau
+   integer                        :: rank
+   double precision               :: tau, tau2
    double precision, pointer      :: L(:,:), L_old(:,:)
 
-
    double precision               :: s
-   double precision, parameter    :: dscale = 1.d0
+   double precision               :: dscale, dscale_tmp
 
    double precision, allocatable  :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
-   integer, allocatable           :: Lset(:), Dset(:), addr(:,:)
+   integer, allocatable           :: addr1(:), addr2(:)
+   integer*8, 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                        :: N, np, nq
+   integer*8                      :: i8, j8, p8, qj8
+   integer                        :: N, np, nq, npmax
 
    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
+   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
+
+   ndim8 = ao_num*ao_num*1_8
    double precision :: wall0,wall1
 
    call wall_time(wall0)
-   ndim = ao_num*ao_num
    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')
@@ -81,13 +91,15 @@ END_PROVIDER
      endif
 
      tau = ao_cholesky_threshold
+     tau2 = tau*tau
 
-     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))
+     allocate(L(ndim8,1))
+!print *, 'allocate : (L(ndim8,1))', memory_of_double8(ndim8)
 
      print *,  ''
      print *,  'Cholesky decomposition of AO integrals'
@@ -100,32 +112,36 @@ END_PROVIDER
 
      rank = 0
 
-     allocate( D(ndim), Lset(ndim), Dset(ndim) )
+     allocate( D(ndim8), Lset(ndim8), Dset(ndim8) )
-     allocate( addr(3,ndim) )
+     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
+         addr1(k) = i
-         addr(2,k) = j
+         addr2(k) = j
-         addr(3,k) = (i-1)*ao_num + j
+         addr3(k) = (i-1)*ao_num + j
        enddo
      enddo
 
      if (do_direct_integrals) then
-       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided)
+       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,16)
-       do i=1,ndim
+       do i8=1,ndim8
-         D(i) = ao_two_e_integral(addr(1,i), addr(2,i),              &
+         D(i8) = ao_two_e_integral(addr1(i8), addr2(i8),              &
-             addr(1,i), addr(2,i))
+             addr1(i8), addr2(i8))
        enddo
        !$OMP END PARALLEL DO
      else
-       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided)
+       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,16)
-       do i=1,ndim
+       do i8=1,ndim8
-         D(i) = get_ao_two_e_integral(addr(1,i), addr(1,i),          &
+         D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8),          &
-             addr(2,i), addr(2,i),                                   &
+             addr2(i8), addr2(i8),                                   &
              ao_integrals_map)
        enddo
        !$OMP END PARALLEL DO
@@ -134,12 +150,26 @@ END_PROVIDER
      Dmax = maxval(D)
 
      ! 2.
+     npmax = huge(1_4)*1_8
+     np = npmax
+     dscale = 1.d0
+     dscale_tmp = Dmax
+     do while (np == npmax)
        np=0
-     do p=1,ndim
+       do p8=1,ndim8
-       if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
+         if ( dscale_tmp*D(p8) > tau2 ) then
            np = np+1
-         Lset(np) = p
+           Lset(np) = p8
+           if (np == npmax) then
+             ! Overflow detected
+             dscale = dscale*0.1d0
+             dscale_tmp = dscale*dscale*Dmax
+!print *, 'Overflow detected '
+!print *, 'dscale = ', dscale
+             exit
            endif
+         endif
+       enddo
      enddo
 
      ! 3.
@@ -149,13 +179,14 @@ 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
 
-       s = 0.01d0
 
        ! Inrease s until the arrays fit in memory
+       s = 0.01d0
+       block_size = max(N,24)
        do while (.True.)
 
          ! b.
@@ -170,10 +201,11 @@ END_PROVIDER
            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)
+             + (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
@@ -184,6 +216,9 @@ END_PROVIDER
 
          if ((s > 1.d0).or.(nq == 0)) then
            call print_memory_usage()
+           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
@@ -191,20 +226,21 @@ END_PROVIDER
        enddo
 
        ! d., e.
-       block_size = max(N,24)
 
        L_old => L
-       allocate(L(ndim,rank+nq), stat=ierr)
+       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)
+       !$OMP PARALLEL DO PRIVATE(k,j8)
        do k=1,rank
-         do j=1,ndim
+         do j8=1,ndim8
-           L(j,k) = L_old(j,k)
+           L(j8,k) = L_old(j8,k)
          enddo
        enddo
        !$OMP END PARALLEL DO
@@ -212,6 +248,8 @@ END_PROVIDER
        deallocate(L_old)
 
        allocate(Delta(np,nq), stat=ierr)
+!print *, 'allocate : Delta(np,nq)', memory_of_double8(np*nq*1_8)
+
        if (ierr /= 0) then
          call print_memory_usage()
          print *,  irp_here, ': allocation failed : (Delta(np,nq))'
@@ -219,6 +257,8 @@ END_PROVIDER
        endif
 
        allocate(Ltmp_p(np,block_size), stat=ierr)
+!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()
          print *,  irp_here, ': allocation failed : (Ltmp_p(np,block_size))'
@@ -226,6 +266,8 @@ END_PROVIDER
        endif
 
        allocate(Ltmp_q(nq,block_size), stat=ierr)
+!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()
          print *,  irp_here, ': allocation failed : (Ltmp_q(nq,block_size))'
@@ -234,32 +276,47 @@ END_PROVIDER
 
 
        allocate(computed(nq))
+!print *, 'allocate : computed(nq)', memory_of_int(nq)
 
-       !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j)
+!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
+
+!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 k=1,N
-         do p=1,np
-           Ltmp_p(p,k) = L(Lset(p),k)
-         enddo
          do q=1,nq
            Ltmp_q(q,k) = L(Dset(q),k)
          enddo
-       enddo
          !$OMP END DO NOWAIT
+       enddo
 
        !$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)
@@ -276,54 +333,67 @@ 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
 
          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
 
          ! ii.
          do dj=1,nq
-           qj = Dset(dj)
+           qj8 = Dset(dj)
-           if (D(qj) == Qmax) then
+           if (D(qj8) == Qmax) then
              exit
            endif
          enddo
 
-         L(1:ndim, 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( addr(1,Lset(k)), addr(1,Dset(m)),&
-                   addr(2,Lset(k)), addr(2,Dset(m)) ) ) then
            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(addr(1,Lset(k)), addr(2,Lset(k)),&
+                         ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),&
-                     addr(1,Dset(m)), addr(2,Dset(m)))
+                         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)
-               endif
                  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
@@ -331,7 +401,8 @@ END_PROVIDER
          ! iii.
          f = 1.d0/dsqrt(Qmax)
 
-         !$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared)
+!print *, 'p4'
+         !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared)
          !$OMP DO
          do p=1,np
            Ltmp_p(p,iblock) = Ltmp_p(p,iblock) * f
@@ -345,7 +416,6 @@ END_PROVIDER
            Ltmp_q(q,iblock) = L(Dset(q), rank)
          enddo
          !$OMP END DO
-
          !$OMP END PARALLEL
 
          Qmax = D(Dset(1))
@@ -357,10 +427,10 @@ END_PROVIDER
 
        print '(I10, 4X, ES12.3)', rank, Qmax
 
-       deallocate(computed)
-       deallocate(Delta)
        deallocate(Ltmp_p)
        deallocate(Ltmp_q)
+       deallocate(computed)
+       deallocate(Delta)
 
        ! i.
        N = rank
@@ -371,17 +441,33 @@ END_PROVIDER
          Dmax = max(Dmax, D(Lset(p)))
        enddo
 
+       np = npmax
+       dscale = 1.d0
+       dscale_tmp = Dmax
+       do while (np == npmax)
          np=0
-       do p=1,ndim
+         do p8=1,ndim8
-         if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
+           if ( dscale_tmp*D(p8) > tau2 ) then
              np = np+1
-           Lset(np) = p
+             Lset(np) = p8
+             if (np == npmax) then
+               ! Overflow detected
+               dscale = dscale*0.5d0
+               dscale_tmp = dscale*dscale*Dmax
+!print *, 'Overflow detected '
+!print *, 'dscale = ', dscale
+               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)', memory_of_double8(ao_num*ao_num*rank*1_8)
+
      if (ierr /= 0) then
        call print_memory_usage()
        print *,  irp_here, ': Allocation failed'
@@ -389,7 +475,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
+           cholesky_ao(1:ao_num,j,k) = L((j-1)*ao_num+1:j*ao_num,k)
+       enddo
      enddo
      !$OMP END PARALLEL DO
      deallocate(L)

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