Introduction

The main goal of this paper is to spread some computational aspects about the implementation of numerical methods for engineering, computer science and mathematics. The original contribution is not related to the mathematical analysis of the numerical methods for partial differential equations (PDE), but rather to explain the aspects about an efficient and clear computational implementation of a numerical method to approximate the solution of a PDE. We hope that the present work will be used for researchers and people interested in learning this kind of topics. It is important to note that there are papers and books which present in detail the methods that we discuss on this work. Nevertheless, there is not literature concerning to the explanation about the implementation of preconditioning techniques for linear systems associated to PDE.

In the first part of this article we discuss some general aspects related to the H1-conforming Finite Element Method (FEM). The FEM approach is considered as one of the well-established and useful technique for the numerical solution of problems in different fields, described with the use of partial differential equations with different types of boundary conditions. One of the most important aspects related to the success of FEM corresponds to the fact that it is based largely on the basic finite element procedures used: the formulation of the problem in variational form, the discretization of this formulation and the effective solution of the resulting finite element equations. These features do not change from different types of problems (see, e.g., Ciarlet (2002) and Johnson (2009), for details).

In general, the discrete problem for FEM requires to solve a linear system of the form . This kind of systems are usually solved by iterative methods, which approximate the exact solution of the system with a certain number of accurate digits. The corresponding selection of the iterative solver depends on the properties of matrix . For example, for symmetric and positive definite matrices the most popular method is the Conjugate Gradient (CG) method. However, in general, iterative solver depends on the condition number of the matrix of the system, which is related to the performance of the method. According to this, preconditioning techniques can be considered to accelerate convergence of the iterative methods. There are several types of preconditioning techniques, for example, one simple preconditioning technique of the system is to consider the equivalent system for some invertible matrix , which improve the conditioning of the system (see, e.g., Chen (2005)).

An important advance about preconditioning of large linear systems is given by multigrid techniques. The main idea in multigrid approaches is that any iterative solver works on a series of nested meshes, usually defined by a coarsening strategy, which depends on each mesh. Geometric grids are preferred for structured meshes, but algebraic multigrid (AMG) methods are better for unstructured meshes (see, e.g, Stüben (2001)). The AMG requires a simple relaxation process, usually Jacobi or Gauss-Seidel relaxation, and a suitable construction of the coarsening and interpolation operators. The goal of this kind of methods is to define the coarse level variables as a subset of the fine level variables. In other words, to split the fine level variables in two disjoint sets, one of them contains the coarse level variables, and the second one the remaining variables. In practice, the algorithms are usually based on heuristics. For example, for Beck’s algorithm, the coarsening strategy are constructed with the information of the graph associated to the matrix of the system (see Beck (1999)).

Finally, this paper is organized as follows. First, we present a classical FEM scheme for a Poisson equation with homogeneous boundary conditions and describe some aspects about its computational implementation. Next, we present an explanation about preconditioning techniques for primal-FEM formulations. In particular, we focus in the algebraic multigrid preconditioner from Beck (1999), since it is an efficient and simple preconditioner for our variational formulation. However, we include more details related to the construction of the meshes and the coarsening operators. In this way, we describe more aspects than Beck (1999) in order to allow more people to understand the main aspects of the computational implementation of multigrid approaches. In next section, we perform some numerical experiments and discuss the results, to describe the advantages and disadvantages of the preconditioner. Here, we present a complete MATLAB® code. Finally, we provide some conclusions and further work related to the preconditioning techniques for FEM schemes.

H1-conforming finite element method

In this section we present an example of a classical boundary value problem in order to justify the construction of the algebraic multigrid preconditioner. More precisely, given a bounded polygonal domain Ω in IR2 with boundary , we seek a scalar field u such that

(1)

where is a volume force, and as usual . Now, given from the Poisson equation note that

where, employing the Green’s formula in the left-hand side, we arrive at

for all . In other words, we obtain a variational formulation of problem (1), which reads: Find such that

. (2)

It is well known that this problem has a unique solution, which continuously depends of the source .

It is important to note from the new problem (2) that it only requires one order derivative and it establishes a function Hilbert space for the solution. On the other hand, the problem (1) does not allow us the explicit calculation of the solution u, because is an infinite dimensional space. A classical way of circumventing this drawback is the introduction of a finite dimensional subspace of the space defined above consisting of piecewise polynomial continuous functions. More precisely, given be a shape-regular triangulation of Ω without the presence of hanging nodes, and an integer k ≥ 1, it follows that

where, we let be the space of polynomials defined in T of total degree at most k. Furthermore, we define a discrete variational formulation (i.e. Galerkin scheme) which reads: Find such that

. (3)

It is easy to check that the problem (3) has a unique solution (see, for example, Ciarlet (2002)), which in fact is an approximation of the continuous solution of problem (2). Furthermore, there is a basis of , and then

where, at least for k = 1 (low order), is called a Lagrange-type basis, because it holds that , for each vertex of . In other words, the unknowns of the problem (3) correspond to the values of the solution at the vertices of the triangulation of Ω.

In general, the finite element method (FEM) for solving partial differential equations, is one of the high-order discretization schemes that has become a very active research area during the last decade. The main advantage of FEM approaches includes numerical solving of big partial differential systems through the analysis of small versions of the problem (on each element of ). There are some variations of the FEM techniques, for example, the Galerkin scheme (3) is known as the classical finite element technique. On the other hand, there is a locally discontinuous version originally introduced in Cockburn (1998), and a virtual version presented in Beirão da Veiga (2013).

Computational implementation

