Working on Cholesky CCSD

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -85,13 +85,23 @@ subroutine run_ccsd_space_orb
 
   do while (not_converged)
 
-    call compute_H_oo(nO,nV,t1,t2,tau,H_oo)
-    call compute_H_vv(nO,nV,t1,t2,tau,H_vv)
-    call compute_H_vo(nO,nV,t1,t2,H_vo)
-
     ! Residue
-    call compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+!    if (do_ao_cholesky) then
-    call compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
+    if (.False.) then
+      call compute_H_oo_chol(nO,nV,t1,t2,tau,H_oo)
+      call compute_H_vv_chol(nO,nV,t1,t2,tau,H_vv)
+      call compute_H_vo_chol(nO,nV,t1,t2,H_vo)
+
+      call compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+      call compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
+    else
+      call compute_H_oo(nO,nV,t1,t2,tau,H_oo)
+      call compute_H_vv(nO,nV,t1,t2,tau,H_vv)
+      call compute_H_vo(nO,nV,t1,t2,H_vo)
+
+      call compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+      call compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
+    endif
     max_r = max(max_r1,max_r2)
 
     ! Update
@@ -839,6 +849,10 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
 
 !  allocate(B1(nV,nV,nV,nV))
 !  call compute_B1(nO,nV,t1,t2,B1)
+!  call dgemm('N','N',nO*nO,nV*nV,nV*nV, &
+!              1d0, tau, size(tau,1) * size(tau,2), &
+!                   B1 , size(B1_gam,1) * size(B1_gam,2), &
+!              1d0, r2, size(r2,1) * size(r2,2))
   allocate(B1_gam(nV,nV,nV))
   do gam=1,nV
     call compute_B1_gam(nO,nV,t1,t2,B1_gam,gam)

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -166,11 +166,9 @@ subroutine four_idx_dgemm
 
   deallocate (a1)
 
+  call map_sort(mo_integrals_map)
   call map_unique(mo_integrals_map)
 
-  integer*8                      :: get_mo_map_size, mo_map_size
-  mo_map_size = get_mo_map_size()
-
 end subroutine
 
 subroutine add_integrals_to_map(mask_ijkl)
@@ -250,7 +248,7 @@ subroutine add_integrals_to_map(mask_ijkl)
 
   call wall_time(wall_1)
 
-  size_buffer = min(ao_num*ao_num*ao_num,8000000)
+  size_buffer = min(ao_num*ao_num,8000000)
   print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+&
       ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
 
@@ -443,11 +441,6 @@ subroutine add_integrals_to_map(mask_ijkl)
   !$OMP END PARALLEL
   call map_merge(mo_integrals_map)
 
-  call wall_time(wall_2)
-  call cpu_time(cpu_2)
-  integer*8                      :: get_mo_map_size, mo_map_size
-  mo_map_size = get_mo_map_size()
-
   deallocate(list_ijkl)
 
 
@@ -463,55 +456,55 @@ subroutine add_integrals_to_map_cholesky
 
   integer :: i,j,k,l,m
   integer :: size_buffer, n_integrals
-  size_buffer = min(mo_num*mo_num*mo_num,16000000)
+  size_buffer = min(mo_num*mo_num,16000000)
 
   double precision, allocatable :: Vtmp(:,:,:)
   integer(key_kind)  , allocatable :: buffer_i(:)
   real(integral_kind), allocatable :: buffer_value(:)
 
-  if (.True.) then
+  call set_multiple_levels_omp(.False.)
-    ! In-memory transformation
 
-    call set_multiple_levels_omp(.False.)
+  !$OMP PARALLEL DEFAULT(SHARED) &
+  !$OMP PRIVATE(i,j,k,l,n_integrals,buffer_value, buffer_i, Vtmp)
+  allocate (buffer_i(size_buffer), buffer_value(size_buffer))
+  allocate (Vtmp(mo_num,mo_num,mo_num))
+  n_integrals = 0
 
-    !$OMP PARALLEL PRIVATE(i,j,k,l,n_integrals,buffer_value, buffer_i, Vtmp)
+  !$OMP DO SCHEDULE(dynamic)
-    allocate (buffer_i(size_buffer), buffer_value(size_buffer))
+  do l=1,mo_num
-    n_integrals = 0
+    call dgemm('T','N',mo_num*mo_num,mo_num,cholesky_ao_num,1.d0, &
+       cholesky_mo_transp, cholesky_ao_num, &
+       cholesky_mo_transp(1,1,l), cholesky_ao_num, 0.d0, &
+       Vtmp, mo_num*mo_num)
 
