In contrast to the standard algebraic multigrid, the Wavelet-based Algebraic Multigrid method relies more strongly on the smoothing method because the coarse spaces are chosen a priori. So, it is very important to develop new smoother methods, especially for those cases where the classical Gauss-Seidel smoothing method does not give good results. This paper proposes a new multilevel smoothing approach based on projection technique. The proposed smoothing method was applied to smoothing the error in a linear systems issued from finite element solutions of the elliptic equation and the results compared with those obtained from the Gauss-Seidel method.

The algebraic multigrid method is a well-grounded numerical technique for solving large sparse linear systems of algebraic equations. The method can be used as an iterative solver or as a preconditioner for one appropriate iterative method. This technique explores the use of some levels to eliminate the high-frequency or oscillatory components of the error, which are not efficiently removed by basic iterative methods, such as Jacobi or Gauss-Seidel [1], [2].

However, the use of the standard multigrid procedure presents some difficulties. The first one is related to the choice of the strength threshold in the coarsening process. The coarsening scheme, for traditional multigrid, can lead to computational complexity growth as the problem size increases, resulting in an elevated memory use and execution time, and in a reduced scalability [3]. Moreover, the coarsening process is inherently sequential in nature that makes difficult its implementation on distributed memory machines [4]. On the other hand, the classical algebraic procedure cannot directly be applied to problems in which the coefficient matrix violates the M-matrix property like, for example, ungauged edge-based finite-element (FE) analysis. In this context, simple approaches such as the use of a shifted coefficient matrix have a limited efficiency [5]. Additionally, the method performs poorly or may break down in elliptic problems with highly discontinuous, oscillatory or anisotropic coefficients [6]-[8].

The investigation of the wavelet decomposition to obtain coarse grid, restriction and prolongation operators, has been the main subject of many multilevel approaches that are aimed to overcome those difficulties [9]-[16]. The similarities between the multigrid methods and wavelets were first brought out by Briggs and Henson [9]. They observed that the space of fine grid vectors in multigrid scheme and the space of highest resolution in the wavelet multiresolution analysis are correlated. Thenceforward, several works have exploited these similarities [10]-[17]. In [15] it is made an exploration of the Briggs's idea to develop wavelet base interpolation and restriction operators in the geometric multigrid context. It compares the results with the geometric multigrid schemes, with conventional intergrid operators, and obtains results similar to the ones obtained with the standard procedure. More recently, other wavelet-based multigrid approaches have shown the efficiency and potentiality of this combination in different areas. Especially in problems related to the computation of electromagnetic fields, the combination of multigrid and wavelets has revealed to be very efficient and promising, presenting in many cases a better performance than some of the most advanced and current traditional multigrid algorithms [19].

As in any multilevel approach, it is very important to ensure an efficient interaction between the smoother and the coarse-grid correction in the wavelet-based algebraic multigrid method (WAMG). In contrast to the standard algebraic multigrid, which uses some simple relaxation scheme and enforces this interaction by choosing the coarser levels and the transfer operators appropriately [1], [21], the WAMG relies more strongly on the smoothing method because the coarse spaces are chosen a priori. So, the WAMG should select a suitable smoothing method to ensure an accurate error representation in the coarsest level. Concerning to this point of view, the WAMG are similar to the traditional geometric multigrid method [1], [21].

In fact, the difficult in the choice of the suitable smoother for geometric multigrid method are not new and many works have already deal with it [22], [23]. In [22], four smoothing methods were compared in the multigrid context to solve the Navier-Stokes equation. Analyses of robustness and efficiency were carried out and the authors concluded that the most robust smoother is also less efficient and has less parallelization potential. Moreover, the paper states that the choice of the smoother method has a strong influence on the robustness and efficiency of a multigrid method and, consequently, further development is necessary. Other works have presented similar conclusions, especially in parallel computation sense [24].

Specifically in the WAMG context and in the computation of electromagnetic fields problems, investigations about the choice of the most appropriate smoother are still very incipient [25]-[27]. In particular, common smoother types do not work properly in the case of the single layer potential [28], [29].

In face of this, it is very important to develop new smoothing methods with suitable parallel properties as an alternative for the classical Gauss-Seidel smoothing method. This paper proposes a new multilevel smoothing approach based on projection techniques which search an approximation in the high frequency subspace created by the high-pass filters in the discrete wavelet transform. The proposed method was applied for smoothing the error in a linear system issued from finite element solutions of the elliptic equation and the results compared with those obtained from the Gauss-Seidel (GS) method. An analysis of the high frequency error after some steps of the smoother is also presented showing the efficiency of the multilevel smoother, mainly for large problems.