A fact of the development of FEM methods lies in the difficulty of its computational implementation. Actually, there is few literature detailing the computational implementation of FEM methods, for example, one of the most famous contribution was presented in Carstensen (2002). In the present paper, we do not present details of the implementation of FEM schemes, but in order to explain a preconditioning technique below, we consider some general aspects. Indeed, let T be an element of and define as the Lagrange basis on T, which holds

where are the local degrees of freedom of v on T given by

1.value of v at the ith vertex of T, for all vertex i of T,

2.values of v at k – 1 uniformly spaced points on e, for each edge e of T, for k ≥ 2,

3.values of v at uniformly spaced points in the interior of T, for k ≥ 3.

It is not difficult to see that the local basis is a restriction of on T. Furthermore, the local discrete operators of T that contribute to the global linear system of problem (3) are given by

.

Using these discrete operators, we can assembly the global system of problem (3) following the algorithm:

1.Define and

2.Define and

3.For each do

4. Set up and

5. Define

6. For each vertex i of T

7. p (i) = global index of i on

8. End For

9. For each edge e of T

10. For each j from 1 to k – 1

11. p = (3 + (k – 1) (e – 1) + j)
(number of nodes of )

12. + (k – 1) ((global index of
e on ) – 1) + j

13. End For

14. End For

15. For each from 1 to

16. p (3k + j) = (number of nodes of ) +
(k – 1) (number of edges of ) + j

17. End For

18.

19.

20.End For

Finally, for the incorporation of boundary conditions, we modify the matrix and the vector on the degrees of freedom corresponding to the boundary edges of . The main idea here, is that each row of the linear system associated to the boundary degrees of freedom becomes exactly the boundary condition of problem (1). The incorporation can be done as we described before. However, in order to define a graph-strategy for the construction of an algebraic multigrid for the global linear system, we actually eliminate the degrees of freedom corresponding to the boundary edges of , and then the main unknowns of the system are the values of the solution in the interior nodes of .

Preconditioning techniques

In practice, the linear system associated to problem (3) has a large number of unknowns, then, solving it by Gaussian elimination is quite slow, because this method has a complexity of (see, for details, Saad (2003)). It is well-known (see, for example, Ciarlet (2002)) that is a symmetric positive definite sparse matrix. Then, the Conjugate Gradient method (CG) is the most appropriate to approximate the system solution. However, for the mth iteration of CG holds that

(4)

where is the condition number of . It follows for k >> 1 that the term behaves like , which converges towards 1. Thus, inequality (4) predicts a slow convergence of the CG. On the other hand, the condition number of the matrix of linear system (3) has an asymptotic behavior of , which establishes that k >> 1 when we use very fine triangulations.

According to the previous discussion, it is necessary to use techniques that reduce the condition number of . In general, these techniques construct a nonsingular matrix , which is called a preconditioner, such that the matrix of the equivalent system has a better condition number than the matrix of the original system. Note that if the system is solved immediately. In the case of CG, which requires symmetry in the matrix, the preconditioned system possess the matrix transpose , which is symmetric. The method of solving a linear system by CG using some preconditioner, is known as a Preconditioned Conjugate Gradient (PCG). Some of the most common preconditioners for symmetric systems are Jacobi, Gauss-Seidel, Cholesky and multigrid techniques (see, for example, Saad (2003)).

Multigrid techniques are the most attractive, because they reduce the size of the matrix of the linear system, by a sequence of nested meshes. Consequently, the number of iterations does not depend on the size of the matrix. More precisely, the matrix of the system is reduced by merging or deleting matrix entries. After the new matrix is generated, the corresponding new system is also recursively reduced with the help of the next mesh. Thus, this process is continued and stops when all meshes are used. This sequence of matrices is known as preconditioner levels. For example, in Figure 1 we show a sequence of meshes that satisfy the nested condition. This multigrid preconditioner is known as a geometric multigrid preconditioner (see Trottenberg (2000)).

Once all the preconditioner levels have been generated, some steps of a relaxation method, such as Jacobi or Gauss-Seidel, are performed to solve the linear system. The approximation obtained is restricted to be used as the initial value of the auxiliary linear system given by the next lower level. In the same way this process is followed to the last level. Then, when the + 1 system is solved, the solution is prolonged to level , where again some steps of a relaxation method are performed before moving on to the next level. This whole process is known as the V-cycle, which is illustrated in Figure 2.

Unfortunately, in practice this sequence of nested meshes is not available, only the finest mesh. For this reason, it is necessary to generate this sequence through the fine mesh by agglomerating cells. In Cartesian meshes this process is simpler than in unstructured meshes, where in most cases the agglomeration of cells is not possible. For example, in the case of a triangle mesh, the agglomeration of triangles does not always result in another triangle. An alternative to solve this problem is to create the next level, using the graph that represents the matrix of the system and not by the graph that represents the mesh. This process is known as algebraic multilevel and was introduced by Ruge and Stüben (1987).

Methodology employed

The research proposed in this work considers a scientific approach to a particular problem in the area of Applied Mathematics. It is done under an exploratory quantitative approach. Also, the methods presented are commonly used in problems related to finding approximations for solutions of boundary problems in partial differential equations.

Protocol

1.A particular algebraic multigrid preconditioner is described in detail for the model problem, introduced in the previous section. We consider the method introduced by Beck (1999), which is an interesting work, but it does not show a clear description, making it not accessible for beginning readers. Thus, we intend to explain in more detail some relevant elements of the preconditioner.

2.A proposal for computational implementation for the Beck’s preconditioner (see Beck (1999)) is presented in MATLAB®. The code follows an object-oriented paradigm and, according to its organization, can even be used for the implementation of other multilevel preconditioners different than Beck’s.

