Adaptive GMRES(m) for the Electromagnetic Scattering Problem

ABSTRACT

In this article, an adaptive version of the restarted GMRES (GMRES(m)) is introduced for the resolution of the finite difference approximation of the Helmholtz equation. It has been observed that the choice of the restart parameter m strongly affects the convergence of standard GMRES(m). To overcome this problem, the GMRES(m) is formulated as a control problem in order to adaptively combine two strategies: a) the appropriate variation of the restarted parameter m, if a stagnation in the convergence is detected; and b) the augmentation of the search subspace using vectors obtained at previous cycles. The proposal is compared with similar iterative methods of the literature based on standard GMRES(m) with fixed parameters. Numerical results for selected matrices suggest that the switching adaptive proposal method could overcome the stagnation observed in standard methods, and even improve the performance in terms of computational time and memory requirements.

Keywords:
iterative method; adaptive GMRES(m); electromagnetic scattering

RESUMO

Este artigo apresenta uma versão adaptativa do GMRES com reinício (GMRES(m)) para a resolução da aproximação por diferenças finitas da equação de Helmholtz. Foi observado que a escolha do parâmetro de reinicialização m afeta fortemente a convergência do GMRES(m). Para contornar este problema, o GMRES(m) é formulado como um problema de controle, que permite combinar adaptativamente duas estratégias: a) a variação apropriada do parâmetro de reinicialização m, se for detectado um estancamento na convergência; e b) o aumento do subespaço de busca, usando vetores obtidos em ciclos anteriores. A proposta é comparada com métodos iterativos semelhantes aos obtidos na literatura. Os resultados numéricos sugerem que o método adaptativo de comutação pode contornar o estancamento observado em métodos conhecidos e até mesmo melhorar o desempenho computacional e os requisitos de memória.

Palavras-chave:
método iterativo; GMRES(m) adaptável; espalhamento eletromagnético

1. INTRODUCTION

This article deals with the numerical resolution of a Helmholtz scattering problem by using an adaptive version of the restarted GMRES. The resolution of the Helmholtz scattering equation using iterative methods is particularly difficult since the problem is ill-posed for a set of frequencies that physically corresponds to the resonance modes of the domain to be solved, the discretization grid has to be refined as a function of the frequency of the operator, and the oscillatory and non-local structure of the solution affects the numerical methods ¹⁰. As a consequence, fast methods (like multigrid) and preconditioners (like incomplete LU) fail to give fast convergence for discretizations of the Helmholtz equation; in fact, the improvement of numerical methods and preconditioners for this kind of problems are an active area of research ^{7), (10), (11}.

The design of a preconditioner is not a trivial task, since it involves the knowledge of the problem, with the possibility of being extremely expensive. In the context of this work, the discretization of the Helmholtz equation yields a linear system of equations that is hard to solve. This occurs because linear systems can be highly indefinite, leading to difficulties when attempting to extract effective preconditioners ⁵. For instance, the standard incomplete LU is ineffective, and in some cases, it does not assure a better rate of convergence ¹⁹ and when it is used with GMRES, the performance deteriorates as the wavenumber becomes larger ¹⁰. For improvement, a static alternative consists in modifying the diagonal of the matrix A by adding purely imaginary values ¹⁹, while a dynamic alternative consists in using a flexible method allowing to change the preconditioner at each step ^{8), (21}. Unfortunately, the aforementioned strategies require the selection of some parameters, which may be hard to tune. Moreover, when solving the Helmholtz equation, it is generally expected that an iterative method with good convergence properties will benefit from a good preconditioner. But a preconditioner can mask the convergence problems of a certain iterative method. Thus to have good matching towards convergence between the iterative method and the preconditioner, the improvement of the iterative method is an important issue in the resolution of the Helmholtz equation.

Generalized Minimal Residual Method (GMRES) is an iterative method frequently selected for its robustness in problems whose discretization results in a large sparse non-Hermitian linear system ²⁰. This method approximates the solution of the linear system $Au=f$ at each iteration, where $A\in {\mathrm{ℂ}}^{n\times n}$ is nonsingular and $u,f\in {\mathrm{ℂ}}^n$. GMRES is a method based on Krylov subspace which obtains at each iteration an approximate solution by minimizing 2-norm of the residual, building an orthogonal basis for the Krylov subspace. To maintain the computational and memory requirement of the orthogonalization process under control, the GMRES is restarted, i.e., the dimension of the Krylov subspace is allowed to grow to a certain maximum m and then, using the obtained approximation of the solution, a new residual is computed and a new Krylov subspace is built. This procedure allows to maintain the dimension of the Krylov subspace at m at the most, and consequently keep the cost under control for the orthogonalization at each cycle of the GMRES(m).

Unfortunately, the convergence to the solution is not guaranteed if the selection of the fixed parameter m is not appropriate, causing slowdown or, even more serious, stagnation in its rate of convergence ²¹. If the rate of convergence presents stagnation, a simple alternative consists of enlarging the maximum allowed dimension m, enlarging the subspace information. However, this strategy does not always ensure faster convergence (9), ²¹. It is necessary to include another kind of information or modify the search subspace. Several adaptive strategies to modify the restart parameter are encountered in the literature (see for instance ¹⁵). In ¹², the use of a stagnation test was proposed; ²⁵ was based on the difference of the Ritz and harmonic Ritz values; ¹ presented a simple strategy based on the angles between consecutive residual vectors for modifying the restart parameter; and ⁶ introduced a proportional-derivative control-inspired strategy for choosing the parameter m adaptively. Other modification strategies are hybrid iterative methods and acceleration techniques ². Our work is framed into the category of augmented methods that is a class of acceleration techniques. Some augmented methods are presented in ^{2), (16), (17), (23}.

In this work, we concentrate our efforts on the improvement of the convergence of the GMRES method itself. To this end, an adaptive control strategy is proposed which includes information from previous cycles for overcoming the stagnation. The proposed strategy is used for the resolution of a linear system of equations obtained from the discretization of an electromagnetic cavity problem ¹⁴. For this problem, the standard GMRES(m) (and modified implementations using fixed parameters ^{2), (16}) have a poor rate of convergence for large wave numbers ⁷. Numerical results for the electromagnetic cavity problem with large wave numbers illustrate the efficiency of the proposed method.