-    allocate (Vtmp(mo_num,mo_num,mo_num))
+    do k=1,l
-
+      do j=1,mo_num
-    !$OMP DO
+        do i=1,j
-    do l=1,mo_num
+          if (dabs(Vtmp(i,j,k)) > mo_integrals_threshold) then
-
+            n_integrals = n_integrals + 1
-      call dgemm('T','N',mo_num*mo_num,mo_num,cholesky_ao_num,1.d0, &
+            buffer_value(n_integrals) = Vtmp(i,j,k)
-         cholesky_mo_transp, cholesky_ao_num, &
+            !DIR$ FORCEINLINE
-         cholesky_mo_transp(1,1,l), cholesky_ao_num, 0.d0, &
+            call mo_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals))
-         Vtmp, mo_num*mo_num)
+            if (n_integrals == size_buffer) then
-      do k=1,l
+              call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
-        do j=1,mo_num
+              n_integrals = 0
-          do i=1,j
-            if (abs(Vtmp(i,j,k)) > mo_integrals_threshold) then
-              n_integrals += 1
-              buffer_value(n_integrals) = Vtmp(i,j,k)
-              !DIR$ FORCEINLINE
-              call mo_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals))
-              if (n_integrals == size_buffer) then
-                call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
-                n_integrals = 0
-              endif
             endif
-          enddo
+          endif
         enddo
       enddo
     enddo
-    !$OMP END DO
+  enddo
+  !$OMP END DO NOWAIT
+
+  if (n_integrals > 0) then
     call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
-    deallocate(buffer_i, buffer_value, Vtmp)
-    !$OMP END PARALLEL
-
-    call map_unique(mo_integrals_map)
-
   endif
+  deallocate(buffer_i, buffer_value, Vtmp)
+  !$OMP BARRIER
+  !$OMP END PARALLEL
+
+  call map_sort(mo_integrals_map)
+  call map_unique(mo_integrals_map)
 
 end
 

src/tools/four_idx_transform.irp.f

@@ -14,6 +14,13 @@ program four_idx_transform
 
   io_mo_two_e_integrals = 'Write'
   SOFT_TOUCH io_mo_two_e_integrals
+  if (.true.) then
+    PROVIDE ao_two_e_integrals_in_map
+  endif
+  if (do_ao_cholesky) then
+    PROVIDE cholesky_mo_transp
+    FREE cholesky_ao
+  endif
   if (.true.) then
     PROVIDE mo_two_e_integrals_in_map
   endif

src/utils/fortran_mmap.c

@@ -22,11 +22,7 @@ void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only)
             perror("Error opening mmap file for reading");
             exit(EXIT_FAILURE);
         }
-        map = mmap(NULL, bytes, PROT_READ, MAP_SHARED | MAP_HUGETLB, fd, 0);
+        map = mmap(NULL, bytes, PROT_READ, MAP_SHARED, fd, 0);
-        if (map == MAP_FAILED) {
-          /* try again without huge pages */
-          map = mmap(NULL, bytes, PROT_READ, MAP_SHARED, fd, 0);
-        }
     }
     else
     {
@@ -53,16 +49,12 @@ void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only)
             exit(EXIT_FAILURE);
         }
 
-        map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_SHARED | MAP_HUGETLB, fd, 0);
+        map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
-        if (map == MAP_FAILED) {
-          /* try again without huge pages */
-          map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
-        }
     }
 
     if (map == MAP_FAILED) {
         close(fd);
-        printf("%s:\n", filename);
+        printf("%s: %lu\n", filename, bytes);
         perror("Error mmapping the file");
         exit(EXIT_FAILURE);
     }

src/utils/mmap.f90

@@ -46,7 +46,14 @@ module mmap_module
       integer(c_size_t)              :: length
       integer(c_int)                 :: fd_
 