3.Once a code has been established in MATLAB® for the Beck’s preconditioner, it is shown its behavior, in order to demonstrate its usefulness. Thus, two model examples of the Poisson problem are considered, whose solution is approximated by FEM scheme. To study the behavior of the preconditioner, only the results related to the resolution of the associated linear systems will be reported. Essentially, the number of iterations performed by the Conjugated Gradient method is shown, when it is used without a preconditioner and with the Beck’s preconditioner.

Additionally, the Incomplete Cholesky method will be considered to compare the behavior between preconditioners. For this purpose, we use meshes that do not exceed 800 000 degrees of freedom. Also, a maximum of 500 iterations and a tolerance relative to the first residue not exceeding 10-6, since these parameters are sufficient to illustrate the desired behavior. On the other hand, it is important to consider one of the examples with a high number of condition to observe that the AMG method retains adequate behavior, which will allow to demonstrate the robustness of this preconditioner.

4.The results are tabulated and plotted to show the good behavior of Beck’s preconditioner over the other methods considered. In particular, the execution time will not be reported because its measurement depends on various factors that are not necessarily controlled by the authors. This should not be construed as a disadvantage of the AMG method, but as a consequence of implementation in software such as MATLAB® and low performance computers.

Beckʼs algebraic multigrid preconditioner

In Ruge and Stüben (1987), Saad (1996), Beck (1999), Castillo and Sequeira (2013), the authors present several algebraic multigrid preconditioners. In particular, in this section we describe an algebraic multigrid technique proposed by Beck (1999), which take advantage of the characteristics of the finite element scheme (3). More precisely, the author identifies the geometric relations of the mesh in the sparse structure of the matrix.

The first algebraic multigrid preconditioner is presented in Ruge and Stüben (1987), which is based on partitioning the variables in each level of the multigrid in two sets by determining strong connection relations imposed by the structure of the matrix associated to the problem. The idea proposed in Beck’s method is a simple coarsening strategy related to the Ruge-Stüben’s method, which define a graph for the matrix using only the sparse structure and ignoring the strength of the connections, i.e. all connections are treated equally (the value of the entries of the matrix are not considered in each partition).

According to the previous discussion, we explain the algebraic multigrid approach using a general linear system:

where and . Now we define and as the data of the first level, where the new right-hand side will be explained later. Thus, the main idea of this method is to reduce to a new matrix with considerably smaller than . Then, this process is repeated in order to generate linear systems with associated matrices such that , for , where N is the maximun number of levels selected such that the last linear system can be solved efficiently.

To perform this procedure for each level we construct matrices and in order to compute . The matrix is known as the restriction operator, whereas the matrix is the prolongation operator. Using these operators, we obtain the matrix from the following algorithm, which define the matrices of each level of the algebraic multigrid preconditioner, summarized in Table 1.

1.For = 1 until N – 1 do

2. Construct and , usually from

3. Define

4.End For

Table 1

Summary of matrices in each level of Beck’s algebraic multigrid preconditioner.

Level	Matrix	Restriction	Prolongation
= 1
= 2
= N – 1
= N		-	-

Note: derived from research.

Finally, in order to describe the construction of each level of the preconditioner (i.e. and ), in the following section we analize a coarsening strategy for sparse matrices, which allow us to construct the restriction and prolongation operators. After that, we introduce the iterative solver based on the V-cycle approach.

Coarsening

We now explain the classical design of the coarsening technique briefly described in Section 3 of Beck (1999). In addition, we incorporate more explanations and examples, including a suggestion for a MATLAB® implementation.

The main idea of the coarsening algorithm is to select a favorably distributed subset of the matrix nodes (i.e. indexes of rows and columns), which are to become the nodes of the new matrix (coarse grid matrix). More precisely, the nodes for the new matrix are known as master nodes, whereas the remaining (which will be dropped) are known as slave nodes. In particular, this coarsening algorithm satisfies that the number of unknowns is reduced substantially in every coarsening step, and the coarse grid matrices preserve the sparsity pattern of the fine grid matrix, that is, coarse grid basis functions retain a local support and a restricted overlap with neighboring ones.

According to Section 3 in Beck (1999), there are three rules for the selection of master nodes. That is:

1.No master node may be directly connected to another master node.

2.There should be as many master nodes as possible.

3.The values of all master nodes are transferred with weight 1. The value for a slave node s is interpolated from the master nodes it is connected to, where each master node contributes with weight . This is used in the construction of the restriction operator below.

Here is important to remember that every Dirichlet node have been eliminated in order that the sparsity pattern of the matrix allows the separation of its nodes into two disjoint sets.

On the other hand, let be a symmetric sparse matrix. Now consider the graph that represents this matrix. More precisely, letting m nodes (numbering from 1 to m), we say that node i is connect to node j, with i ≠ j, if entry . For example, consider matrix

(5)

which graph is given in Figure 3.

Figure 3. Graph that represents the matrix A. Note: derived from research.

Note that the graph follows from the sparsity pattern of A. Now, we select one node to be our first master node. For example, we take node 1 and then, every neighboring node of 1 will be mark as slave node (i.e. node 2, 7 and 8). From the remaining nodes (3, 4, 5, 6 and 9) take the next master node. We take node 3, and then its neighbors (2, 4 and 9) are now slave nodes. Following this procedure, we take node 5 as master node and mark its neighbors as slaves. In Figure 4, we summarize the previous classification strategy.