This paper is organized as follows. In §2, the formulation of GMRES(m) is introduced and the Electromagnetic Scattering Problem is characterized. In §3, the strategies for overcoming the stagnation are described. The numerical results (presented in §4 with conclusion in §5), show that the adaptive strategy improves the convergence of GMRES(m).

Throughout this paper, ℂ^n×m is the set of all n × m complex matrices. I _n is the n × n identity matrix, and e _j is its j-th column. Given a matrix M, M ^T denotes its transpose and M* its conjugate transpose or Hermitian transpose. Notation ║ ║ denotes the 2-norm for vectors and the induced norm for matrices. The inner product is denoted as 〈·,·〉.

2 GMRES(M) AND ITS MODIFICATIONS

GMRES(m) approximates the solution to the linear system at the j-th restart cycle using the previous residual, $r_{j-1}=f-A{u}_{j-1}$, for constructing a Krylov subspace of ${\mathcal{H}}_m(A,r_{j-1})=span\left\{r_{j-1},A{r}_{j-1},...,A^{m-1}{r}_{j-1}\right\}$ of dimension m. The j-th approximation is built as

$u_j=u_{j-1}+{\mathcal{H}}_m(A,r_{j-1}),$ (2.1)

where the index m denotes that the restarting parameter was set to the value m. GMRES(m) obtains an approximate solution which minimizes the 2-norm of the residual r _j , i.e.,

$\underset{u_j\in u_{j-1}+{\mathcal{H}}_m\left(a,r_{j-1}\right)}{min}||f-A{u}_j||.$ (2.2)

To solve this problem, the Arnoldi process is normally used for obtaining an orthonormal basis for the Krylov subspace. At the j-th cycle, the first m steps of this procedure can be expressed as:

$A{V}_m=V_{m+1}{\tilde{H}}_m,$ (2.3)

where $V_m\in {\mathrm{ℂ}}^{n\times m}$ and $V_{m+1}:=\left[V_m{v}_{m+1}\right]\in {\mathrm{ℂ}}^{n\times (m+1)}$ have orthonormal columns and ${\tilde{H}}_m\in {\mathrm{ℂ}}^{(m+1)\times m}$ is the upper Hessenberg matrix formed by an upper matrix H _m of dimension m × m and an entry h _m+1,m placed at position $(m+1,m)$. If the Arnoldi process starts with $v_1=\left(\frac{1}{\beta }\right)r_{j-1}$, where $\beta =||r_{j-1}||$, then by construction the columns of V _m are an orthogonal basis of the subspace ${\mathcal{K}}_m(A,r_{j-1})$.

The approximate solution u _j is obtained by $u_j=u_{j-1}+V_m{y}_j$ ²⁰ and the expression (2.2) becomes,

$\underset{y_j}{min}||r_{j-1}-A{V}_m{y}_j||.$ (2.4)

Therefore, the approximate solution u _j minimizes the residual norm. GMRES(m) does not maintain orthogonality between approximation spaces generated at successive restarts ^{2), (20}. As a result, GMRES(m) exhibits a slower rate of convergence than its counterpart the GMRES. The problem arises when the residual norm between consecutive cycles does not have an adequate decrease, i.e., when an abrupt slowering occurs in its rate of convergence. A more extreme situation appears when the rate of convergence exhibits a stagnation, i.e. when the reduction in the residual norms at consecutive cycles is very small, i.e., $||r_j||\approx ||r_{j-1}||$. Therefore, the residual vectors point in nearly the same direction at the end of every restart cycle, i.e., $\angle (r_j,r_{j-1})\approx 0$, meaning that the GMRES(m) did not introduce a large decrease in the residual norm between consecutive cycles ^{2), (6}.

There are several strategies for modifying GMRES(m) to accelerate and avoid stagnation. In this paper, we focus on the so-called augmented methods. The general idea consists in determining a subspace of dimension s, as the direct sum of two spaces of smaller dimension, denoted as $span\left\{{\nu }_1,...,{\nu }_m,{\omega }_1,...,{\omega }_k\right\}$ where the first m vectors are determined using the standard Arnoldi procedure with the current residual normalized, i.e., $r_{j-1}/||r_{j-1}||$; while the remaining vectors, which consist of the augmented part, contain relevant information saved from outer cycles. Two alternative methods are explored. The first one is the GMRES-E(m, d), proposed in ¹⁶. This method computes $\left\{{\omega }_1,...,{\omega }_d\right\}$ as harmonic Ritz vectors associated with the smallest harmonic Ritz value. This strategy seems to be particularly effective when a priori information on the problem confirms the presence of a group of relative small eigenvalues, which occurs in Electromagnetic Scattering Problem ⁷. The second method is the LGMRES(m, l) proposed in ². In this case, at each restart cycle an approximate solution is constructed by using the first m vectors ν ₁ , . . . , ν _m and an additional basis ω ₁ , . . . , ω _l , in which each of the vectors contains error information of the each previously built l subspaces. Next, the above methods are described.

GMRES-E(m, d): Including approximate eigenvectors. The goal consists in the elimination of the components that supposedly slow down convergence ¹⁶. The strategy consists of enriching the subspace by introducing eigenvectors associated to the problematic eigenvalues. In practice, the approximate eigenvectors are the harmonic Ritz vectors associated to the harmonic Ritz values per cycle ^{3), (16), (18}.

At the j-th cycle, the harmonic Ritz value ${\tilde{\lambda }}_k$ with the associated harmonic Ritz vector ${\vartheta }_k=W_s{g}_k$, where $g_k\in {\mathrm{ℂ}}_s$, with respect to the subspace $A{𝒦}_s(A,r_{j-1})$ satisfies the following expression

$(A{\vartheta }_k-{\tilde{\lambda }}_k{\vartheta }_k)\perp A{\mathcal{K}}_s(A,r_0),$ (2.5)

which implies,

$$\begin{gathered} {\left(A{W}_s\right)}^H\left(A{W}_s{g}_k-{\tilde{\lambda }}_k{W}_s{g}_k\right)=0, \\ {W}_s^H{A}^{H}A{W}_s{g}_k={\tilde{\lambda }}_k{W}_s^H{A}^H{W}_s{g}_k. \end{gathered}$$ (2.6)

Using the Arnoldi relation,

$A{W}_s=V_s+{}_1\tilde{H}_s$ (2.7)

