# A Twisted Block Tangential Filtering Decomposition Preconditioner

Authors: Ruirui Wang; Qiang Niu; Linzhang Lu
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

For block-tridiagonal linear system of equations, a variant of tangential filtering preconditioners is proposed in this paper. The new variant is based on a twisted block factorization along with certain filtering property. For practical usage, a class of composite preconditioners
tested, which are constructed by combining the twisted tangential filtering decomposition
preconditioner with the classicalILU(0) preconditioner in a multiplicative way. The performance
of the new preconditioners is compared with other classical preconditioners; the superiority and
the weakness of the preconditioners are pointed out.


---

## Body

## 1. Introduction

In many scientific and engineer applications, for example, simulation of laser propagation in a plasma [1] and study of transport in highly heterogeneous porous media [2], we have to numerically solve certain partial differential equations in 2D or 3D. The discretization of these PDEs by finite difference or finite volume schemes usually leads to large block-tridiagonal linear system:(1.1)Ax=b,
with
(1.2)A=(D1U1L1D2⋱⋱⋱Um-1Lm-1Dm)∈ℛn×n,b∈ℛn.Due to the size of the problem, preconditioned Krylov iterative methods have become one of the most popular choices. It is generally recognized that the efficiency of the linear systems solver heavily depends on the property of the preconditioners. Therefore, the construction of robust and efficient preconditioners has become an interesting research topic. Several incomplete block factorization preconditioners have been proposed by many researchers, see, for example,  Axelsson [3–5] and Meurant [6–8]. Frequency filtering preconditioners [9–16] advocated by G. Wittum and his successor are a special kind of incomplete block factorization preconditioner. This class of preconditioning techniques has been illustrated particularly efficient for linear systems arising from the discretization of partial differential equations with discontinuous coefficients.With the development of techniques of parallel computing, developing high-performance preconditioners that are suitable for parallel computing environment is becoming an important topic. In this paper, we propose a tangential filtering preconditioner constructed by the framework of twisted block factorization. Firstly, the constructed preconditioner has a filtering property. Secondly, the construction and solving procedures of the twisted factorization preconditioner are carried out from two sides, which can be done in parallel. The performance of the newly built preconditioner is compared with the tangential filtering preconditioner proposed in [10]. For practical applications, Achdou and Nataf [10] propose to combine the tangential filtering preconditioner with the ILU(0) preconditioner. In this paper, we also consider to combine the twisted tangential filtering preconditioner with the ILU(0) preconditioner in the following way:
(1.3)MITF=(MILU-1+MTBTD-1-MILU-1AMTBTD-1)-1.
The performance of several different preconditioners is compared on some linear systems generated from the discretization of boundary value problems with discontinuous coefficients. The results show that the twisted block factorization preconditioner and its corresponding preconditioner output other preconditioners on some problems.In Section2, we give a brief introduction of twisted block tangential filtering decomposition and then introduce the twisted block tangential filtering decomposition preconditioner. In Section 3, we analyze the properties of the preconditioner. In Section 4, we give numerical experiments to compare the performances of different types of preconditioners.


## 2. A Twisted Block Tangential Filtering Decomposition