Figure 4. Master (black) and slave (white) nodes for matrix A. Note: derived from research.

Observe that the coarsening process operates in a geometric way by sequentially choosing a master node and eliminating the corresponding neighboring nodes of the graph. Moreover, a primary criterion for selecting the next master node is to take the node with the minimum number of connections (taking into account the eliminations). And a secondary criterion is to follow the original numbering. We use the first choice by listing the nodes in ascending order, according to the number of non-zero elements of each row of A. We now present a MATLAB® implementation for previous coarsening strategy.

% --------------------------------------------------------

% This function marks the master and slave nodes for a

% given matrix A of order m.

% --------------------------------------------------------

% marked: is a m-size vector which contains 0 on its ith

% component if node i is slave. Otherwise, node

% i is master, and the ith component of the vector

% contains a new numbering for node i.

% sizeM: is the number of master nodes.

% --------------------------------------------------------

function [marked, sizeM] = coarsening(A)

m = length(A); % Order of the matrix

% ----------------------------------------------------

nnzRow = zeros(1, m); % Number of non-zero elements of

% A by row

for i = 1 : m

nnzRow(i) = nnz( A(i,:) );

end

[~, idxRow] = sort(nnzRow); % Sort indices according to

% the number of non-zero elements

% ----------------------------------------------------

% First, every node is marked with a negative number,

% in order to indicate that the node is unmarked.

marked = -1 * ones(1, m);

sizeM = 0;

% Using the symmetry of A, we follow the columns of A:

for j = idxRow

c = A(:,j); % Define the jth column of A

I = find(c); % Find indices of non-zero elements

if marked(j) == -1 % Unmarked node

% Node j has not been marked yet. That’s

% why it will be marked as master and its

% neighboring nodes will be slaves.

% -------------------------------------------

% Mark its neighbors as slaves:

marked(I) = 0;

% Mark as master:

sizeM = sizeM + 1; % A master node was found

marked(j) = sizeM; % New numbering for node j

end

end

end

Using function coarsening() with the matrix defined in (5), we obtain the following results

marked = 1 0 2 0 3 0 0 0 0

sizeM = 3

which establish that the current node 1, 3 and 5 will be, respectively, the node 1, 2 and 3 of the matrix in the next level of the preconditioner.

Restriction and prolongation operators

After choosing master nodes, we can assemble the restriction matrix . Indeed, given the matrix of level ℓ, let be the number of master nodes obtained by the previous coarsening algorithm, such that . Then, letting be the number of neighboring master nodes of slave node s, we define as

Observe that the values of all master nodes are transferred with weight 1, whereas the values of all slave nodes are interpolated from its neighboring master nodes. This selection is strongly related to the definition of the degrees of freedom of scheme (3), at least for k = 1. For example, the restriction matrix for graph in Figure 4 is given by

On the other hand, Beck (1999) suggests a similar construction for prolongation operator . However, in order to preserve symmetry in the matrix of the next level, we define

Finally, using the restriction and prolongation operators we define as the matrix of next level. Then, we continue the process employing to construct and . The stopping criterion corresponds to a maximum order for the last matrix. In addition, the last matrix presents a considerably reduction in the order of the original matrix, then direct methods (for example LU factorization) can be used for solver effects.

According to the above discussion, we now propose the following computational implementation for the construction of all the matrices required by the algebraic multigrid preconditioner.

% ---------------------------------------------------------

% This function constructs the levels of the algebraic

% multigrid preconditioner.

% ---------------------------------------------------------

% A: is a mxm symmetric positive definite sparse matrix.

% maxSize: is the maximum size of the matrix in the last

% level.

% MA: is a list of matrices for coarse matrices

% MR: is a list of matrices for restriction operators

% MP: is a list of matrices for prolongation operators

% L: is the lower triangular matrix of the LU factorization

% of the last matrix A

% U: is the upper triangular matrix of the LU factorization

% of the last matrix A

% steps: is a list of the Smoothing steps for Gauss-Seidel

% ---------------------------------------------------------

function [MA, MR, MP, L, U, steps] = createLevels(A, maxSize)

MA = {}; % List of matrices

MR = {}; % List of restriction operators

MP = {}; % List of prolongation operators

% ----------------- Create the levels -----------------

numLevels = 0; % Number of levels

m = length(A); % Order of the matrix

while m >= maxSize

% Mark nodes as Master and Slave

[marked, sizeM] = coarsening(A);

% Construct Restriction matrix

R = sparse(sizeM, m);

for j = 1 : m

if marked(j) == 0

% -------------- SLAVE NODE ---------------

c = A(:,j); % Define the jth column of A

I = find(c); % Find indices of non-zero

% elements of the jth column of A

[~, ~, I] = find( marked(I) ); % Find indices

% of master nodes

% arround node j

nnz = length(I); % Number of neighboring master

% nodes of node j

R(I, j) = 1 / nnz;

else

% -------------- MASTER NODE --------------

R(marked(j), j) = 1;

end

end

% Store matrices on lists

numLevels = numLevels + 1;

MR{numLevels} = R;

MA{numLevels} = A;

MP{numLevels} = R’;

% Preparations for next level

A = R * A * R’;

m = sizeM;

end

% ---- Compute LU Factorization for the last matrix ---

[L,U] = lu(A);

% ---------- Create list of smoothing steps -----------

steps = zeros(numLevels, 1);

for i = 1 : numLevels

steps(i) = i + 1;

end

end

The vector steps, defined by function createLevels(), contains the number of relaxation steps for the V-cycles of the preconditioner. For simplicity, we perform ℓ + 1 relaxation steps for level ℓ.