II. THE SMOOTHING PROPERTY OF THE STATIONARY ITERATIVE METHODS

The stationary iterative methods like Damped Jacobi and Gauss-Seidel have a smoothing property which is very well exploited by the multigrid technique. Those iterative methods are very efficient for removing the high frequency components of the error while keeping unchanged the smooth components [21].

The convergence of those iterative approaches is defined by the eigenvalues of its iteration matrix. For the Damped Jacobi, in a one-dimensional elliptic problem discretized by Finite Difference Methods, for example, the iteration matrix Pω is given by [1]

in which D is the main diagonal part of A, I is the identity matrix and M = A - D.

The eigenvalues of the matrix Pωfor this example are

where n is the dimension of A.

Then, if h=1/n denotes the mesh spacing:

It means that the first and largest eigenvalue of Pωwill be very close to one when h→0,

As the eigenvalue is related to λ1 an eigenvector in the form

the convergence for Damped Jacobi is slow for the lowest frequency. This property is shared with the others stationary iterative methods [21].

As the smoothing components look like more oscillatory in a coarse grid, the transfer of the problem to a grid with larger mesh size (or larger value of h in equation (3) either a lower value of n in equation (5)) will improve the convergence of the stationary iterative method. Then, if we are sure that the error of our approximation has become smooth after some iteration steps, we can approximate this error by a suitable procedure on a coarser grid [1].

The recursion of this two grid process produces the multigrid method.

III. THE REQUIREMENTS OVER THE SMOOTHER IN WAMG

In order to analyze the requirements over the smoothing methods in WAMG, we have to understand the DWT behavior and the effects of the use of this technique in the algebraic multigrid context. Usually in the literature, the discrete wavelet transform is defined using the concept of filter bank [30]. In such cases low-pass and high-pass filters are chosen to obtain the low frequency and the high frequency terms of the input signal, respectively. In this procedure, if the discrete input signal x = {x[n]}n∈ is sufficiently smooth then the high frequency terms are close to zero and the output low frequency component will be a good approximation of the signal.

To better understand the effects of this technique in the WAMG context, we will suppose that there are only two grids: the original grid and a coarse grid.

The WAMG method uses only the low-pass filters to produce an approximation of the original coefficient matrix. This approximation will represent the system of equation in the coarse grid and it will be used to solve the residual equation on the coarse grid. Then, the residual equation solution is used in the correction of the error in the original grid.

As the WAMG coarse grid is created by the low-pass filters, the coarse grid error correction will only be accurate if the error is smooth. The fact that the error is smooth means that the high frequency components become small after a few stationary iterative iterations whereas the low frequency components hardly change [21]. Therefore, the smoothing method should to be able to remove or, at least, to lessen the high frequency components of the error.

IV. THE PROPOSED MULTILEVEL SMOOTHER METHOD

In a projection process an approximation for the solution of a linear system is extracted from a subspace which is the subspace of candidate approximations, or search subspace. If m is the dimension of , then the residual vector b  Ax is constrained to be orthogonal to m linearly independent vectors in order to obtain such an approximation. This orthogonal condition defines another subspace , of dimension m, called the subspace of constraints [32], [33].

Mathematically, if we know an initial guess, x0∈, the approximate solutions can be defined as

in which r0 = b  Ax0.

As an ideal smoother should remove the high frequency error components it is interesting that this smoother searches an approximation to the linear system solution in a high frequency subspace.

In other hand, the multiresolution analysis in discrete wavelet transform decompose the original grid (space) V in two subspaces V1 and W1 such that [30]

such as

where Φ and ψ are the scaling and wavelets function, respectively. Therefore, we can develop a smoothing method based on a projection approach with,

which will search the approximations to the solution in the high frequency subspace W1. As result the obtained error should be smooth in the sense of the Discrete Wavelet Transform.

In a matrix representation, denoting by W = [w1 , ..., wm] an n x m matrix whose column-vectors form a basis of W1, a two-level smoothing method based on projection technique can be represented by the following algorithm:

Algorithm: Two-level smoothing method

1. Get the current approximation xi

2. rh = b - Axi

