4.2 Smoothed Aggregation

In order to define the prolongator Pk, used to compute the coarse-level matrix Ak+1, MLD2P4 uses the smoothed aggregation algorithm described in [2, 28]. The basic idea of this algorithm is to build a coarse set of indices Ωk+1 by suitably grouping the indices of Ωk into disjoint subsets (aggregates), and to define the coarse-to-fine space transfer operator Pk by applying a suitable smoother to a simple piecewise constant prolongation operator, with the aim of improving the quality of the coarse-space correction.

Three main steps can be identified in the smoothed aggregation procedure:

  1. aggregation of the indices of Ωk to obtain Ωk+1;
  2. construction of the prolongator Pk;
  3. application of Pk and Rk = (Pk)T to build Ak+1.

In order to perform the coarsening step, the smoothed aggregation algorithm described in [28] is used. In this algorithm, each index j Ωk+1 corresponds to an aggregate Ωjk of Ωk, consisting of a suitably chosen index i Ωk and indices that are (usually) contained in a strongly-coupled neighborood of i, i.e.,

               {                 ∘ -------}
Ωkj ⊂ N ki (θ) =  r ∈ Ωk : |akir| > θ |akiiakrr|  ∪ {i} ,
(3)

for a given threshold θ [0, 1] (see [28] for the details). Since this algorithm has a sequential nature, a decoupled version of it is applied, where each processor independently executes the algorithm on the set of indices assigned to it in the initial data distribution. This version is embarrassingly parallel, since it does not require any data communication. On the other hand, it may produce some nonuniform aggregates and is strongly dependent on the number of processors and on the initial partitioning of the matrix A. Nevertheless, this parallel algorithm has been chosen for MLD2P4, since it has been shown to produce good results in practice [5, 7, 27].

The prolongator Pk is built starting from a tentative prolongator Pk nk×nk+1, defined as

                   {              k
P¯k = (¯pkij),  p¯kij =    1    if i ∈ Ω j,
                      0    otherwise,
(4)

where Ωjk is the aggregate of Ωk corresponding to the index j Ωk+1. Pk is obtained by applying to Pk a smoother Sk nk×nk:

Pk = Sk ¯Pk,
in order to remove nonsmooth components from the range of the prolongator, and hence to improve the convergence properties of the multilevel method [2, 26]. A simple choice for Sk is the damped Jacobi smoother:
  k        k   k -1  k
S  =  I - ω (D  )  A F,

where Dk is the diagonal matrix with the same diagonal entries as Ak, A F k = (a ijk) is the filtered matrix defined as

      {   k          k                             ∑
a¯kij =    aij  if j ∈ N i (θ), (j ⁄= i),   ¯akii = akii -   (akij - ¯akij),
         0   otherwise,                            j⁄=i
(5)

and ωk is an approximation of 4(3ρk), where ρk is the spectral radius of (Dk)-1A F k [2]. In MLD2P4 this approximation is obtained by using AF k as an estimate of ρk. Note that for systems coming from uniformly elliptic problems, filtering the matrix Ak has little or no effect, and Ak can be used instead of A F k. The latter choice is the default in MLD2P4.