V-cycles iteration

According to the previous sections, we now aim to present an iterative solver for linear system based on the iterative refinement method (see, e.g, Wilkinson (1994) or Burden, Faires and Burden (2015)). Indeed, given, , and and an initial guess, we consider the following algorithm

1.Define

2.For k = 1 until number of V-cycles do

3. For = 1 until N – 1 do

4. Define

5. Perform some steps of a relaxation
method to the system

6. Define

7. End For

8. Solve by a direct method

9. For = N – 1 until 1 do

10. Define

11. Perform some steps of a relaxation
method to the system

12. End For

13.End For

Here N is the total number of levels defined by createLevels() in the previous section. Notice that, after the calculation of the first residual and according to the iterative refinement method, we need to solve alternative linear system in order to make the correction . Thus, for solving we perform some steps of a relaxation method, such as the Jacobi iteration or Gauss-Siedel iteration. Next, we restrict the new residual to the following level. We continue this process until level N – 1. In last level, we compute the solution of the corresponding system employing a direct solver. In particular, for simplicity we use LU Factorization, but Cholesky Factorization is also a good choice. Finally, through each level, the solution is prolonged and updated until a new approximation to the original system is obtained.

It is important to observe that the original residual is not changed. In addition, in order to preverse symmetry (see, e.g., Chen (2005)), we utilize Forward Gauss-Seidel iteration when each residual is restricted. On the other hand, we use Backward Gauss-Seidel iteration when each residual is prolonged. We end this section by illustrating the definition of every V-cycle in Figure 5. The MATLAB® code for this solver will be presented in function apply() in a section below.

Preconditioned Conjugate Gradient Method

The V-cycle iteration is a solver for the linear system . However, this is used as a preconditioner for the usual preconditioned conjugate gradient (PCG) method. More precisely, consider the PCG method (see, e.g., Saad (2003)) given by

1.Define

2.Apply preconditioner to solve

3.Define

4.Define

5.For k = 0 until convergence do

6. Define

7. Define

8. Define

9. Update

10. Update

11. Apply preconditioner to solve

12. Define

13. Update

14.End For

Next, note from steps 2 and 11 that it is necessary to solve linear systems for every preconditioner. In this case, we will use the iterative solver defined in previous section.

Complete MATLAB® code

In this section, we propose a simplified computational implementation for the algebraic multigrid preconditioner introduced by Beck (1999). All procedures described in the foregoing sections are unified in a MATLAB® class called amgBR, which is presents at end of this paper (see Annex 1).

On the other hand, consider the following implementation for PCG method:

% ---------------------------------------------------------

% Preconditioned conjugate gradient method.

% ---------------------------------------------------------

% P: is a preconditioner

% A: is the matrix of the linear system

% b: is the right-hand side of the linear system

% x: is an initial guess

% tol: tolerance

% iterMax: is the maximum number of iterations

% iterUsed: is the number of iterations used

% ---------------------------------------------------------

function [x, iterUsed] = pcg(P, A, b, x, tol, iterMax)

r = b - A*x;

z = P.apply(r);

p = z;

nrm0Tol = tol * norm(r);

iterUsed = 0;

while iterUsed < iterMax && norm(r) >= nrm0Tol

v = A * p;

w = r’*z;

alpha = w / (v’*p);

x = x + alpha * p;

r = r - alpha * v;

z = P.apply(r);

beta = (r’*z) / w;

p = z + beta * p;

iterUsed = iterUsed + 1;

end

end

Then, the class amgRB, defined at the beginning of this section, can be used as in the following example:

P = amgRB(A, maxSize, numVCycles); % Define preconditioner

x0 = zeros(length(A),1); % Define initial guess

[x, k] = pcg(P, A, b, x0, 1.0e-6, 1000)

where A, b, maxSize, and numVCycles must be previously defined.

Numerical experiments

In this section, we carry out two numerical experiments to validate our code and illustrate the behavior of the algebraic multigrid preconditioner proposed in Beck (1999) at least for k = 1. For all the examples, triangular meshes were generated using the software TRIANGLE developed by Shewchuk (1996). In addition, we modify the continuous problem (1) to

where is a uniformly positive definite tensor describing the material properties of Ω. The idea of introducing the tensor is to consider a problem with high condition number (see Example 2 below).

The numerical results were obtained using the MATLAB® code introduced in previous section. The corresponding linear systems are solved by the Conjugate-Gradient method with a tolerance of 10–6, 500 as the maximum number of iterations, and taking as initial guess the vector which each entry is equal to one. Here, we use ndofs for the number of unknown of each linear system, that is, the value of . In addition, we introduce the integer in order to define the number of smoothing steps used by Gauss-Seidel methods. More precisely, for a level , the number of relaxation steps are given by . Moreover, we only use one V-Cycle for the AMG preconditioner. Finally, we also consider the classical Incomplete Cholesky preconditioner (see, e.g, Saad (2003)) as an alternative solver.

In Example 1, we consider , the identity matrix, and choose the data so that the exact solution is given by

for all . In Figure 6, we show a triangular decomposition for our domain obtained with TRIANGLE, where it is important to note that only unstructured meshes are used.

Figure 6. Mesh for Example 1 with 1531 triangles. Note: derived from research.

