Optimization

src/CMakeLists.txt

@@ -17,6 +17,7 @@ set(SRC
     ddx.f90
     ddx_legacy.f90
     ddx_multipolar_solutes.f90
+    ddx_profiling.f90
     llgnew.f
     cbessel.f90
     ddx_cinterface.f90

src/ddx.f90

@@ -390,33 +390,42 @@ subroutine ddrun(ddx_data, state, electrostatics, psi, tol, esolv, &
         return
     end if
 
+    call time_push()
     call setup(ddx_data % params, ddx_data % constants, &
         & ddx_data % workspace, state, electrostatics, psi, ddx_error)
+    call time_pull("setup")
     if (ddx_error % flag .ne. 0) then
         call update_error(ddx_error, "ddrun: setup returned an error, exiting")
         return
     end if
 
     ! solve the primal linear system
     if (do_guess) then
+        call time_push()
         call fill_guess(ddx_data % params, ddx_data % constants, &
             & ddx_data % workspace, state, tol, ddx_error)
+        call time_pull("guess")
         if (ddx_error % flag .ne. 0) then
             call update_error(ddx_error, &
                 & "ddrun: fill_guess returned an error, exiting")
             return
         end if
     end if
+
+    call time_push()
     call solve(ddx_data % params, ddx_data % constants, &
         & ddx_data % workspace, state, tol, ddx_error)
+    call time_pull("solve")
     if (ddx_error % flag .ne. 0) then
         call update_error(ddx_error, "ddrun: solve returned an error, exiting")
         return
     end if
 
     ! compute the energy
+    call time_push()
     call energy(ddx_data % params, ddx_data % constants, &
         & ddx_data % workspace, state, esolv, ddx_error)
+    call time_pull("energy")
     if (ddx_error % flag .ne. 0) then
         call update_error(ddx_error, "ddrun: energy returned an error, exiting")
         return
@@ -425,16 +434,20 @@ subroutine ddrun(ddx_data, state, electrostatics, psi, tol, esolv, &
     ! solve the primal linear system
     if (ddx_data % params % force .eq. 1) then
         if (do_guess) then
+            call time_push()
             call fill_guess_adjoint(ddx_data % params, ddx_data % constants, &
                 & ddx_data % workspace, state, tol, ddx_error)
+            call time_pull("guess adjoint")
             if (ddx_error % flag .ne. 0) then
                 call update_error(ddx_error, &
                     & "ddrun: fill_guess_adjoint returned an error, exiting")
                 return
             end if
         end if
+        call time_push()
         call solve_adjoint(ddx_data % params, ddx_data % constants, &
             & ddx_data % workspace, state, tol, ddx_error)
+        call time_pull("solve adjoint")
         if (ddx_error % flag .ne. 0) then
             call update_error(ddx_error, &
                 & "ddrun: solve_adjoint returned an error, exiting")
@@ -445,8 +458,10 @@ subroutine ddrun(ddx_data, state, electrostatics, psi, tol, esolv, &
     ! compute the forces
     if (ddx_data % params % force .eq. 1) then
         force = zero
+        call time_push()
         call solvation_force_terms(ddx_data % params, ddx_data % constants, &
             & ddx_data % workspace, state, electrostatics, force, ddx_error)
+        call time_pull("solvation forces")
         if (ddx_error % flag .ne. 0) then
             call update_error(ddx_error, &
                 & "ddrun: solvation_force_terms returned an error, exiting")

src/ddx_constants.f90

@@ -223,6 +223,19 @@ module ddx_constants
     !> Whether the diagonal of the matrices has to be used in the mvp for
     !! ddCOSMO, ddPCM or inner ddLPB iterations
     logical  :: dodiag
+    !> List of exposed points in a CSR format.
+    integer, allocatable :: iexposed(:)
+    integer, allocatable :: exposed(:)
+    !> List of buried points in a CSR format.
+    integer, allocatable :: iburied(:)
+    integer, allocatable :: buried(:)
+    !> List of overlapped points in a CSR format.
+    integer, allocatable :: ioverlap(:)
+    integer, allocatable :: overlap(:)
+    !> Adjoint list of overlapped points (the relation is not self
+    ! adjoint) in a CSR format.
+    integer, allocatable :: adj_ioverlap(:)
+    integer, allocatable :: adj_overlap(:)
 end type ddx_constants_type
 
 contains
@@ -423,7 +436,7 @@ subroutine constants_init(params, constants, ddx_error)
         end do
     end if
 
-    ! Generate geometry-related constants (required by the LPB code)
+    ! Generate geometry-related constants
     call constants_geometry_init(params, constants, ddx_error)
     if (ddx_error % flag .ne. 0) then
         call update_error(ddx_error, "constants_geometry_init returned an " // &
@@ -773,9 +786,9 @@ subroutine constants_geometry_init(params, constants, ddx_error)
     !! Local variables
     real(dp) :: swthr, v(3), maxv, ssqv, vv, t
     integer :: i, isph, jsph, inear, igrid, iwork, jwork, lwork, &
-        & old_lwork, icav, info
+        & old_lwork, icav, info, ij, adj_iwork
     integer, allocatable :: work(:, :), tmp_work(:, :)
-    real(dp) :: start_time
+    real(dp) :: start_time, thigh, vij(3), vvij, tij
     !! The code
     ! Prepare FMM structures if needed
     start_time = omp_get_wtime()
@@ -992,6 +1005,98 @@ subroutine constants_geometry_init(params, constants, ddx_error)
             end if
         enddo
     enddo
+
+    allocate(constants % iburied(params % nsph + 1), &
+         & constants % iexposed(params % nsph + 1), &
+         & constants % buried(params % nsph * params % ngrid), &
+         & constants % exposed(params % nsph * params % ngrid), stat=info)
+    if (info .ne. 0) then
+        call update_error(ddx_error, "Buried exposed lists allocation failed")
+        return
+    endif
+    iwork = 1
+    do isph = 1, params % nsph
+        constants % iburied(isph) = iwork
+        do igrid = 1, params % ngrid
+            if (constants % ui(igrid, isph) .lt. one) then
+                !write(7, *) isph, igrid
+                constants % buried(iwork) = igrid
+                iwork = iwork + 1
+            end if
+        end do
+    end do
+    constants % iburied(params % nsph + 1) = iwork
+
+    !!do isph = 1, params % nsph
+    !!    write(6,*) isph, constants % iburied(isph), constants % iburied(isph + 1) - 1
+    !!    do iwork = constants % iburied(isph), constants % iburied(isph + 1) - 1
+    !!        igrid = constants % buried(iwork)
+    !!        write(8, *) isph, igrid
+    !!    end do
+    !!end do
+    !!stop 0
+
+    allocate(constants % ioverlap(constants % inl(params % nsph+1)), &
+        & constants % adj_ioverlap(constants % inl(params % nsph+1)), &
+        & constants % overlap(constants % inl(params % nsph+1)*params % ngrid), &
+        & constants % adj_overlap(constants % inl(params % nsph+1)*params % ngrid), &
+        & stat=info)
+    if (info .ne. 0) then
+        call update_error(ddx_error, "Overlapped grid lists allocation failed")
+        return
+    end if
+    thigh = one + pt5*(params % se + one)*params % eta
+    iwork = 1
+    adj_iwork = 1
+    do isph = 1, params % nsph
+        do ij = constants % inl(isph), constants % inl(isph+1) - 1
+            jsph = constants % nl(ij)
+            constants % ioverlap(ij) = iwork
+            constants % adj_ioverlap(ij) = adj_iwork
+            do igrid = 1, params % ngrid
+                if (constants % ui(igrid, isph) .lt. one) then
+                    vij = params % csph(:,isph) &
+                        & + params % rsph(isph)*constants % cgrid(:,igrid) &
+                        & - params % csph(:,jsph)
+                    vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) &
+                        & + vij(3)*vij(3))
+                    tij  = vvij / params % rsph(jsph)
+                    if (tij.lt.thigh) then
+                        constants % overlap(iwork) = igrid
+                        iwork = iwork + 1
+                    end if
+                end if
+                if (constants % ui(igrid, jsph) .lt. one) then
+                    vij = params % csph(:,jsph) &
+                        & + params % rsph(jsph)*constants % cgrid(:,igrid) &
+                        & - params % csph(:,isph)
+                    vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) &
+                        & + vij(3)*vij(3))
+                    tij  = vvij / params % rsph(isph)
+                    if (tij.lt.thigh) then
+                        constants % adj_overlap(adj_iwork) = igrid
+                        adj_iwork = adj_iwork + 1
+                    end if
+                end if
+            end do
+        end do
+    end do
+    constants % ioverlap(ij) = iwork
+    constants % adj_ioverlap(ij) = adj_iwork
+
+    !do isph = 1, params % nsph
+    !    do ij = constants % inl(isph), constants % inl(isph+1) - 1
+    !        jsph = constants % nl(ij)
+    !        !do ijgrid = constants % ioverlap(ij), constants % ioverlap(ij+1) - 1
+    !        !    igrid = constants % overlap(ijgrid)
+    !        do ijgrid = constants % adj_ioverlap(ij), constants % adj_ioverlap(ij+1) - 1
+    !            igrid = constants % adj_overlap(ijgrid)
+    !            write(6,*) isph, jsph, igrid
+    !        end do
+    !    end do
+    !end do
+    !stop 0
+
     ! Build cavity array. At first get total count for each sphere
     allocate(constants % ncav_sph(params % nsph), stat=info)
     if (info .ne. 0) then
@@ -2144,6 +2249,34 @@ subroutine constants_free(constants, ddx_error)
                 & "deallocation failed!")
         end if
     end if
+    if (allocated(constants % iburied)) then
+        deallocate(constants % iburied, stat=istat)
+        if (istat .ne. 0) then
+            call update_error(ddx_error, "`iburied` " // &
+                & "deallocation failed!")
+        end if
+    end if
+    if (allocated(constants % buried)) then
+        deallocate(constants % buried, stat=istat)
+        if (istat .ne. 0) then
+            call update_error(ddx_error, "`buried` " // &
+                & "deallocation failed!")
+        end if
+    end if
+    if (allocated(constants % iexposed)) then
+        deallocate(constants % iexposed, stat=istat)
+        if (istat .ne. 0) then
+            call update_error(ddx_error, "`iexposed` " // &
+                & "deallocation failed!")
+        end if
+    end if
+    if (allocated(constants % exposed)) then
+        deallocate(constants % exposed, stat=istat)
+        if (istat .ne. 0) then
+            call update_error(ddx_error, "`exposed` " // &
+                & "deallocation failed!")
+        end if
+    end if
 end subroutine constants_free
 
 end module ddx_constants