3. r2h = WT rh

4. e2h = (WT AW)-1r2h

5. xi+1 = xi + We2h

6. Return xi+1

The recursion of this two-level process produces a new multilevel smoothing method (MS) which has shown to be efficient in smoothing the high frequency error.

V. DESCRIPTION OF THE TEST PROBLEMS

In order to analyze the efficiency of the proposed smoother the method was applied to solve linear systems issued from finite element solutions of the elliptic equation (12), defined in a n × n  point uniform grid on the square domain, with Dirichlet boundary conditions:

Three different discretizations were applied and the resulting matrices are presented in Table I.

For analyzing the WAMG convergence behavior we assume that the exact solution is known (13), and we use it to supply the right hand side vector b.

VI. NUMERICAL RESULTS AND COMMENTS

First of all, the error is analyzed after one, two and three smoother iterations, in order to identify the efficiency of the proposed method in smoothing the high frequency error. So, the smoother was used alone, without the WAMG method. As in the WAMG context is supposed to there is no predefined geometric information about the problem, the error is analyzed in an algebraic sense, what means the error is slowly reduced by the smoothing method [21]. The exact solution presented in (13) also was defined in the algebraic point of view. The corresponding results are presented in Fig. 1. The classical Gauss-Seidel was used for comparison.

It is possible to see from Fig. 1 that the proposed multilevel smoother is more affective in the reduction of the high frequency error for this case.

As a zero initial guess was adopted the initial error is equal to the exact solution .

The numerical reduction of the high frequency components of the error was also evaluated. In this case, the projection of the error in the high frequency subspace W, denoted as We, was analyzed and the results are reported in Fig. 2. Again here, only the application of the smoother has been considered.

As the WAMG method accomplishes the coarse grid correction on the low frequency subspace V, an ideal smoother method should remove the high frequency projection We. The Euclidian norm of the We vector presented in Fig. 2 shows the efficiency of the proposed approach. In this case, the smaller the norm of the We vector the better is the smoothing result.

Finally, the multilevel smoother and the Gauss-Seidel method were compared within the original WAMG algorithm. The WAMG method using the two smoothers was applied to solve the linear systems with the matrices presented in Table I. The variations forward, backward and symmetric of the smoothing method Gauss-Seidel were tested for the WAMG. In the forward (backward) approach the smoothing process is realized in the increasing (decreasing) order of the matrix rows. For the symmetric version, the smoother is applied in an alternate form in even and odd rows. The number of WAMG iterations needed for convergence, and the setup and solver times (in seconds) for each configuration are presented in Table II. The objective of this test was to observe the effect of the reduction of the high frequency error components in the WAMG convergence behavior.

The number of iterations corresponds to the number of V-cycles necessary to reduce the Euclidean norm of the residual vector to the order of 10-6.

From the results in Table II, two main aspects must be highlighted: the multilevel smoother (MS) is computationally more expensive than the Gauss-Seidel (GS) method. It is an expected result because in MS all the levels are visited even though the method was applied only on the first level. Moreover, the WAMG setup time is about two times higher, because it involves the construction of a coarse matrix for the smoother, as defined in step 4 of the algorithm for the high frequency subspace W. However, the reduction number of WAMG iterations means we can get an overall reduction in computational time for large systems. The second aspect is related to the WAMG convergence which preserves the typical mesh-independent convergence behavior of the multigrid cycles.

VII. CONCLUSIONS

The proposed method seems to be very efficient in smoothing the high frequency error. The results in Fig. 1 show that the errors produced by the multilevel approach appear less oscillatory and smaller than those created by the classical Gauss-Seidel method in the three first iterations.

Although being computationally more expensive than the Gauss-Seidel method, the multilevel smoother produces a reduction in the number of iterations that means we could get an overall reduction in computational time. In the tested cases, the WAMG convergence was better when the magnitude of the high frequency components of the error was reduced by the smoother, as presented in Fig. 2. In this aspect the multilevel approach produced the best results, mainly for the larger problem.

Additionally, the multilevel approach has suitable characteristics for parallelization in distributed memory computers, since it can be applied independently in the rows and columns of the matrix exactly in the same way presented in [19] for the WAMG method. So, it can be a good alternative to inherently sequential approaches such as the smoothers based on incomplete LU decompositions, which are robust but has less parallelization potential [22]. The application of the proposed method in the parallel context should be addressed in further work.