In Table 2, the results of the AMG preconditioner studied in this work are presented and compared to the classical CG method. We also obtain results for the Incomplete Cholesky method. The AMG method shows an improvement with respect to the classical methods, in terms of the number of iterations and the obtained residuals. For instance, note that the CG method almost never converges after 500 iterations for any of the experiments, while the Incomplete Cholesky method, although converging, the number of iterations is significantly increasing compare to AMG. In fact, the most important aspect to observe is that the AMG method always converges, and its number of iterations does not exceed 22 (considering all experiments), which allow us to see the robustness of the preconditioner. In addition, these iterations decrease as the parameter increases, which establishes that several relaxation steps carry out fewer iterations, but will require more execution time, which is not shown in this paper because the measurement of time depends on multiple reasons related to particular situations.

On the other hand, for the last mesh of 782939 degrees of freedom, we show in Table 3 a description of each level created by the AMG preconditioner. There, it is possible to see a significant reduction in the unknowns of the system. Finally, in Figure 7, we show that the residuals of the AMG method decreases faster.

Table 3

Example 1: Description of the levels of the AMG preconditioner for the mesh of 782939 degrees of freedom.

Level	Order of the matrix	Number of non-zero entries
1	782939	5472891
2	207841	2001147
3	42911	478627
4	7873	94105
5	1364	16372
6	237	6010

Note: derived from research.

In Example 2, we consider the model problem on a domain with two different materials and an isotropic diffusion. More precisely, Ω consists of two squares, the first one is the unit square (0,1)2, which has diffusion equal to the identity matrix , and within this square, a second square with diffusion . An example of the meshes for this domain is shown in Figure 8. Here the source function ℓ is equal to zero.

Figure 8. Mesh for Example 2 with 1543 triangles. Note: derived from research.

In Table 3, we note that the AMG method performs a behavior similar to the obtained in Example 1, although the condition number of the system matrix in Example 2 is 109 greater than the condition number of Example 1. Thus, the AMG preconditioner does not show a significant change by increasing the conditioning of the matrix. Finally, in Table 4, a level description is presented, whereas in Figure 9, we show the residuals history, all this for the last mesh used in Table 4.

Table 5

Example 2: Description of the levels of the AMG preconditioner for the mesh of 781506 degrees of freedom.

Level	Order of the matrix	Number of non-zero entries
1	781506	5458246
2	208021	2003307
3	42978	477048
4	7909	93343
5	1404	16444
6	250	4927

Note: derived from research.

Figure 9. Example 2: Residual history for 100 iterations (left) and 20 iterations (right). Note: derived from research.

Conclusions and future directions

In this article, we have described an object oriented implementation of the algebraic multigrid preconditioner introduced in Beck (1999), applied to elliptic problems using unstructured triangular meshes. The numerical experiments presented illustrate the appropriate behavior of the preconditioner. In particular, there is not a considerable increase in the number of iterations of PCG method even when the mesh size is decreased.

In the future we would like to consider a new coarsening strategy, which allow us a better selection for master nodes. In addition, in this work we used only H1-conforming elements in order to explain the preconditioner in a cleaner way. However, the results presented here can be extended to the case of H(div)-conforming elements. Finally, we are interesting to develop an extension of the Beck’s preconditioner to polynomial bases of degree higher than one (that is, k > 1), which implies that the degrees of freedom are not only in the vertices, and then the coarsening technique would not be accurate.

Acknowledgements

The work of Jeremías Ramírez was partially supported by Escuela de Matemática of Universidad Nacional, Costa Rica, through the project 0161-17. The work of Esteban Segura was partially supported by Escuela de Matemática and CIMPA, Universidad de Costa Rica, Costa Rica, through the project 821-B7-254. The work of Filánder Sequeira was partially supported by Escuela de Matemática of Universidad Nacional, Costa Rica, through the project 0103-18.

References

Beck, R. (1999). Graph-based algebraic multigrid for Lagrange-type finite elements on simplicial meshes. Berlin, Germany: Konrad Zuse Zentrum für Informationstechnik. Recuperado de https://pdfs.semanticscholar.org/c719/beb475ff33cab3f235c361f7e8d76e9f820c.pdf

Beirão da Veiga, L.; Brezzi, F.; Cangiani, A.; Manzini, G.; Marini, L. D. & Russo, A. (2013). Basic principles of virtual element methods. Mathematical Models and Methods in Applied Sciences, 23(01), 199-214. doi: https://doi.org/10.1142/S0218202512500492

Burden, R. L.; Faires, J. D. & Burden, A. M. (2015). Numerical Analysis (10th ed.). United States: Cengage Learning.

Carstensen, C. & Klose, R. (2002). Elastoviscoplastic finite element analysis in 100 lines of Matlab. Journal of Numerical Mathematics, 10(3), 157-192. doi: https://doi.org/10.1515/JNMA.2002.157

Castillo, P. E. & Sequeira, F. A. (2013). Computational aspects of the Local Discontinuous Galerkin method on unstructured grids in three dimensions. Mathematical and Computer Modelling, 57(9-10), 2279-2288. doi: https://doi.org/10.1016/j.mcm.2011.07.032

Chen, K. (2005). Matrix preconditioning techniques and applications. United Kingdom: Cambridge University Press. doi: https://doi.org/10.1017/CBO9780511543258

Ciarlet, P. G. (2002). The finite element method for elliptic problems. United States: SIAM. doi: https://doi.org/10.1137/1.9780898719208

Cockburn, B. & Shu, C. W. (1998). The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6), 2440-2463. doi: https://doi.org/10.1137/S0036142997316712

Johnson, C. (2009). Numerical solution of partial differential equations by the finite element method. United States: Courier Corporation.

Ruge, J. W. & Stüben, K. (1987). Algebraic multigrid. In: S.F., McCormick. (Ed.), Multigrid Methods (pp. 73-130). United States: SIAM. doi: https://doi.org/10.1137/1.9781611971057.ch4