The blockLDU factorization of A is
(2.1)A=[T1L1T2⋱⋱Lm-1Tm][T1-1T2-1⋱Tm-1][T1U1T2⋱⋱Um-1Tm],
where Ti, (i=1,…,m) are square invertible ni×ni matrices. Matrices Li-1 (resp. Ui-1) are ni×ni-1 (resp. ni-1×ni ) matrices. The matrices Ti satisfy the induction formula
(2.2)Ti={D1,i=1;Di-Li-1Ti-1-1Ui-1,2≤i≤m.
The block LDU factorization can be written as
(2.3)A=(L+T)T-1(U+T).
Similar to the block LDU factorization, the block UDL factorization of A has the form
(2.4)A=(U+T̂)T̂-1(L+T̂),
where T̂=Blockdiag(T̂1,T̂2⋯T̂m) with T̂i satisfies
(2.5)T̂i={Di-UiT̂i+1-1Li,1≤i≤m-1,Dm,i=m.
Then the twisted block factorization can be written as
(2.6)A=(E+T̃)T̃-1(F+T̃),
where(2.7)E=(0L10⋱⋱Lj-10Uj⋱⋱0Um-10),F=(0U10⋱⋱Uj-10Lj⋱⋱0Lm-10).
The index j satisfies 1<j<m, the matrix T̃=Blockdiag(T̃1,T̃2···T̃m) with the diagonal block T̃i satisfies the following relationship:
(2.8)T̃i={D1,i=1Di-Li-1T̃i-1-1Ui-1,2≤i≤j-1,Dj-Li-1T̃i-1-1Ui-1-UiT̃i+1-1Li,i=j,Di-UiT̃i+1-1Li,j+1≤i≤m-1Dm,i=m.Different from diagonal block matrixDi of A, the matrices T̃i becomes dense quickly. Therefore, factorization (2.6) cannot be used for large problems in practice. However, the framework can be used to build an incomplete twisted block-factorization preconditioner for A. Precisely, we can replace the blocks T̃i by suitably chosen sparse or block-sparse approximations T̅i, (i=1,2,…,m). Then an incomplete factorization preconditioner M is constructed, which has the following form:
(2.9)M=(E+T̅)T̅-1(F+T̅),
with T̅=Blockdiag(T̅1,T̅2⋯T̅m).From (2.9) it is easy to see that solving linear system Mx=f is equivalent to solving the following two linear systems
(2.10)(E+T̅)y=f,(T̅-1F+I)x=y.
By exploiting the structure, both of the linear systems can be solved by the forward and backward sweeps. Suppose y=(y1T,y2T,…,ymT)T, f=(f1T,f2T,…,fmT)T, and x=(x1T,x2T,…,xmT)T according to the block structure of T̅. Then the process of solving Mx=f can be described in Algorithm 1.Remark 2.1.
Each of the solvers for(E+T̅)y=f and (T̅-1F+I)x=y described in Algorithm 1 involves forward and backward sweeps, and the two sweeps have no relationship with each other, so the forward and backward sweeps can be run in parallel.The procedure of constructing of T̅ is consistent with the idea presented in [10]. Suppose we have approximation βi-1 of T̅i-1-1 which satisfies
(2.11)∥T̅i-1βi-1-I∥≤α<1
then
(2.12)∥(T̅i-1βi-1-I)2∥≤∥(T̅i-1βi-1-I)∥2≤α2,
which implies
(2.13)∥T̅i-1(2βi-1-βi-1T̅i-1βi-1)∥≤α2.
It means that 2βi-1-βi-1T̅i-1βi-1 is a better approximation of T̅i-1-1 than βi-1. According to (2.8), we have the following formula for T̅(2.14)T̅i={D1,i=1,Di-Li-1(2βi-1-βi-1T̅i-1βi-1)Ui-1,2≤i≤j-1,Dj-Li-1(2βi-1-βi-1T̅i-1βi-1)Ui-1-Ui(2βi+1-βi+1T̅i+1βi+1)Li,i=j,Di-Ui(2βi+1-βi+1T̅i+1βi+1)Li,j+1≤i≤m-1,Dm,i=m.
Then the new block factorization preconditioner M based on the twisted factorization can be constructed by choosing βi-1 properly. Following the tangential filtering condition proposed in [10], a diagonal approximation βi-1 can be determined such that
(2.15)(M-A)t=0,
where t is a filtering vector.Lemma 2.2.
If the matricesT̅i-1(1≤i-1≤m) are invertible, then one has
(2.16)M-A=Blockdiag(N1,N2,…,Nm)
with
(2.17)Ni={0,i=1,Li-1(βi-1T̅i-1-I)T̅i-1-1(T̅i-1βi-1-I)Ui-1,2≤i≤j-1,Li-1(βi-1T̅i-1-I)T̅i-1-1(T̅i-1βi-1-I)Ui-1+Ui(βi+1T̅i+1-I)T̅i+1-1(T̅i+1βi+1-I)Li,i=j,Ui(βi+1T̅i+1-I)T̅i+1-1(T̅i+1βi+1-I)Li,j+1≤i≤m-1,0,i=m.Proof.
Consider the matrixM-A and observe that
(2.18)2βi-1-βi-1T̅i-1βi-1-T̅i-1-1=-(βi-1T̅i-1-I)T̅i-1-1(T̅i-1βi-1-I),
thus (2.17) holds.Now we consider how to form a diagonal matrixβi-1. Let t=(t1T,t2T,···,tmT)T be a given vector. If there are no zero entries in the vectors Ui-1ti(2≤i≤j) and Liti(j≤i≤m-1), then it is possible to find diagonal matrices βi-1 such that M produces the same effect with A when operating on the filtering vector t, that is,
(2.19)(M-A)t=0.
From (2.17), we can see that it is sufficient to make
(2.20)(T̅i-1βi-1-I)Ui-1ti=0,2≤i≤j,(2.21)(T̅i+1βi+1-I)Liti=0,j≤i≤m-1.
These requirements can be satisfied by setting βi-1 as follows:
(2.22)βi-1=Diag(T̅i-1-1Ui-1fi./Ui-1ti),2≤i≤j,(2.23)βi+1=Diag(T̅i+1-1Liti./Liti),j≤i≤m-1,
where ./ designs the pointwise vector division, and Diag(v) is the diagonal matrix constructed from the vector v. We refer to the preconditioner constructed by the above procedure as twisted block tangential filtering decomposition preconditioner.


## 3. Analysis of the Twisted Tangential Filtering Preconditioner

In this section, we restrictA to be symmetric positive definite, and use A≻B (A≽B) to denote that  A-B is symmetric positive definite (semidefinite). Consider the twisted tangential filtering preconditioner M formed by (2.9), which ensures the filtering property (2.19). Furthermore, the following lemma holds and it has been established in [10].Lemma 3.1.
IfA≻0, then matrices T̅i≽T̃i, 1≤i≤m. Moreover, M≻0 and M-A≽0 hold. The proof is similar to the proof of Lemma  2.1 of [10], so it is omitted here. From Lemma 3.1, one has the following result. The proof can be found in [17, 18].Theorem 3.2.
Let(3.1)A=M-N,
be the splitting of coefficient matrix A induced by the twisted tangential filtering decomposition preconditioner M, then ρ(M-1N)<1.


## 4. Numerical Expriments

In this section, we present some numerical results to test the performance of preconditioners discussed in this paper. The performance of composite preconditioners is compared withMILU. Two kinds of approaches of constructing the filtering preconditioner M are considered. The combination approach (1.3) is used for all the composite preconditioners.Consider the boundary value problem used in [10]
(4.1)η(x)u+div(a(x)u)-div(κ(x)∇u)=finΩ,u=0on∂ΩD,∂u∂n=0on∂ΩN,
where Ω=[0,1]n (n=2, or 3), ∂ΩN=∂Ω∖∂ΩD. The function η, the vector field a, and the tensor κ are the given coefficients of the partial differential operator. In 2D case, we have ∂ΩD=[0,1]×{0,1}, and in 3D case, we have ∂ΩD=[0,1]×{0,1}×[0,1]. Due to the discontinuous coefficients in the PDE equation and the size of A, an efficient preconditioner plays an important role in solving (4.1) by preconditioned iterative methods.Several types of preconditioners tested in our numerical experiments, we outline the notations as follows,MILU: the ILU(0) preconditioner;MTFFD: the tangential frequency filtering decomposition preconditioner;MTBTD: the twisted block tangential filtering decomposition preconditioner;MITF: the composite preconditioner generated by MILU and MTFFD;MITB: the composite preconditioner generated by MILU and MTBTD.Two filtering vectors are tested, the Ritz vector ofA used in [10], and [1,1,…,1]T which is used as a filtering vector in [19]. The index  j is set to be j=m/2, where α denotes the largest integer not exceeding α.The linear systems are solved byFGMRES [20] method preconditioned by the previously mentioned preconditioners. The algorithm is unrestarted and the maximum Krylov subspace is set to be 200. For comparison reasons, the number of iterations of the ILU(0) preconditioned FGMRES method for constructing the filtering test-vector has been chosen to be 20. The algorithm is stopped whenever the relative norm ∥b-Axk∥/∥b∥ is less than 10-12. The exact solution is generated randomly. In the following tables, iters denotes the number of iterations, error denotes the infinite norm of the difference between the final approximate solution and the exact solution. We use “×" to denote that the method fails to converge within 200 iterations. For preconditioners MILU, MTFFD, and MTBTD, every iteration requires only one preconditioner solve, so the total preconditioner solves are equal to the iteration number. For the composite preconditioners, assuming that the ILU(0) preconditioner has the same cost with the filtering preconditioner, so the costs for composite preconditioner is twice of the iteration. All the experiments are performed in  MATLAB [21]. The codes have not been optimized for the highest efficiency and therefore we do not report the time, but we outline the number of iterations.The considered boundary value problems (4.1) are discretized on a regular Cartesian grid with a cell-centred finite volume scheme. Full up-winding is used for the convective term in the partial differential equation. The following five different cases are considered.Case 1.
The advection-diffusion problem with a rotating velocity in two dimensions.The tensorκ is the identity, and the velocity is a=(2π(x2-0.5),2π(x1-0.5))T. The function η is zero. The uniform grid with n×n nodes (n=50,100,200,300) are tested respectively. Table 1 displays the results obtained by using different preconditioners. The ILU(0) preconditioner needs the most iterations to converges. The preconditioners MTFFD and MTBTD have better performances compared with ILU(0). The composite preconditioners are more efficient and MITB works a little better than MITF. When changing the filtering vector, the iteration numbers have a small change, but the Ritz vector needs additional 20 steps to calculate. Figure 1 depicts the convergence curves of FGMRES method preconditioned by several different preconditioners. The filtering vector is set to be t=(1,1,…,1)T. We can see that the composite preconditioners are efficient, the FGMRES method preconditioned by MTFFD and MTBTD produces nearly the same convergence curve.Table 1
Results for  Case1: the advection-diffusion problem with a rotating velocity in two dimensions. Top resuls are using t=(1,1,…,1)T, bottom resuls are using Ritz vector corresponding to the smallest eigenvalue of A.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50601.31e-10436.50e-12441.10e-11184.18e-11162.94e-111001085.40e-10638.00e-12641.69e-11261.57e-10237.52e-112001861.15e-9901.10e-11912.81e-11375.67e-10323.44e-10300××1101.50e-111112.72e-11456.88e-10398.56e-1050601.31e-10435.50e-12441.26e-11185.22e-11163.43e-111001085.40e-10629.30e-12641.49e-11261.82e-10236.79e-112001861.15e-9901.02e-11913.05e-11375.72e-10331.57e-10300××1101.83e-111113.00e-11456.70e-10397.10e-10Figure 1
Convergence history for FGMRES when1/h=100 for the advection-diffusion problem with a rotating velocity in two dimensions.Case 2.
Nonhomogenous problems with large jumps in the coefficients in two dimensions.The coefficientsη and a are both zero. The tensor κ is isotropic and discontinuous. It jumps from the constant value 103 in the ring 1/22≤|x-c|≤1/2, c=(1/2,1/2)T, to 1 outside. We tested uniform grids with n×n nodes, n=100,200,300,400. Table 2 displays the results obtained by using different preconditioners.Table 2
Results for  Case2: Nonhomogenous problems with large jumps in the coefficients in two dimensions. Top results are using t=(1,1,…,1)T as filtering vector, bottom results are using Ritz vector corresponding to the smallest eigenvalue of A.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50616.62e-11435.67e-12441.07e-11182.27e-11169.90e-111001084.71e-10626.91e-12642.23e-11261.43e-10231.01e-102001871.54e-9899.23e-12922.14e-11373.10e-10331.92e-10300××1091.81e-111123.33e-11456.06e-10407.31e-1050616.62e-11435.44e-12458.40e-12183.14e-11163.29e-111001084.71e-10629.16e-12641.85e-11261.07e-10231.08e-102001871.54e-9891.03e-11921.75e-11372.73e-10331.76e-10300××1091.41e-111123.28e-11455.94e-10404.80e-10Case 3.
Skyscraper problems.The tensorκ is isotropic and discontinuous. The domain contains many zones of high permeability which are isolated from each other. Let [x] denote the integer value of x. In 2D, we have
(4.2)κ(x)={103*([10*x2]+1),if[10*xi]=0mod(2),i=1,2,1,otherwise,
and in 3D
(4.3)κ(x)={103*([10*x2]+1),if[10*xi]=0mod(2),i=1,2,3,1,otherwise.
The coefficients η and a are both zero. Tables 3 and 4 display the results obtained by using different preconditioners for 2D and 3D problems.Table 3
Results for  Case3: skyscraper problems. Top results are in 2D and bottom results are in 3D, t=(1,1,…,1)T.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror501522.01e-7××××164.20e-8153.93e-8100××××××262.63e-7244.86e-7200××××××399.50e-7371.90e-6300××××××472.77e-6442.64e-610194.57e-12162.40e-12171.86e-1287.08e-1371.66e-11151901.07e-8××××139.91e-9258.92e-9201322.57e-8××××115.13e-9206.21e-930××××××146.91e-8138.84e-8Table 4
Results for  Case3: skyscraper problems. Top results are in 2D and bottom results are in 3D, t is set to be the Ritz vector corresponding to the smallest eigenvalue of A.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror501522.01e-7××××168.91e-8151.67e-7100××××××262.29e-7253.13e-7200××××××381.74e-6372.48e-6300××××××473.52e-6443.28e-610194.57e-12164.77e-12172.56e-1281.77e-01287.86e-13151901.07e-8××××142.16e-9241.60e-8201322.57e-8××××116.43e-9203.61e-930××××××145.60e-8142.37e-8Case 4.
Convective skyscraper problems.The same  happens  with the Skyscraper  problems except that the velocity field is changed to bea=(1000,1000,1000)T. The tested results are displayed in Tables 5 and 6.Table 5
Results for  Case4: convective skyscraper problems. Top results are in 2D, bottom results are in 3D, t=(1,1,…,1)T.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50953.84e-91393.43e-9××132.96e-10132.92-101001732.81e-8××××181.26e-9215.08e-9200××××××262.31e-9283.05e-8300××××××273.22e-8348.32e-810132.80e-9121.24e-9121.82e-964.65e-1056.61e-915923.03e-10544.61e-10××308.45e-11176.55e-1120712.60e-10347.25e-101171.67e-997.99e-12189.51e-10301161.59e-91052.15e-9××331.0e-10141.30e-9Table 6
Results for  Case4: convective skyscraper problems. Top results are in 2D, bottom results are in 3D, t is set to be the ritz vector corresponding to the smallest eigenvalue of A.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50953.84e-91393.43e-9××135.13e-10134.99e-101001732.81e-8××××181.67e-9217.82e-9200××××××254.41e-9282.13e-8300××××××281.85e-8345.74e-810132.80e-9122.39e-9123.83e-962.31e-1054.56e-915923.03e-10563.17e-10××294.51e-10174.83e-1120712.60e-10346.41e-101218.30e-1092.78e-11186.38e-10301161.59e-91081.78e-9××333.19e-10153.96e-10Case 5.
Anisotropic layers.The domain is made of 10 anisotropic layers with jumps of up to four orders of magnitude and an anisotropy ratio of up to103 in each layer. For 3D problem, the cube is divided in to 10 layers parallel to z=0, of size 0.1, in which the coefficients are constant. The coefficient κx in the ith layer is given by v(i), the latter being the ith component of the vector v=[α,β,α,β,α,β,γ,α,α], where α=1, β=102 and γ=104. We have κy=10κx and κz=1000κx. The velocity field is zero. Numerical results are shown in Tables 7 and 8.Table 7
Results for  Case5: anisotropic layers. Top results are in 2D, bottom results are in 3D, t=(1,1,…,1)T.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50995.88e-8539.38e-9471.61e-8112.44e-8103.78e-81001906.73e-7766.22e-8721.22e-8176.23e-7167.38e-8200××1103.30e-81033.45e-8292.30e-6271.05e-7300××1367.95e-91272.08e-8401.04e-7371.22e-720271.43e-7242.22e-8222.66e-8101.43e-891.22e-830342.55e-7273.31e-8263.97e-8114.87e-8105.98e-840408.14e-7281.98e-8294.97e-8111.75e-7117.12e-8Table 8
Results for  Case5: anisotropic layers. Top results are in 2D, bottom results are in 3D, t is set to be the Ritz vector corresponding to the smallest eigenvalue of A.
MILUMTFFDMTBTDMITFMITB1/hiterserroriterserroriterserroriterserroriterserror50995.88e-8534.91e-9485.58e-9111.85e-8103.78e-81001906.73e-7771.26e-8721.53e-8182.77e-7167.84e-7200××1101.13e-81033.03e-8293.85e-8278.38e-8300××1361.69e-81261.55e-8401.09e-7377.40e-820271.43e-7251.54e-8221.73e-8104.75e-991.09e-830342.55e-7272.71e-8272.65e-8117.30e-8105.71e-840408.14e-7282.48e-8293.11e-8112.82e-7111.15e-7From the tests results presented in this paper, we can see that the composite preconditioners have better performance than using just a single preconditioner. The FGMRES method preconditioned byMTBTD produces nearly the same results as by preconditioner MTFFD, and also for the composite preconditioners MITB and MITF. However, the preconditioner proposed in this paper has the advantage of parallel computation. For Advection-diffusion and nonhomogeneous problems, there is a little difference between using ones or Ritz vector as the filtering vector. Considering the additional costs for Ritz vector, it is reasonable to use the ones as filtering vector.


## 5. Conclusion

In this paper, we introduce a new variant of tangential filtering decomposite preconditionerMTBTD, which is based on the twisted factorization of the coefficient matrix A. The new one is comparable to the preconditioner MTFFD presented in [10]. Considering the process of preconditioning with MTBTD described in Algorithm 1, the preconditioner MTBTD is superior to MTFFD for its parallel property. And for the same reason, the composite preconditioner MITB surpasses MITF.Algorithm 1:Solving Mx=f.
(1)  Solving(E+T̅)y=fy1=T̅1-1f1,ym=T̅m-1fm.for i=2:j-1yi=T̅i-1(fi-Li-1yi-1),end.for i=m-1:-1:j+1yi=T̅i-1(fi-Uiyi+1),end.yj=T̅j-1(fj-Ujyj+1-Lj-1yj-1).(2)  Solving(T̅-1F+I)x=yxj=yj,for i=j-1:-1:1xi=yi-T̅i-1Uiyi+1,endfor i=j+1:mxi=yi-T̅i-1Li-1yi-1,end.



---

*Source: 282307-2009-06-14.xml*