where W _s is a n × s matrix, its first m columns are Arnoldi’s vectors and the last d corresponds to the approximate eigenvectors ϑ _k and V _s+1 is the $n\times (s+1)$ orthonormal matrix whose first $m+1$ columns are the Arnoldi vectors and last d columns result from orthogonalizing the d approximation eigenvectors $({\vartheta }_k,k=1,...,d)$ against the previous columns of Arnoldi’s vectors. ${\tilde{H}}_s$ is an $(s+1)\times s$ upper Hessenberg matrix with its elements constructed by the new orthogonalization process. At the j-th cycle, the reduced generalized eigenvalue problem is defined using the expression (2.6) and the new Arnoldi relation (2.7).

${\tilde{H}}_s^{*}{\tilde{H}}_s{g}_k={\tilde{\lambda }}_k{H}_s^{*}{g}_k,$ (2.8)

where ${\tilde{H}}_s^{*}$ is the hermitian of the upper Hessenberg matrix whose dimension is $s\times (s+1)$ and $H_s^{*}$ is the hermitian of the Hessenberg matrix whose dimension is s × s. Usually the value of s is much less than n. These eigenvalues are called harmonic Ritz values and are the roots of the GMRES residual polynomial ²¹. The eigenvectors associated with these harmonic Ritz values are called harmonic Ritz vectors. Some g _k associated with the smallest ${\tilde{\lambda }}_k$ are needed to deflates the smallest eigenvalues and thus improves the convergence ^{16), (17}. Since the harmonic Ritz vectors are useful only at an specific cycle, it needs to be recomputed at each cycle. For this reason, only the residuals are stored for being reused in the next cycle ¹⁶.

At the end of the j-th restart cycle, GMRES-E(m, d) seeks the approximate solution u _j of the form

$u_j=u_{j-1}+W_s{y}_j,$ (2.9)

such that y _j is obtained by solving the following minimization problem

$||r_j||=||f-A{u}_j||=\underset{y_j\in {\mathrm{ℂ}}^{s+1}}{min}||\beta e_1-{\tilde{H}}_s{y}_j||.$ (2.10)

where $\beta =||r_{j-1}||$. The sequence of the residual norm of the GMRES-E(m, d) has the property of being non-increasing but can not guarantee convergence ^{16), (18}.

LGMRES(m, l): Including the error approximations in the search subspace. The motivation of LGMRES(m, l) is based on preventing an alternating behavior observed in the GMRES(m) residual at consecutive cycles which results in deteriorating the convergence ². LGMRES(m, l) includes approximations to the error in the current search subspace. The error approximation at the j-th restart cycle is defined by using the approximate solution at previous cycles as

${\phi }_{j-1}=u_{j-1}-u_{j-2}$ (2.11)

and ${\phi }_j\equiv 0$ for $j<1$. This error approximation vector is used for augmenting the search subspace K _m (A, r _j−1 ) at the next cycle. Note that ${\phi }_j\in {𝒦}_m(A,r_{j-2})$. Therefore, this error approximation φ_j in some sense represents the space K _m (A, r _j−2 ) generated in the previous cycle and subsequently discarded in the restarting procedure. Hence it is a natural choice for enriching the next approximation space 𝒦_m (A, r _j−1 ).

The augmented approximation space $\mathcal{I}={𝒦}_m(A,r_{j-1})\cup span\left\{{\phi }_k\right\},k=\left\{\right(j-l),...,(j-1\left)\right\}$ has dimension $s=m+l$. This method, instead of using 𝒦 _m (A, r _j−1 ) in equation (2.1), uses the subspace ℐ. The matrix V _s+1 is the $n\times (s+1)$ orthonormal matrix whose first $m+1$ columns are the Arnoldi vectors and the last l columns result from orthogonalizing the l error approximation vectors $({\phi }_k,k=(j-l),...,(j-1\left)\right)$ against the previous columns of Arnoldi vectors. W _s is an n × s matrix whose first m columns are the orthogonalized Arnoldi vectors and the last l columns are the l most recent error approximations (typically normalized so that all columns are of unit length). Then the new relationship at the j-th cycle is

$A{W}_s=V_{s+1}{\tilde{H}}_s$ (2.12)

holds for LGMRES(m, l), ${\tilde{H}}_s$ denotes an $(s+1)\times s$ upper Hessenberg matrix. In practice, $l<<m$, following ² $l\le 3$ is a good choice for for LGMRES(m, l). Similar to expressions (2.9) and (2.10), the approximate solution at the j-th cycle is $u_j=u_{j-1}+W_s{y}_j$ with $W_s=[v_1{v}_2...v_m{\phi }_{j-l}...{\phi }_{j-1}]$ and y _j minimize the residual norm ║r _j ║.

Remark. It is important to remark that LGMRES is not helpful when one of the following situations occurs:

• (a) when GMRES(m) skip angles $(\angle (r_j,r_{j-2}\left)\right)$ are not small;

• (b) when GMRES(m) sequential angles $(\angle (r_j,r_{j-1}\left)\right)$ vary greatly from cycle to cycle;

• (c) when GMRES(m) converges in a small number of iterations; or

• (d) when GMRES(m) skip angles and sequential angles are near zero (indicating stagnation).

LGMRES is not typically a substitute for preconditioning and does not help when a problem stagnates for a given restart parameter. Possible improvements to the algorithm include a robust adaptive variant ².

A-LGMRES-E(m, d, l): Including approximate eigenvectors and errors approximations simultaneously with adaptive restart parameter. The proposed method, an adaptive version of GMRES(m) denoted as A-LGMRES-E(m, d, l), is inspired by the augmented method presented in ¹⁸ that combine the GMRES-E(m, d) and LGMRES(m, l) with fixed parameters. In problems with stagnation, according to remark (d) of LGMRES, the error approximation vectors do not help to improve the rate of convergence, i.e., ${\phi }_{j-1}=u_{j-1}-u_{j-2}\approx 0$. Hence, these errors vector are discarded and only the harmonic Ritz vectors are maintained to enrich the search subspace. As can be seen in the numerical results of the LGMRES-E, keeping constant m and enriching the search subspace, it does not necessarily avoid stagnation. To solve this problem a controller to augment the size of the search subspace is proposed for enlarging the Arnoldi basis, since decreasing the restart parameter does not contribute to an improvement in the convergence ⁴. This is done by adding a positive integer value α to the value m. Thus the search subspace is formed by the $m_{j-1}+\alpha$ Arnoldi vectors and d harmonic Ritz vectors, where m _j−1 is the restart parameter of the previous cycle. The dimension of the search subspace at the j-th cycle is defined by the following rule:

$s_j=\left\{\begin{array}{l}{m}_{j-1}+l+d\text{if}||y_{j-1}||\ge \delta \\ m_{j-1}+\alpha +d\text{if}||y_{j-1}||\ge \delta ,\end{array}\right.$ (2.13)

where δ is the stagnation threshold parameter. A slow convergence is considered to occur when $||y_{j-1}||<\delta$. The α and δ values are constant and are selected according to before works and numerical experiments (see Section 4). According to the previous rule, the matrix W _s in the j-th cycle is formed in the following way,

$W_s=\left\{\begin{array}{l}\left[v_1\cdots v_m{\vartheta }_1\cdots {\vartheta }_d{\phi }_{j-l}\cdots {\phi }_{j-1}\right]if||y_{j-1}||\ge \delta \\ \left[v_1\cdots v_m{v}_{m+1}\cdots v_{m+\alpha }{\vartheta }_1\cdots {\vartheta }_d\right]if||y_{j-1}||<\delta .\end{array}\right.$ (2.14)

Finally, the approximate solution at the j-th cycle is,

$u_j=u_{j-1}+W_s{y}_j,$ (2.15)

where y _j is the vector that minimizes the residual norm ║r _j ║ and W _s is a n × s _j matrix containing the Krylov subspace enriched with information from previous cycles. Baker ² and Morgan ¹⁶ suggest between 1 and 3 as the number of error approximations l and the harmonic Ritz vectors d, respectively; since the increase of these values does not significantly improve the decrease in the number of cycles necessary for convergence. In this work, the values $l=1$ and $d=3$ are chosen considering the above suggestions and giving more importance to the harmonic Ritz vectors, since the error approximations allow to accelerate the convergence but do not avoid stagnation. When the matrix is non-Hermitian and non-normal, as in the case of the problem addressed in this work, it is observed that the best selection of the aforementioned parameters is difficult to obtain since an a priori behavior of the iterative method with respect to the parameters is not completely understood.

The pseudocode for the j-th cycle of the proposed method denoted as A-LGMRES-E(m, d, l) is presented in the Algorithm 1.

Algorithm 1
The j-th cycle of A-LGMRES-E(m, l, d)

Require: Given $A,u_{j-1},r_{j-1},m_{j-1},||y_{j-1}||,{\left\{\phi \right\}}_{j-l}^{j-1},{\left\{\vartheta \right\}}_1^d,\delta ,l,d,\alpha ,m_{max}$
1:	if $||y_{j-1}||\ge \delta$ then
2:		$s_j=m_{j-1}+l+d$
3:		Generate Arnoldi basis and matrix ${\tilde{H}}_s$ using ϑ 1 . . . ϑ d and φ j−l . . . φj−1
4:	else
5:		if $m_j\le m_{max}$ then
6:			$s_j=m_{j-1}+\alpha +d,$
7:			Generate Arnoldi basis and matrix ${\tilde{H}}_s$ using ϑ 1 . . . ϑ d
8:		else
9:			$s_j=m_{max}+d,$
10:			Generate Arnoldi basis and matrix ${\tilde{H}}_s$ using ϑ 1 . . . ϑ d
11:		end if
12:	end if
13:	Find $y_j=argmi{n}_{y\in {\mathrm{ℂ}}^s}||\beta e_1-{\tilde{H}}_{s}y||$, compute u j and r j ;
14:	if $||r_j||<$ tolerance then
15:		stop;
16:	else
17:		$j=j+1$
18:		Compute the error approximations vectors, φ j−l , . . ., φj−1 ,
19:		Compute the harmonic Ritz vectors, ϑ 1, . . ., ϑ d ,
20:	end if

3 AN ELECTROMAGNETIC CAVITY PROBLEM

The electromagnetic problem observed at Figure 1 is focused on a 2-D geometry by assuming that the medium is invariant in the z-direction and nonmagnetic with constant magnetic permeability $\mu (x,y)={\mu }_0$. The ground plane (x-axis) and the wall of the cavity are perfect electric conductors, and the interior of the cavity is filled with inhomogeneous material characterized with its relative permittivity ε _r (x, y) ⁷.

Figure 1:
The scattering geometry. ${\epsilon }_r=1en\Omega {\mathrm{\backslash \Omega }}_1$ and ${\epsilon }_r=2in{\Omega }_1$ ^{7), (24}.

For a transverse magnetic (TM) polarization, in which the magnetic field is transverse to the invariant direction and the electric field is $E=(0,0,u(x,y\left)\right)$, the modeling of the cavity problem yields the Helmholtz equation (3.1) together with a Sommerfeld’s radiation condition imposed at infinity (equation (3.3)):

$\Delta u+k_0^2{\epsilon }_{r}u=f,in\Omega =[0,a]\times [b,0],$ (3.1)

$u=0,onS,$ (3.2)

${\partial }_{n}u=\mathcal{T}\left(u\right)+g,on\Gamma$ (3.3)

where k ₀ is the free space wave number, $\Omega =[0,a]\times [-b,0]\in {\mathrm{ℝ}}^2$ is the problem domain, f is the source term and $f=0$ in the free space, S denotes the walls of cavity, 𝒯 is a non-local boundary operator, Γ is the aperture between the cavity and the free space and $g\left(x\right)=-2i\beta e^{i\alpha x}$.

The discretization by finite differences of equation (3.1) gives:

$\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2}+k_0^2{\epsilon }_r(x_i,y_j)u_{i,j}=f(x_i,y_j)$ (3.4)

for $i=1,2,...,M,j=1,2,...,N$. The expression (3.2) is discretized by

$u_{0,j}=u_{M+1,j}=u_{i,0}=0,i=1,2,...,M,j=1,...,N.$ (3.5)

The discrete form of the non-local boundary condition (3.3) is given by

$\frac{u_{i,N+1}-u_{i,N}}{h_y}=\sum _l^M{g}_{il}{u}_{l,N+1}+g\left(x_i\right),i=1,...,M.$ (3.6)

For more details about the system of equations, see ⁷.

The discretization of Helmholtz problem (3.1)-(3.3) by finite differences yields a linear system of equations $Au=f$, in which the coefficient matrix A is non Hermitian and highly indefinite for large values of the wave number k ₀¹¹. In addition, the obtained linear system of equations is ill-conditioned and the growth of the mesh size (M and N) leads to very large matrices, and hence to high computational costs. Usually, direct methods do not perform well, and a general iterative method is required. GMRES is normally used in this context (non-Hermitian and indefinite matrices ²⁰) but due to the high computational and memory requirements the GMRES(m) is normally used. Unfortunately, as mentioned at the introduction section, §1, the possible problems of convergence of the GMRES(m) are exacerbated by the particularities of the Helmholtz problem described above. This is the reason way the A-LGMRES-E(m, d, l) is designed and numerically analyzed for the resolution of this problem.

The mesh size is determined by accuracy requirements on the discretization. The quantities,

$\frac{\lambda }{h_x}=\frac{2\pi }{k_0{h}_x}=\frac{2\pi }{k}(N+1),$ (3.7)

$\frac{\lambda }{h_y}=\frac{2\pi }{k_0{h}_y}=\frac{2\pi }{k}(M+1),$ (3.8)

are the numbers of mesh points per wavelength in both directions. A commonly employed engineering rule [8] states that, for a second order finite difference,

$\frac{\lambda }{h_x}\le 10or,k_0{h}_x\le \pi /5,$ (3.9)

$\frac{\lambda }{h_y}\le 10or,k_0{h}_y\le \pi /5,$ (3.10)

is required to obtain satisfactory results. This rule is used in this work, hence the number of points per wavelength is maintained as a constant, which means that the grid is refined as the wave number is increasing, i.e., when the wave number is increased, the values of M and N is also increased.

4 NUMERICAL EXPERIMENTS

The numerical experiments consider $\Omega =[0,1]\times [-0.25,0]$ filled with the non-homogeneous medium:

${\epsilon }_r=\left\{\begin{array}{l}2in{\Omega }_1=\left(0.2,0.8\right)\times \left(-0.25,-0.20\right),\\ 1in\Omega \backslash {\Omega }_1.\end{array}\right.$ (4.1)

The experiments were run on a computer with Intel Core i7-6700T with 8GB RAM, using the software MATLAB R2016a for Windows 10. The resulting linear system $Au=f$ was solved with the following iterative methods: GMRES(m) ²⁰; GMRES-E(m, d) ¹⁶ using d harmonic Ritz vectors for augmenting the search subspace, whose dimension at each cycle is $s=m+d$; LGMRES(m, l) ² using l errors approximation vectors for augmenting the search subspace, whose dimension at each cycle is $s=m+l$; LGMRES-E(m, l) ¹⁸ using l errors approximation vectors and d harmonic Ritz vectors for augmenting the search subspace, whose dimension at each cycle is $s=m+l+d$. For enriching the subspace in the proposed A-LGMRES-E(m _j , l, d) method, if $||y_j||\ge \delta$ it is used an adaptive restart parameter $m_j=m_{j-1},d$ harmonic Ritz vectors and l error approximation vectors; while if $||y_j||<\delta$ it is used a restart parameter $m_j=m_{j-1}+\alpha$ and d harmonic Ritz vectors only. In the latter case, the dimension at each cycle is according to rule (2.13).

It has also implemented a modified version of standard GMRES (m), that uses the rule (2.13) to modify adaptively the restart parameter and without any enrichment for the search subspace, i.e., $l=d=0$. This method is denoted as GMRES(m _j ) and it is used for comparison purposes. The strategies that modify the restart parameter include a minimum m _min and a maximum m _max for it, i.e., $m_{min}\le m_j\le m_{max}$. The initial restart parameter is denoted as m ₀ and is chosen $m_0=m_{min}$ for these strategies.

In Table 1 is presented the considered wave numbers k ₀ and size of the grid for the matrices tested of the problem introduced in §3. In this table, the size, the number of nonzero matrix elements and the condition number of the matrix A are defined as size(A), nnz(A) and cond(A) respectively; M, N are discretization parameters. The minimum and the maximum eigenvalue (in magnitude) are represented by λ ₁ and λ _n , respectively. The initial solution is $x_0=0$ for all numerical experiments in this section, the stopping criterion on the relative residual norm is $||r_j||/||r_0||<{10}^{-6}$ or a maximum amount of 3000 restart cycles. The parameters considered for each method are summarized in Table 2. The harmonic Ritz vectors are obtained using the eigs MATLAB function with an initial vector of ones instead of using a random vector as it does by default. This allows to keep the same number of restarts cycles when a problem is running several times. The following parameters are those defined by default in the eigs function. Some of these principal parameters are tolerance of 10⁻⁶ and a maximum number of 300 iterations ²². For the cases where the method does not reach convergence before 3000 restart cycles, the method is stopped and the time is denoted as NC.

Table 1:
List of test problems in the cavity problem with different wave numbers and grids.

A	M	N	k 0	size(A)	nnz(A)	cond(A)	λ 1	λ n
cavity01	39	9	2π	390	3258	1.197E+03	23.16	12595.62
cavity02	59	14	2π	885	7583	1.238E+03	67.02	28593.93
cavity03	99	24	4π	2475	21633	12.469E+03	14.10	79792.69
cavity04	99	24	4π	2475	21633	4.531E+03	27.02	79668.00
cavity05	199	49	4π	9950	88258	43.281E+03	14.21	319668.00
cavity06	299	74	4π	22425	199883	1.120E+05	10.11	719668.01
cavity07	199	49	8π	9950	88258	1.498E+05	3.53	319179.12
cavity08	299	74	8π	22425	199883	1.323E+05	8.88	719179.76
cavity09	399	99	8π	39900	356508	1.761E+05	12.03	1279180.09
cavity10	199	49	10π	9950	88258	14.795E+03	58.84	318816.55
cavity11	299	74	10π	22425	199883	33.825E+03	50.84	718817.42
cavity12	399	99	10π	39900	356508	66.058E+03	42.17	1278817.86

Table 2:
Parameters considered for each method.

Methods	Parameters
GMRES(m)	m=30.
LGMRES(m, l)	m=27, l=3.
GMRES-E(m, d)	m=27, d=3.
LGMRES-E(m, l, d)	m=26, l=1, d = 3.
GMRES(m j )	mmin =30, m max =100, δ=0.5, α=4.
A-LGMRES-E(m j , l, d)	mmin =26, m max =100, l=1, d = 3, δ =0.5, α=4.

Experimentally, it is observed that the values ║y _j ║ are between 1E-01 and 1E-07 for the numerical tests of the GMRES(m) (see Table 3). Note that only four problems converged to the prescribed tolerance, which shows that the problem of the cavity is difficult for the GMRES(m) with m fixed and without subspace enrichment. The problems where the GMRES(m) exhibits stagnation (from cavity05 to cavity12) have in general a mean(║y _j ║) lower than the problems where it converges faster (cavity01 and cavity02). This gives an idea of how susceptible is the method to suffer stagnation in case of the value of m remains constant. In this work, the stagnation threshold parameter is considered to be $\delta =0.5$ and it is of the order of max(║y _j ║) for problems without stagnation (see Table 3). The value of δ is chosen heuristically. Smaller values of δ could modify the values of m unnecessarily (recall that the modification of the value of m is given by the rule (2.13) and it is directly affected by the value of δ). With reference to the value of α, which is the increment of the restart parameter m when there is stagnation, previous works have used very small increments such as ${\Delta }_m=2$ ¹³ and ${\Delta }_m=3$ ⁶ when a stagnation is detected. In this work, we considered the same order of increase, but with a value of 4, i.e., according to our incremental rule, we consider $\alpha =4$.

Table 3:
Results for GMRES(m) for each problem in Table 1.

Problems	min(||y j ||)	mean(||y j ||)	max(||y j ||)	Cycles
cavity01	1.45E-06	4.76E-03	1.94E-01	106
cavity02	1.70E-06	4.98E-03	2.37E-01	77
cavity03	8.07E-07	2.95E-04	5.01E-02	2073
cavity04	1.99E-07	3.95E-04	5.67E-02	1280
cavity05	1.73E-04	5.65E-04	4.37E-02	NC
cavity06	1.82E-05	9.52E-05	1.97E-02	NC
cavity07	5.40E-04	1.78E-03	9.47E-02	NC
cavity08	2.31E-04	7.25E-04	2.95E-02	NC
cavity09	6.05E-05	3.19E-04	1.52E-02	NC
cavity10	1.60E-04	1.90E-03	1.13E-01	NC
cavity11	5.49E-04	1.75E-03	4.69E-02	NC
cavity12	3.32E-04	1.13E-03	2.64E-02	NC

Table 3 presents the behavior of the norm of the vector y _j , when GMRES(m) method is used for all the tested problems. It is observed that the norm of the vector y _j (see equation (2.15)) is small (less than 10−2) indicating a slowdown rate of convergence; i.e., $||r_j||\approx ||r_{j-1}||$. Table 4 shows results on problems from Table 1 for the proposed method A-LGMRES-E. It is listed the numbers of cycles required for converging to $||r_j||/||r_0||<{10}^{-6}$, as well as, min(||y _j ||), mean(||y _j || and max(║y _j ║) which are the minimal, mean and maximal values of ║y _j ║, respectively. A-LGMRESE converged to the prescribed tolerance in all cases, with exception of three problems: the cavity08, cavity09 and cavity12. For these problems, the proposed method achieves the lowest relative norm with respect to the other iterative methods tested (see Table 6).

Table 4:
Results for A-LGMRES-E(m _j ,3,1) for each problem in Table 1.

Problems	min(||y j ||)	mean(||y j ||)	max(||y j ||)	Cycles
cavity01	6.23E-04	1.00E+00	6.51E+00	27
cavity02	3.35E-04	8.62E-01	9.15E+00	31
cavity03	7.79E-05	1.80E-01	5.15E+00	156
cavity04	2.97E-05	3.48E-01	4.65E+00	161
cavity05	4.17E-05	1.54E-01	1.47E+01	587
cavity06	2.54E-05	4.99E-02	2.27E+01	1879
cavity07	3.18E-04	1.27E-01	2.00E+01	2559
cavity08	2.15E-03	1.42E-01	2.55E+01	NC
cavity09	6.34E-03	1.44E-01	5.02E-04	NC
cavity10	1.16E-04	6.11E-01	3.71E+01	753
cavity11	8.92E-06	2.27E-01	2.76E+01	2076
cavity12	9.27E-03	2.64E-01	3.31E+01	NC

Table 5:
Metrics for selected matrices and iterative methods: time in seconds and number of cycles required for $||r_j||/||r_0||<{10}^{-6}$ (Cycles).

	GMRES(m)	GMRES(m j )	LGMRES	GMRES-E	LGMRES-E	A-LGMRES-E
Problems	Time	Time	Time	Time	Time	Time
	(Cycles)	(Cycles)	(Cycles)	(Cycles)	(Cycles)	(Cycles)
cavity01	0.68	0.69	0.47	0.91	0.65	1.02
	(106)	(19)	(68)	(74)	(43)	(27)
cavity02	0.90	1.11	1.14	1.59	1.44	1.44
	(77)	(21)	(101)	(101)	(60)	(31)
cavity03	47.24	30.92	34.38	34.97	15.31	31.84
	(2073)	(141)	(1522)	(1255)	(549)	(156)
cavity04	28.28	31.84	18.68	17.70	10.98	28.98
	(1280)	(149)	(833)	(602)	(385)	(161)
cavity05	NC	390.11	NC	NC	NC	413.64
	(3000)	(563)	(3000)	(3000)	(2000)	(587)
cavity06	NC	2865.62	NC	NC	NC	2309.83
	(3000)	(2662)	(3000)	(3000)	(3000)	(1879)
cavity07	NC	1900.67	NC	NC	NC	1821.93
	(3000)	(2662)	(3000)	(3000)	(3000)	(2559)
cavity08	NC	NC	NC	NC	NC	NC
	(3000)	(3000)	(3000)	(3000)	(3000)	(3000)
cavity09	NC	NC	NC	NC	NC	NC
	(3000)	(3000)	(3000)	(3000)	(3000)	(3000)
cavity10	NC	706.96	NC	224.16	225.09	421.43
	(3000)	(1086)	(3000)	(2626)	(2637)	(753)
cavity11	NC	NC	NC	NC	NC	2510.73
	(3000)	(3000)	(3000)	(3000)	(3000)	(2076)
cavity12	NC	NC	NC	NC	NC	NC
	(3000)	(3000)	(3000)	(3000)	(3000)	(3000)

Table 6:
Relative residual norm for problems that do not converge before 3000 cycles for each method of Table 2.

Problems	GMRES(m)	GMRES(m j )	LGMRES	GMRES-E	LGMRES-E	A-LGMRES-E
	$||r_j||/||r_0||$	$||r_j||/||r_0||$	$||r_j||/||r_0||$	$||r_j||/||r_0||$	$||r_j||/||r_0||$	$||r_j||/||r_0||$
cavity08	7.45E-04	3.06E-04	6.82E-04	4.86E-04	2.68E-04	3.00E-06
cavity09	3.49E-03	1.40E-03	3.11E-03	3.27E-03	2.72E-03	9.45E-05
cavity12	1.21E-03	4.66E-04	1.23E-03	6.86E-04	6.52E-04	2.68E-05

Figures 2-(a) and 2-(b) show the convergence curves and the value ||y _j || versus the number of restart cycles for the problems cavity03 and cavity04. It is observed that the A-LGMRES-E method, which includes information from previous cycles and updates the restart parameter according to rule (2.13), presents a similar rate of convergence to the GMRES(m _j ) method, which updates only the restart parameter without including information from previous cycles. Also, in Figure 2, the search subspace dimension of the adaptive methods is compared in the second sub-figure of columns (a) and (b), that is, the value m _j for GMRES(m _j ) and the value s _j for A-LGMRES-E(s _j ). The value s _j has lower growth than m _j , and this allows a smaller number of matrix-vector multiplications for the A-LGMRES-E(s _j ) in the first fifty cycles. In the third sub-figure, the controlled variable ||y _j || are presented for the two methods. The values ||y _j || of A-LGMRES-E(s _j ) are relatively larger than the corresponding values provided by GMRES(m _j ) but lower than the threshold δ in some cases, allowing the increase of the value s _j more slowly than the m _j of GMRES(m _j ).

Figure 2:
Comparative numerical results for adaptive methods. Relative residual norm $\left(||r_j||/||r_0||\right)$; search subspace dimension (m _j for GMRES(m _j ) and s _j for A-LGMRES-E(s _j )) and controlled variable ║y _j ║ for (a) cavity03 and (b) cavity04.

Comparing the augmented methods with information of previous cycles (see Figure 3), i.e. the LGMRES(m, l), GMRES-E(m, d), LGMRES-E(m, l, d) and A-LGMRES-E(m _j , l, d); the A-LGMRES-E has better rate of convergence for the showed problems, having lowest execution times and number of restart cycles. Furthermore, the LGMRES(m, l) does not converge for the problems cavity06 and cavity11 (Figures 3-(c) and 3-(d)). This is also true for the LGMRESE(m, l, d), which combines error approximations and harmonic Ritz vectors ¹⁸, when constant values of m, l and d are used. This shows that the addition of information vectors from previous cycles with fixed restart parameter is not enough to get the convergence, so an adjustment of m is needed to improve the rate of convergence.

Figure 3:
Convergence curves for the augmented methods A-LGMRES-E, LGMRES-E, LGMRES and GMRES-E: relative residual norm vs. number of restart cycles, for (a) cavity03, (b) cavity04, (c) cavity06 and (d) cavity10. Parameters for the adaptive strategy A-LGMRES-E: $\delta =0.5,\alpha =4$. NC means that the method does not reach convergence before 3000 cycles.

All implemented methods are compared in Table 5. It is observed that the method A-LGMRESE(m _j , l, d) has lowest values for the number of cycles necessary for converging and produces steepest decrease in the rate of convergence with respect to the methods that use fixed parameter. For problems that do not converge, the proposed method achieves the best reduction of the relative residual norm ║r _j ║/║r ₀║ with respect to the other methods tested. The best reduction of the relative residual norm for difficult problems are indicated by boldface in Table 6. It is observed that the methods with fixed parameters achieved lower reduction of the relative residual norm in cycles carried out with respect to the other tested methods.

5 CONCLUSION

An adaptive method based on a threshold criterion was introduced for identifying stagnation in the GMRES(m). The criterion yields the expansion of the search subspace for both improving the speed and overcoming the stagnation. The proposed method, as well as several standard methods, were implemented for the resolution of the finite difference approximation of the Helmholtz equation. Numerical experiments for different discrete domain sizes and values of wave number k ₀ were compared, and the results show that the proposal is good enough to improve convergence when comparing with other iterative methods with either fixed parameters and enriched subspace, or adaptive parameters without subspace enrichment. The computation is especially challenging when k ₀ is increased. In this case, a more exhaustive research is necessary for linear systems with large k ₀ in order to identify what is more convenient; either to modify the restart parameter m or an appropriate augmentation of the search subspace for GMRES(m).

ACKNOWLEDGEMENTS

GEE thanks the technical support given by NIDTEC, Polytechnic School, UNA. JCC acknowledges the financial support given by CONACyT through scholarship Prociencia-2015 and the Vinculation Program for Scientists and Technologists PVCT18-32. CES thanks to the PRONII-CONACyT program.

REFERENCES

• ¹
  A. Baker, E. Jessup & T. Kolev. A simple strategy for varying the restart parameter in GMRES(m). Journal of Computational and Applied Mathematics, 230(2) (2009), 751-761. doi:10.1016/J.CAM.2009.01.009. URL https://www.sciencedirect.com/science/article/pii/S0377042709000132
  » https://doi.org/10.1016/J.CAM.2009.01.009» https://www.sciencedirect.com/science/article/pii/S0377042709000132
• ²
  A. Baker, E. Jessup & T. Manteuffel. A Technique for Accelerating the Convergence of Restarted GMRES. SIAM Journal on Matrix Analysis and Applications, 26(4) (2005), 962-984. doi:10.1137/ S0895479803422014.
  » https://doi.org/10.1137/ S0895479803422014
• ³
  J.C. Cabral, P. Torres & C.E. Schaerer. On control strategy for stagnated GMRES(m). J. Comp. Int. Sci., 14 (2017), 109-118. URL http://epacis.net/jcis/PDF_JCIS/JCIS11-art.0127.pdf
  » http://epacis.net/jcis/PDF_JCIS/JCIS11-art.0127.pdf
• ⁴
  J.C. Cabral F. & C.E. Schaerer. On adaptive strategy for overcome stagnation in LGMRES(m,l). In “XXXVIII Iberian Latin-American Congress on Computational Methods in Engineering” (2017).
• ⁵
  E. Chow & Y. Saad. Experimental study of ILU preconditioners for indefinite matrices. Journal of Computational and Applied Mathematics, 86(2) (1997), 387-414.
• ⁶
  R. Cuevas Núñez, C.E. Schaerer & A. Bhaya. A proportional-derivative control strategy for restarting the GMRES(m) algorithm. Journal of Computational and Applied Mathematics, 337 (2018), 209-224. doi:10.1016/J.CAM.2018.01.009. URL https://www.sciencedirect.com/science/article/pii/S037704271830030X
  » https://doi.org/10.1016/J.CAM.2018.01.009» https://www.sciencedirect.com/science/article/pii/S037704271830030X
• ⁷
  K. Du. GMRES with adaptively deflated restarting and its performance on an electromagnetic cavity problem. Applied Numerical Mathematics, 61(9) (2011), 977 - 988. doi: https://doi.org/10.1016/j.apnum.2011.04.003. URL http://www.sciencedirect.com/science/article/pii/S0168927411000705
  » https://doi.org/10.1016/j.apnum.2011.04.003» http://www.sciencedirect.com/science/article/pii/S0168927411000705
• ⁸
  H.C. Elman, O.G. Ernst & D.P. O’leary. A multigrid method enhanced by Krylov subspace iteration for discrete Helmholtz equations. SIAM Journal on scientific computing, 23(4) (2001), 1291-1315.
• ⁹
  M. Embree. The Tortoise and the Hare Restart GMRES. SIAM Review, 45(2) (2003), 259-266. doi:10.1137/S003614450139961. URL http://epubs.siam.org/doi/10.1137/S003614450139961
  » https://doi.org/10.1137/S003614450139961» http://epubs.siam.org/doi/10.1137/S003614450139961
• ¹⁰
  O.G. Ernst & M.J. Gander. Why it is difficult to solve Helmholtz problems with classical iterative methods. In “Numerical analysis of multiscale problems”. Springer (2012), pp. 325-363.
• ¹¹
  M.J. Gander, I.G. Graham & E.A. Spence. Applying GMRES to the Helmholtz equation with shifted Laplacian preconditioning: what is the largest shift for which wavenumber-independent convergence is guaranteed? Numerische Mathematik, 131(3) (2015), 567-614.
• ¹²
  M. Habu & T. Nodera. GMRES(m) algorithm with changing the restart cycle adaptively. In “Proceedings of Algorithmy 2000 Conference on Scientific Computing” (2000), pp. 254-263.
• ¹³
  M. Habu & T. Nodera. GMRES(m) algorithm with changing the restart cycle adaptively. In “Proceedings of Algorithmy 2000 Conference on Scientific Computing” (2000), pp. 254-263.
• ¹⁴
  J. Jin. “The Finite Element Method in Electromagnetics”. Wiley (2015).
• ¹⁵
  W. Joubert. On the convergence behavior of the restarted GMRES algorithm for solving nonsymmetric linear systems. Numerical Linear Algebra with Applications, 1(5) (1994), 427-447. doi:10.1002/nla.1680010502. URL http://doi.wiley.com/10.1002/nla.1680010502
  » https://doi.org/10.1002/nla.1680010502» http://doi.wiley.com/10.1002/nla.1680010502
• ¹⁶
  R. Morgan. A Restarted GMRES Method Augmented with Eigenvectors. SIAM Journal on Matrix Analysis and Applications, 16(4) (1995), 1154-1171. doi:10.1137/S0895479893253975.
  » https://doi.org/10.1137/S0895479893253975
• ¹⁷
  R.B. Morgan. GMRES with Deflated Restarting. SIAM Journal on Scientific Computing, 24(1) (2002), 20-37. doi:10.1137/S1064827599364659. URL http://epubs.siam.org/doi/10.1137/S1064827599364659
  » https://doi.org/10.1137/S1064827599364659» http://epubs.siam.org/doi/10.1137/S1064827599364659
• ¹⁸
  Q. Niu & L. Lu. Restarted GMRES method Augmented with the Combination of Harmonic Ritz Vectors and Error Approximations. International Journal of Mathematical and Computational Sciences, 6 (2012). URL https://pdfs.semanticscholar.org/111c/4cbc91b183467e57b60b5485872646c5d028.pdf
  » https://pdfs.semanticscholar.org/111c/4cbc91b183467e57b60b5485872646c5d028.pdf
• ¹⁹
  D. Osei-Kuffuor & Y. Saad. Preconditioning Helmholtz linear systems. Applied numerical mathematics, 60(4) (2010), 420-431.
• ²⁰
  Y. Saad & M. Schultz. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7(3) (1986), 856-869. doi:10.1137/0907058.
  » https://doi.org/10.1137/0907058
• ²¹
  V. Simoncini & D.B. Szyld. Recent computational developments in Krylov subspace methods for linear systems. Numerical Linear Algebra with Applications, 14(1) (2007), 1-59.
• ²²
  G.W. Stewart. A Krylov-Schur algorithm for large eigenproblems. SIAM Journal on Matrix Analysis and Applications, 23(3) (2002), 601-614.
• ²³
  Z. Teng & X. Wang. Heavy Ball Restarted CMRH Methods for Linear Systems. Mathematical and Computational Applications, 23(1) (2018), 10.
• ²⁴
  Y. Wang, K. Du & W. Sun. Preconditioning iterative algorithm for the electromagnetic scattering from a large cavity. Numerical Linear Algebra with Applications, 16(5) (2009), 345-363. doi:10.1002/nla.615. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nla.615
  » https://doi.org/10.1002/nla.615» https://onlinelibrary.wiley.com/doi/abs/10.1002/nla.615
• ²⁵
  L. Zhang & T. Nodera. A new adaptive restart for GMRES(m) method. ANZIAM Journal, 46 (2005), 409-425. doi:10.21914/anziamj.v46i0.968. URL http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/968
  » https://doi.org/10.21914/anziamj.v46i0.968» http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/968

Publication Dates

• Publication in this collection
  30 Apr 2020
• Date of issue
  Jan-Apr 2020

History

• Received
  10 Dec 2018
• Accepted
  10 Jan 2020