Saad, Y. (1996). ILUM: a multi-elimination ILU preconditioner for general sparse matrices. SIAM Journal on Scientific Computing, 17(4), 830-847. doi: https://doi.org/10.1137/0917054

Saad, Y. (2003). Iterative methods for sparse linear systems (2th ed.). United States: SIAM. Recuperado de https://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf

Shewchuk, J. R. (1996). Triangle: engineering a 2D quality mesh generator and Delaunay triangulator. In Workshop on Applied Computational Geometry (pp. 203-222). Heidelberg, Germany: Springer. doi: https://doi.org/10.1007/BFb0014497

Stüben, K. (2001). A review of algebraic multigrid. In: C., Brezinski, L., Wuytack. (Eds.), Numerical Analysis: historical Developments in the 20th Century (pp. 331-359). Elsevier. doi: https://doi.org/10.1016/B978-0-444-50617-7.50015-X

Trottenberg, U.; Oosterlee, C. W. and Schuller, A. (2000). Multigrid. United States: Academic Press.

Wilkinson, J. H. (1994). Rounding errors in algebraic processes. United States: Courier Corporation.

ANNEX 1: MATLAB Class for algebraic multigrid preconditioner

% -------------------------------------------------------------

% Class for the Algebraic Multigrid Preconditioner

% introduced by Rudolf Beck

% -------------------------------------------------------------

classdef amgRB

properties (GetAccess=private)

maxSize = 100; % Default values

numVCycles = 5; % Default values

A = []; % Matrix for first level

MA = {}; % List of matrices

MR = {}; % List of restriction operators

MP = {}; % List of prolongation operators

L = []; % For LU factorization of the last matrix

U = []; % For LU factorization of the last matrix

tol = 1e-6; % Tolerance for the iterative solver

steps = []; % Smoothing steps for Gauss-Seidel (GS)

ML = {}; % List of lower triangular matrices for GS

MU = {}; % List of upper triangular matrices for GS

MFGS = {}; % List of inverse matrices for Forward-GS

MBGS = {}; % List of inverse matrices for Backward-GS

end

% ---------------------------------------------------------

% PUBLIC METHODS

% ---------------------------------------------------------

methods

% ---------------------------------------------------------

% Construct method

% This function constructs the levels of the algebraic

% multigrid preconditioner.

% ---------------------------------------------------------

% Receive:

% A: is a mxm symmetric positive definite sparse matrix.

% maxSize: is the maximum size of the matrix in the last

% level.

% numVCycles: is the number of V-cycles.

% Update:

% MA: is a list of matrices for coarse matrices

% MR: is a list of matrices for restriction operators

% MP: is a list of matrices for prolongation operators

% L: is the lower triangular matrix of the LU factorization

% of the last matrix A

% U: is the upper triangular matrix of the LU factorization

% of the last matrix A

% steps: is a list of the smoothing steps for Gauss-Seidel

% ML: is a list of lower triangular matrices for GS

% MU: is a list of upper triangular matrices for GS

% MFGS: is a list of inverse matrices for Forward-GS

% MBGS: is a list of inverse matrices for Backward-GS

% ---------------------------------------------------------

function obj = amgRB(A, maxSize, numVCycles)

obj.A = A;

obj.maxSize = maxSize;

obj.numVCycles = numVCycles;

% ----------------- Create the levels -----------------

numLevels = 0; % Number of levels without the last one

m = length(A); % Order of the matrix

while m >= maxSize

% Mark nodes as Master and Slave

[marked, sizeM] = coarsening(obj, A);

% Construct Restriction matrix

R = sparse(sizeM, m);

for j = 1 : m

if marked(j) == 0

% -------------- SLAVE NODE ---------------

c = A(:,j); % Define the jth column of A

I = find(c); % Find indices of non-zero

% elements of the jth column of A

[~,~,I] = find( marked(I) ); % Find indices

% of master nodes

% around node j

nnz = length(I); % Number of neighboring

% master nodes of node j

R(I, j) = 1 / nnz;

else

% -------------- MASTER NODE --------------

R(marked(j), j) = 1;

end

end

% Store matrices on lists

numLevels = numLevels + 1;

obj.MR{numLevels} = R;

obj.MA{numLevels} = A;

obj.MP{numLevels} = R’;

% Preparations for next level

A = R * A * R’;

m = sizeM;

end

% ---- Compute LU Factorization for the last matrix ---

[L,U] = lu(A);

obj.L = L;

obj.U = U;

% ---------- Create list of smoothing steps -----------

obj.steps = zeros(numLevels, 1);

for i = 1 : numLevels

obj.steps(i) = i + 1;

end

% ---------------- Setup Gauss-Seidel -----------------

obj.ML{numLevels} = []; % Store memory

obj.MU{numLevels} = []; % Store memory

obj.MFGS{numLevels} = []; % Store memory

obj.MBGS{numLevels} = []; % Store memory

% Create Gauss-Seidel matrices

for l = 1 : numLevels

A = obj.MA{l}; % Get matrix A_l

I = speye(size(A)); % Create identity matrix

obj.ML{l} = tril(A,-1); % Lower triangular part of A

obj.MU{l} = triu(A, 1); % Upper triangular part of A

obj.MFGS{l} = tril(A) \ I; % Forward-GS matrix

obj.MBGS{l} = triu(A) \ I; % Backward-GS matrix

end

end

% ---------------------------------------------------------

% Get the number of levels of the preconditioner

% ---------------------------------------------------------

% k: is the number of levels.

% ---------------------------------------------------------