-      length = PRODUCT( shape(:) ) * bytes
+      integer :: i
+
+      length = int(bytes,8)
+      do i=1,size(shape)
+        length = length * shape(i)
+      enddo
+print *,  'map_length: ', length
+
       if (read_only) then
           map = c_mmap_fortran( trim(filename)//char(0), length, fd_, 1)
       else
@@ -66,7 +73,13 @@ module mmap_module
       integer(c_size_t)              :: length
       integer(c_int)                 :: fd_
 
-      length = PRODUCT( shape(:) ) * bytes
+      integer :: i
+
+      length = int(bytes,8)
+      do i=1,size(shape)
+        length = length * shape(i)
+      enddo
+print *,  'map_length: ', length
       fd_ = fd
       call c_munmap_fortran( length, fd_, map)
   end subroutine
@@ -82,7 +95,13 @@ module mmap_module
       integer(c_size_t)              :: length
       integer(c_int)                 :: fd_
 
-      length = PRODUCT( shape(:) ) * bytes
+      integer :: i
+
+      length = int(bytes,8)
+      do i=1,size(shape)
+        length = length * shape(i)
+      enddo
+print *,  'map_length: ', length
       fd_ = fd
       call c_msync_fortran( length, fd_, map)
   end subroutine

src/utils_cc/mo_integrals_cc.irp.f

@@ -53,33 +53,8 @@ subroutine gen_v_space(n1,n2,n3,n4,list1,list2,list3,list4,v)
     allocate(v1(cholesky_ao_num,n1,n3), v2(cholesky_ao_num,n2,n4))
     allocate(buffer(n1,n3,n2,n4))
 
-    !$OMP PARALLEL PRIVATE(i1,i2,i3,i4,idx1,idx2,idx3,idx4,k)
+    call gen_v_space_chol(n1,n3,list1,list3,v1,cholesky_ao_num)
-    !$OMP DO
+    call gen_v_space_chol(n2,n4,list2,list4,v2,cholesky_ao_num)
-    do i3=1,n3
-      idx3 = list3(i3)
-      do i1=1,n1
-        idx1 = list1(i1)
-        do k=1,cholesky_ao_num
-          v1(k,i1,i3) = cholesky_mo_transp(k,idx1,idx3)
-        enddo
-      enddo
-    enddo
-    !$OMP END DO NOWAIT
-
-    !$OMP DO
-    do i4=1,n4
-      idx4 = list4(i4)
-      do i2=1,n2
-        idx2 = list2(i2)
-        do k=1,cholesky_ao_num
-          v2(k,i2,i4) = cholesky_mo_transp(k,idx2,idx4)
-        enddo
-      enddo
-    enddo
-    !$OMP END DO NOWAIT
-
-    !$OMP BARRIER
-    !$OMP END PARALLEL
 
     call dgemm('T','N', n1*n3, n2*n4, cholesky_ao_num, 1.d0, &
          v1, cholesky_ao_num, &
@@ -129,6 +104,30 @@ subroutine gen_v_space(n1,n2,n3,n4,list1,list2,list3,list4,v)
 
 end
 
+subroutine gen_v_space_chol(n1,n3,list1,list3,v,ldv)
+
+  implicit none
+
+  integer, intent(in)           :: n1,n3,ldv
+  integer, intent(in)           :: list1(n1),list3(n3)
+  double precision, intent(out) :: v(ldv,n1,n3)
+
+  integer                       :: i1,i3,idx1,idx3,k
+
+  !$OMP PARALLEL DO PRIVATE(i1,i3,idx1,idx3,k)
+  do i3=1,n3
+    idx3 = list3(i3)
+    do i1=1,n1
+      idx1 = list1(i1)
+      do k=1,cholesky_ao_num
+        v(k,i1,i3) = cholesky_mo_transp(k,idx1,idx3)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+end
+
 ! full
 
 BEGIN_PROVIDER [double precision, cc_space_v, (mo_num,mo_num,mo_num,mo_num)]
@@ -345,6 +344,38 @@ BEGIN_PROVIDER [double precision, cc_space_v_vvvv, (cc_nVa, cc_nVa, cc_nVa, cc_n
 
 END_PROVIDER
 
+BEGIN_PROVIDER [double precision, cc_space_v_vv_chol, (cholesky_ao_num, cc_nVa, cc_nVa)]
+
+  implicit none
+
+  call gen_v_space_chol(cc_nVa, cc_nVa, cc_list_vir, cc_list_vir, cc_space_v_vv_chol, cholesky_ao_num)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, cc_space_v_vo_chol, (cholesky_ao_num, cc_nVa, cc_nOa)]
+
+  implicit none
+
+  call gen_v_space_chol(cc_nVa, cc_nOa, cc_list_vir, cc_list_occ, cc_space_v_vo_chol, cholesky_ao_num)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, cc_space_v_ov_chol, (cholesky_ao_num, cc_nOa, cc_nVa)]
+
+  implicit none
+
+  call gen_v_space_chol(cc_nOa, cc_nVa, cc_list_occ, cc_list_vir, cc_space_v_ov_chol, cholesky_ao_num)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, cc_space_v_oo_chol, (cholesky_ao_num, cc_nOa, cc_nOa)]
+
+  implicit none
+
+  call gen_v_space_chol(cc_nOa, cc_nOa, cc_list_occ, cc_list_occ, cc_space_v_oo_chol, cholesky_ao_num)
+
+END_PROVIDER
+
 ! ppqq
 
 BEGIN_PROVIDER [double precision, cc_space_v_ppqq, (cc_n_mo, cc_n_mo)]