Parallelized all the adjoint operations and the FMM l2l and m2m operations.

src/ddx_constants.f90

@@ -219,6 +219,11 @@ module ddx_constants
     !! computation of forces (gradients). Allocated and used only if
     !! fmm=1.
     integer :: grad_nbasis
+    !> Number of layers in the bisection tree. Defined only if fmm=1.
+    integer :: nlayers
+    !> The first and the last clusters of each layer of the tree. The
+    !! dimension is (2, nlayers). Allocated and used only if fmm=1.
+    integer, allocatable :: layers(:, :)
     !> Inner tolerance for microiterations done when using ddLPB
     real(dp) :: inner_tol
     !> Whether the diagonal of the matrices has to be used in the mvp for
@@ -837,7 +842,11 @@ subroutine constants_geometry_init(params, constants)
             & constants % parent, constants % cnode, constants % rnode, &
             & constants % snode, constants % error_message, &
             & constants % error_flag)
-        if (params % error_flag .ne. 0) return
+        if (constants % error_flag .ne. 0) return
+
+        call compute_layers(constants)
+        if (constants % error_flag .ne. 0) return
+
         ! Get number of far and near admissible pairs
         iwork = 0
         jwork = 1
@@ -1429,6 +1438,73 @@ subroutine mkprec(lmax, nbasis, nsph, ngrid, eps, ui, wgrid, vgrid, &
     error_message = ""
 end subroutine mkprec
 
+!> Compute and store information about the layers of the tree used in
+!! the fmm operations.
+!!
+!! @param[inout] constants: Object containing all constants.
+!!
+subroutine compute_layers(constants)
+    implicit none
+    type(ddx_constants_type), intent(inout) :: constants
+    integer :: i, distance, info, max_distance
+    integer, allocatable :: distance_from_root(:)
+
+    allocate(distance_from_root(constants % nclusters), stat=info)
+    if (info.ne.0) then
+        constants % error_flag = 1
+        constants % error_message = "Allocation failed in compute_layers"
+        return
+    end if
+
+    ! first we compute the distances from the root of the tree and
+    ! the number of layers
+    distance_from_root(1) = 0
+    do i = 1, constants % nclusters
+        distance = distance_from_root(i)
+        if (constants % children(1, i) .ne. 0) &
+            & distance_from_root(constants % children(1, i)) = distance + 1
+        if (constants % children(2, i) .ne. 0) &
+            & distance_from_root(constants % children(2, i)) = distance + 1
+        if (distance .gt. max_distance) max_distance = distance
+    end do
+    constants % nlayers = max_distance + 1
+
+    ! now we can allocate the layer information using the actual number
+    ! of layers
+    allocate(constants % layers(2, constants % nlayers), stat=info)
+    if (info.ne.0) then
+        constants % error_flag = 1
+        constants % error_message = "Allocation failed in compute_layers"
+        return
+    end if
+
+    ! finally we can populate the layer information
+    distance = 0
+    constants % layers(1, 1) = 1
+    do i = 1, constants % nclusters
+        if (distance_from_root(i) .lt. distance) then
+            constants % error_flag = 1
+            constants % error_message = "Error in compute_layers"
+            return
+        end if
+        if (distance_from_root(i) .ne. distance) then
+            constants % layers(2, distance + 1) = i
+            if (distance_from_root(i) .lt. constants % nlayers) &
+                & constants % layers(1, distance + 2) = i + 1
+            distance = distance_from_root(i)
+        end if
+    end do
+    constants % layers(2, constants % nlayers) = constants % nclusters
+
+    deallocate(distance_from_root, stat=info)
+    if (info.ne.0) then
+        constants % error_flag = 1
+        constants % error_message = "Deallocation failed in compute_layers"
+        return
+    end if
+
+end subroutine compute_layers
+
 !> Build a recursive inertial binary tree
 !!
 !! Uses inertial bisection in a recursive manner until each leaf node has only
@@ -2249,6 +2325,14 @@ subroutine constants_free(constants)
             return
         end if
     end if
+    if (allocated(constants % layers)) then
+        deallocate(constants % layers, stat=istat)
+        if (istat .ne. 0) then
+            constants % error_message = "`layers` deallocation failed!"
+            constants % error_flag = 1
+            return
+        end if
+    end if
 end subroutine constants_free
 
 end module ddx_constants