function [k] = getNumLevels(obj)

k = length(obj.MR) + 1;

end

% ---------------------------------------------------------

% Apply preconditioner to a vector, used by PCG method

% ---------------------------------------------------------

% x, y: are m-size vectors

% ---------------------------------------------------------

function [y] = apply(obj, x)

y = zeros( size(x) );

[y, ~, ~] = solve(obj, x, y);

end

% ---------------------------------------------------------

% Iterative solver for Ax = b

% ---------------------------------------------------------

% b: is the right-hand side of the system

% x0: is an initial guess for system Ax = b

% x: is the new approximation for system Ax = b

% cyclesUsed: is the number of V-cycles used

% normResidual: is a list for the norm-2 of the

% residual for each V-cycle

% ---------------------------------------------------------

function [x, cyclesUsed, normResidual] = solve(obj, b, x0)

numLevels = length(obj.MR); % Number of levels without

% the last one, that is, N-1

% ---------- When the number of levels is 1 -----------

if numLevels == 0

% Solve system using LU Factorization

x = obj.U \ ( obj.L \ b );

cyclesUsed = 0;

normResidual(1) = norm(b - obj.A*x);

return % Leave the function

end

% ------------- Compute the first residual ------------

A = obj.A; % Matrix of the first level

r = b - A*x0; % First residual

nrmR = norm(r); % Norm-2 of the residual

normResidual(1) = nrmR;

nrm0Tol = obj.tol * nrmR; % Tolerance relative to the

% first residual

% ------------------- Star V-cycles -------------------

cyclesUsed = 0;

while cyclesUsed < obj.numVCycles && nrmR >= nrm0Tol

% ------------ Go downward in V-cycle -------------

for l = 1 : numLevels

A = obj.MA{l}; % Get matrix A_l

x = zeros( size(r) ); % Create solution x_l

% ---------------------------------------------

x = forwardGaussSeidel(obj, r, x, l);

% ---------------------------------------------

% Stores new vectors

approxim{l} = x;

residual{l} = r;

% Restrict residual

r = obj.MR{l} * ( r - A*x );

end

% --- Perform LU Factorizacion in the last level --

approxim{numLevels+1} = obj.U \ ( obj.L \ r );

% ------------- Go upward in V-cycle --------------

for l = numLevels : -1 : 1

A = obj.MA{l}; % Get matrix A_l

x = approxim{l}; % Get vector x_l

r = residual{l}; % Get vector r_l

% Prolonging residual

x = x + obj.MP{l} * approxim{l + 1};

% ---------------------------------------------

x = backwardGaussSeidel(obj, r, x, l);

% ---------------------------------------------

approxim{l} = x; % Stores new approximation

end

% ----------------- End of V-cycle ----------------

r0 = b - A*x; % Create new residual for the first

% level, but this does not modify the

% original

nrmR = norm(r0);

normResidual(cyclesUsed + 2) = nrmR;

cyclesUsed = cyclesUsed + 1;

end

end

end % End Public Methods

% ---------------------------------------------------------

% PRIVATE METHODS

% ---------------------------------------------------------

methods (Access = private)

% ---------------------------------------------------------

% This function marks the master and slave nodes for a

% given matrix A of order m.

% ---------------------------------------------------------

% marked: is a m-size vector which contains 0 on its ith

% component if node i is slave. Otherwise, node

% i is master, and the ith component of the vector

% contains a new numbering for node i.

% sizeM: is the number of master nodes.

% ---------------------------------------------------------

function [marked, sizeM] = coarsening(obj, A)

m = length(A); % Order of the matrix

% -----------------------------------------------------

nnzRow = zeros(1, m); % Number of non-zero elements of

% A by row

for i = 1 : m

nnzRow(i) = nnz( A(i,:) );

end

[~, idxRow] = sort(nnzRow); % Sort indices according to

% the number of non-zero elements

% -----------------------------------------------------

% First, every node is marked with a negative number,

% in order to indicate that the node is unmarked.

marked = -1 * ones(1, m);

sizeM = 0;

% Using the symmetry of A, we follow the columns of A:

for j = idxRow

c = A(:,j); % Define the jth column of A

I = find(c); % Find indices of non-zero elements

if marked(j) == -1 % Unmarked node

% Node j has not been marked yet. That’s

% why it will be marked as master and its

% neighboring nodes will be slaves.

% ---------------------------------------------

% Mark its neighbors as slaves:

marked(I) = 0;

% Mark as master:

sizeM = sizeM + 1; % A master node was found

marked(j) = sizeM; % New numbering for node j

end

end

end

% ---------------------------------------------------------

% Forward Gauss-Seidel iteration.

% ---------------------------------------------------------

% b: is the right-hand side vector.

% x: is the initial vector.

% l: is the corresponding level.

% ---------------------------------------------------------

function x = forwardGaussSeidel(obj, b, x, l)

numSteps = obj.steps(l);

for k = 1 : numSteps

x = obj.MFGS{l} * ( b - obj.MU{l}*x );

end

end

% ---------------------------------------------------------

% Backward Gauss-Seidel iteration.

% ---------------------------------------------------------

% b: is the right-hand side vector.

% x: is the initial vector.

% l: is the corresponding level.

% ---------------------------------------------------------

function x = backwardGaussSeidel(obj, b, x, l)

numSteps = obj.steps(l);

for k = 1 : numSteps

x = obj.MBGS{l} * ( b - obj.ML{l}*x );

end

end

end % End Private Methods

end

% ------------------------------------------------------------

% END OF THE FILE

% ------------------------------------------------------------