# Approximate Schur-Block ILU Preconditioners for Regularized Solution of Discrete Ill-Posed Problems

Authors: Suzan C. Buranay; Ovgu C. Iyikal
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

High order iterative methods with a recurrence formula for approximate matrix inversion are proposed such that the matrix multiplications and additions in the calculation of matrix polynomials for the hyperpower methods of orders of convergencep
=
4
k
+
3
, where k
≥
1 is integer, are reduced through factorizations and nested loops in which the iterations are defined using a recurrence formula. Therefore, the computational cost is lowered from κ
=
4
k
+
3 to κ
=
k
+
4 matrix multiplications per step. An algorithm is proposed to obtain regularized solution of ill-posed discrete problems with noisy data by constructing approximate Schur-Block Incomplete LU (Schur-BILU) preconditioner and by preconditioning the one step stationary iterative method. From the proposed methods of approximate matrix inversion, the methods of orders p
=
7,11,15,19 are applied for approximating the Schur complement matrices. This algorithm is applied to solve two problems of Fredholm integral equation of first kind. The first example is the harmonic continuation problem and the second example is Phillip’s problem. Furthermore, experimental study on some nonsymmetric linear systems of coefficient matrices with strong indefinite symmetric components from Harwell-Boeing collection is also given. Numerical analysis for the regularized solutions of the considered problems is given and numerical comparisons with methods from the literature are provided through tables and figures.


---

## Body

## 1. Introduction

The numerical solution of many scientific and engineering problems requires the solution of large linear systems of equations in the form(1)
A
x
=
b
, where x
,
b
∈
R
n and A
∈
R
n
×
n is nonsingular and is usually unsymmetric and unstructured matrix. When A is large and sparse, the solution of (1) is generally obtained through Krylov subspace methods such as simple iteration, generalized minimal residual (GMRES), biconjugate gradient (BICG), Quasi-Minimal Residual (QMR), conjugate gradient squared (CGS), biconjugate gradient stabilized (BICGSTAB), and for Hermitian positive-definite problems conjugate gradient (CG) and for Hermitian indefinite problems minimum residual (MINRES). These methods converge very rapidly when the coefficient matrix A is close to the identity matrix [1]. Unfortunately, in many practical applications A is not close to the identity matrix such as the discretization of the Fredholm integral equation of first kind [2]. The robustness and efficiency of Krylov subspace methods can be improved dramatically by using a suitable preconditioner [1, 3]. Basically there are two classes of preconditioners: implicit and explicit preconditioners [4].In implicit preconditioning methods typically first an approximate factorization ofA that is A
=
M
+
D is computed where D is the defect matrix and M is used by an implicit way, as the incomplete LU (ILU) factorization of A, [4]. Experimental study of ILU preconditioners for indefinite matrices is given in [5] to gain a better practical understanding of ILU preconditioners and help improve their reliability, because besides fatal breakdowns due to zero pivots, the major causes of failure are inaccuracy and instability of the triangular solves. Another example for implicit preconditioning is the algorithm BBI (Band Block Inverse) which was presented to compute the blocks of the exact inverse of block band matrix B within a band consisting of p
1
≥
p blocks to the left and q
1
≥
q blocks to the right of the diagonal, in each block row, and is practically viable only when the bandwidths p, q are small and the entries B
i
,
j are scalars or if block matrices have small order (see Chapter 8 of [4]).The explicit preconditioning methods compute an approximationG on explicit form of the inverse A
-
1 and use it to form a preconditioning matrix on an explicit form. That is, for solving (1), the system G
A
x
=
G
b is considered by iteration when A and G are explicitly available [4]. As an example an approximate inverse preconditioner for Toeplitz systems with multiple right-hand sides is given in [6]. The study about factorized sparse approximate inverse preconditioning (FSAI) for solving linear algebraic systems with symmetric positive-definite coefficient matrices is proposed theoretically in [7], whereas the corresponding iterative construction of preconditioners is given in [8]. “For matrices with irregular structure, it is unclear how to choose the nonzero pattern of the FSAI preconditioner at least in a nearly optimal way and the serial costs of constructing this preconditioner can be very high if its nonzero pattern is sufficiently dense” [9] and the authors proposed two new approaches to rise the efficiency. Most recently, in [10] two classes of iterative methods for matrix inversion are proposed and these methods are used as approximate inverse preconditioners to precondition BICG method for solving linear systems (1) with the same coefficient matrix and multiple right sides.There are also hybrid preconditioning methods that form in a two-stage way, where an incomplete factorization (implicit) method is used for a matrix partitioned in blocks and a direct method (explicit) is used to approximate the inverse of the pivot blocks that occur at each stage of the factorization. As an example, in Chapter 7 of [4] when A is an M matrix the methods are constructed by recursively applying a 2
×
2 block partitioning of matrix A, where elimination is to take place and then computing an approximation of its Schur complement. Using an approximate version of Young’s “Property A” for symmetric positive-definite linear systems with a 2
×
2 block matrix it is proved in [11] that the condition number of the Schur complement is smaller than the condition number obtained by the block-diagonal preconditioning.When the matrixA is large and sparse matrix, a technique based on exploiting the idea of successive independent sets has been used to develop preconditioners that are akin to multilevel preconditioners [12, 13], called ILU with multielimination (ILUM). It is well known that the block ILU preconditioners usually perform better than their point counterparts and the block versions of the ILUM factorization preconditioning technique are proposed in [14]. Robustness of BILUM preconditioner by adding interlevel acceleration in the form of solving interlevel Schur complement matrices approximately and using two implementation strategies to utilize Schur complement technique in multilevel recursive ILU preconditioning techniques (RILUM) is proposed in [15]. For problems arising from practical applications, such as those from computational fluid dynamics, the coefficient matrices are often irregularly structured, ill-conditioned, and of very large size. If the underlying PDEs are discretized by high order finite elements on unstructured domains, the coefficient matrices may have many nonzeros in each row. Furthermore, some large blocks may be ill-conditioned or near singular and the standard techniques used to invert these blocks may produce unstable or inaccurate LU causing less accurate preconditioners. Because of these drawbacks simple strategies used in the standard BILUM technique proposed in [14] may become inefficient [16].In order to determine a meaningful approximation of the solution of (1) when the coefficient matrix A is severely ill-conditioned or near singular one typically replaces the linear system (1) by a nearby system(2)
A
~
y
=
b
~
,that is less sensitive to the perturbations of the right-hand side b and the coefficient matrix A
. This replacement is commonly referred to as regularization. Discrete ill-posed problems where both the coefficient matrix A and the right-hand side b are contaminated by noise appear in a variety of engineering applications [17–19]. Tikhonov regularization method [20–23], Singular Value Decomposition (SVD) method, [24, 25], well-posed stochastic extensions method [26], extrapolation techniques [27], the iterative method yielding best accessible estimate of the solution of the integral equation [28], and a technique for the numerical solution of certain integral equations of first kind [29] are the most popular regularization methods.The motivation of this study is firstly to propose high order iterative methods for approximate matrix inversion of a real nonsingular matrixA such that the matrix multiplications and additions in the calculation of matrix polynomials for the hyperpower methods of orders of convergence p
=
4
k
+
3, where k
≥
1 is integer, are reduced through factorizations and nested loops of which the iterations are defined using a recurrence formula, and secondly to give an algorithm that constructs 2
×
2 block incomplete LU decomposition based on approximate Schur complement for the nonsingular coefficient matrix A
~
∈
R
n
×
n (when n is even) of the algebraic linear system of equations arising from the ill-posed discrete problems with noisy data. In this algorithm the methods of orders p
=
7,11,15,19 are applied for approximating the Schur complement matrices and the obtained preconditioners are used to precondition the one step stationary iterative method. Section 6 of the study is devoted for numerical examples. Algorithm 6, based on Algorithm 1 in [10], and the new proposed algorithm called Algorithm 11 are applied to find regularized solutions of algebraic linear systems obtained by discretization of two problems of Fredholm integral equation of first kind with noisy data of which the first example is the harmonic continuation problem [2, 26], and the second example is Phillip’s problem [29]. For numerical comparisons the same problems are solved by Algorithms 6 and 11 with method of order p
=
7 proposed in [30] and methods of orders p
=
11,15,19, proposed in [31]. Furthermore, experiments of both algorithms are conducted on some nonsymmetric linear systems of coefficient matrices with strong indefinite symmetric components from Harwell-Boeing collection. Analysis and investigations are provided by the obtained numerical results. In last section concluding remarks are given based on theoretical and numerical analysis.


## 2. Hyperpower Iterative Methods for Approximate Matrix Inversion

Exhaustive studies are published investigating iterative methods for computing the approximate inverse of a nonsingular square matrixA (see [10, 32–34] and references therein), generalized inverse A
T
,
S
2 [30, 35], outer inverses [31], and Moore-Penrose inverses [36], based on matrix multiplication and addition. Let I denote the n
×
n unit matrix and A
∈
R
n
×
n be nonsingular matrix. We denote the approximate inverse of A at the m
-
t
h iteration by V
m, and the residual matrix by T
m
=
I
-
A
V
m. The most known iterative methods for obtaining approximate inverse of A are the p
-
t
h order hyperpower method [35] for p
≥
2
(3)
V
m
+
1
=
V
m
∑
j
=
0
p
-
1
I
-
A
V
m
j
=
V
m
∑
j
=
0
p
-
1
T
m
j
,
m
=
0,1
,
…
,which requires p matrix by matrix multiplications for each iteration.WhenA is m
×
n matrix of rank r with real or complex elements, equation (3) may be used to find generalized inverses as Moore-Penrose inverse based on the choise of initial inverse see [35]. One of the recent studies One of the recent studies includes the following factorization of (3) for p
=
7 given in [30] for computing generalized inverse A
T
,
S
(
2
)
(4)
V
m
+
1
=
V
m
I
+
T
m
+
T
m
4
I
+
T
m
+
T
m
2
,
m
=
0,1
,
…
.The method in (4) performs 5 matrix by matrix multiplications. Another study is [31] in which several systematic algorithms for factorizations of the hyperpower iterative family of arbitrary orders are proposed for computing outer inverses. Among these methods, 11
-
t
h order method(5)
V
m
+
1
=
V
m
I
+
T
m
I
+
T
m
+
T
m
2
+
T
m
3
I
+
T
m
3
+
T
m
6
,
m
=
0,1
,
…
,
15
-
t
h order method(6)
V
m
+
1
=
V
m
I
+
T
m
+
T
m
2
I
+
T
m
2
+
T
m
4
I
+
T
m
4
+
T
m
8
,
m
=
0,1
,
…
,and the 19
-
t
h order method(7)
V
m
+
1
=
V
m
I
+
T
m
+
T
m
2
I
+
T
m
2
+
T
m
4
I
+
T
m
6
+
T
m
12
,
m
=
0,1
,
…
,requiring 7,7, and 8 matrix by matrix multiplications, respectively, are stated by the authors. In [10] two classes of iterative methods for matrix inversion are proposed and these methods are used as approximate inverse preconditioners to precondition BICG method for solving linear systems (1) of the same coefficient matrix with multiple right sides. For a given integer parameter k
≥
1; Θ
q
k
(
T
m
), q
k
=
3
∗
2
k
+
1, and Θ
~
q
~
k
(
T
m
), q
~
k
=
5
∗
2
k
-
1, are matrix-valued functions, where T
m
=
I
-
A
V
m and V
m is the approximate inverse of A at the m
-
t
h iteration. Class 1 methods are generated by Θ
q
k
(
T
m
) and have order of convergence p
=
q
k
=
3
∗
2
k
+
1 and Class 2 methods have orders p
=
q
~
k
=
5
∗
2
k
-
1, which are generated by Θ
~
q
~
k
(
T
m
). These classes are given by(8)
C
l
a
s
s
1
:
Θ
q
1
T
m
=
I
+
T
m
I
+
T
m
2
I
+
T
m
3
-
T
m
3
,
k
=
1
,
q
1
=
7
Θ
q
k
T
m
=
Θ
q
k
-
1
T
m
I
+
T
m
3
∗
2
k
-
1
-
T
m
3
∗
2
k
-
1
,
q
k
=
3
∗
2
k
+
1
,
k
>
1
V
m
+
1
=
V
m
Θ
q
k
T
m
,
m
=
0,1
,
…
,
(9)
C
l
a
s
s
2
:
Θ
~
q
~
1
T
m
=
I
+
T
m
I
+
T
m
2
I
+
T
m
5
+
T
m
4
,
k
=
1
,
q
~
1
=
9
Θ
~
q
~
k
T
m
=
Θ
~
q
~
k
-
1
T
m
I
+
T
m
5
∗
2
k
-
1
+
T
m
5
∗
2
k
-
1
-
1
,
q
~
k
=
5
∗
2
k
-
1
,
k
>
1
V
m
+
1
=
V
m
Θ
~
q
~
k
T
m
,
m
=
0,1
,
…
.


## 3. A New Family of Methods with Recurrence Formula

LetI denote the n
×
n unit matrix and A
∈
R
n
×
n be nonsingular matrix. We define(10)
ℵ
T
m
=
T
m
T
m
,
(11)
Ω
T
m
=
T
m
+
ℵ
T
m
,
(12)
Γ
T
m
=
ℵ
T
m
ℵ
T
m
,
(13)
Ψ
T
m
=
ℵ
T
m
+
Γ
T
m
,the matrix-valued functions consisting of matrix multiplications and additions and call Family Generator Function to the matrix-valued function Φ
p
(
T
m
), given as(14)
Φ
p
T
m
=
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
-
1
Γ
j
T
m
,
p
=
4
k
+
3
,
k
≥
1
,where Γ
0
T
m
=
I
. We say that I
+
Γ
T
m is 1 nested loop; I
+
Γ
T
m
I
+
Γ
T
m is 2 nested loop. Using this presentation we express (14) in Horner’s form(15)
Φ
p
T
m
=
Ω
T
m
I
+
Ψ
T
m
I
+
Γ
T
m
I
+
Γ
T
m
I
+
Γ
T
m
…
, with k
-
1 nested loops, which can be given by the following recurrence formula:(16)
P
0
=
I
,
P
j
=
I
+
Γ
T
m
P
j
-
1
,
j
=
1,2
,
…
,
k
-
1
,
Φ
^
p
T
m
=
Φ
p
T
m
=
Ω
T
m
I
+
Ψ
T
m
P
k
-
1
. For an integer k
≥
1 and the corresponding Family Generator Function Φ
^
p, p
=
4
k
+
3, we present a family of iterative methods for approximate matrix inversion of A as(17)
V
m
+
1
=
V
m
I
+
Φ
^
p
T
m
,
m
=
0,1
,
…
.Theorem 1.
LetA
∈
R
n
×
n be nonsingular matrix. The family of methods (17) is a factorization of (3) with k
-
1 nested loops for the orders p
=
4
k
+
3, k
≥
1 integer.Proof.
Whenk
=
1, p
=
7 from (16), Φ
^
7
(
T
m
)
=
Ω
T
m
I
+
Ψ
(
T
m
), m
=
0,1
,
…, formula (17) can be written as(18)
V
m
+
1
=
V
m
I
+
T
m
+
ℵ
T
m
I
+
ℵ
T
m
+
Γ
T
m
(19)
=
V
m
∑
j
=
0
6
T
m
j
for
m
=
0,1
,
…
,which is the p
=
7
-
t
h order method (3). Assume that(20)
V
m
+
1
=
V
m
I
+
Φ
^
4
k
+
3
T
m
,
m
=
0,1
,
…
,
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
-
1
Γ
j
T
m
=
V
m
∑
j
=
0
4
k
+
2
T
m
j
=
V
m
∑
j
=
0
p
-
1
T
m
j
, holds true for k, then for k
+
1 it follows that p
=
4
(
k
+
1
)
+
3 and (21)
V
m
+
1
=
V
m
I
+
Φ
^
4
k
+
1
+
3
T
m
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
Γ
j
T
m
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
-
1
Γ
j
T
m
+
Γ
k
T
m
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
-
1
Γ
j
T
m
+
Ω
T
m
Ψ
T
m
Γ
k
T
m
=
V
m
∑
j
=
0
4
k
+
2
T
m
j
+
T
m
+
ℵ
T
m
ℵ
T
m
+
Γ
T
m
Γ
k
T
m
=
V
m
∑
j
=
0
4
k
+
2
T
m
j
+
T
m
4
k
+
3
+
T
m
4
k
+
4
+
T
m
4
k
+
5
+
T
m
4
k
+
6
=
V
m
∑
j
=
0
4
k
+
6
T
m
j
=
V
m
∑
j
=
0
4
k
+
1
+
2
T
m
j
=
V
m
∑
j
=
0
p
-
1
T
m
j
,and this is the p
-
t
h order method (3) for p
=
4
(
k
+
1
)
+
3
.Theorem 2 (see [34]).
LetA
∈
R
n
×
n be a nonsingular matrix. For a given k
≥
1 integer if the method (17) is used, the necessary and sufficient condition for the convergence of (17) to A
-
1 is that ρ
T
0
<
1 holds, where ρ is spectral radius, T
0
=
I
-
A
V
0, and V
0 is the initial approximation. The corresponding sequence of residuals satisfies T
m
+
1
=
T
m
p, p
=
4
k
+
3.Proof.
Using (21), on the basis of Theorem 4 in Section 7 of Chapter 7 given in [35] and using (17), at m
+
1
-
t
h iteration we have(22)
T
m
+
1
=
I
-
A
V
m
+
1
=
I
-
A
V
m
I
+
Φ
^
p
T
m
,
m
=
0,1
,
…
.
=
I
-
A
V
m
∑
j
=
0
p
-
1
T
m
j
,
=
T
m
p
.Denoting the error by E
m
=
A
-
1
-
V
m and using (21),(23)
A
-
1
-
E
m
+
1
=
V
m
+
1
=
V
m
I
+
Φ
^
p
T
m
,
m
=
0,1
,
…
.
for
17
=
V
m
∑
j
=
0
p
-
1
T
m
j
=
A
-
1
-
E
m
I
+
T
m
+
T
m
2
+
⋯
+
T
m
p
-
1
=
A
-
1
-
E
m
I
+
A
E
m
+
A
E
m
2
+
⋯
+
A
E
m
p
-
1
=
A
-
1
-
E
m
A
E
m
p
-
1
.Now using (22) and (23) we obtain(24)
E
m
+
1
=
E
m
A
E
m
p
-
1
=
A
-
1
T
m
p
=
A
-
1
T
m
-
1
p
2
=
…
=
A
-
1
T
0
p
m
+
1
.If ρ
T
0
<
1 then E
m
+
1 converges to zero matrix O
∈
R
n
×
n as m
→
∞
. That is, the proposed family of methods (17) converge to A
-
1 with p
=
4
k
+
3 order for k
≥
1 integer.Lemma 3.
LetA
∈
R
n
×
n be a nonsingular matrix. For a given k
≥
1 integer if the method from the family (17) is applied to find the approximate inverse of A with an initial choice V
0 satisfying A
V
0
=
V
0
A and ρ
T
0
<
1 where T
0
=
I
-
A
V
0 then A
V
m
=
V
m
A for all m
=
0,1
,
…
.Proof.
Proof follows from induction using the proposed family of methods (17).


## 4. Computational Complexity

Let(25)
x
m
=
V
m
b
,
(26)
r
m
=
b
-
A
x
m
,be the corresponding approximate solution of (1) and the residual error, respectively [10]. Let κ be the number of matrix by matrix multiplications (M
M
s) per iteration of p
-
t
h order hyperpower method given in factorized form. Let ν be the number of matrix by matrix additions (M
a
s) other then addition by identity and γ be the number of matrix additions by identity per step. For obtaining an error r
m
/
b
≤
ε by an iterative method obtained by factorization of hyperpower method (3) using nested loops the authors of [10] showed that total number m
=
ln
⁡
ln
⁡
ε
/
ln
⁡
α
(
1
/
ln
⁡
p
) iterations are required where T
0
=
α
<
1. Therefore, as criteria to compare iterative methods originated from (3) the authors of [10] defined Asymptotic Convergence Factor A
C
F
=
κ
/
ln
⁡
p where this factor occurs in the product of κ with iteration number m as κ
m
=
l
n
ln
⁡
ε
/
ln
⁡
α
(
κ
/
ln
⁡
p
)
. When the two factorizations of p
-
t
h order hyperpower method have the same A
C
F values then it is necessary to use another quantity which measures the efficiency of the hyperpower method both with respect to matrix by matrix multiplications and matrix by matrix additions.Definition 4.
We define the asymptotic convergence values by(27)
A
C
V
=
κ
ln
⁡
p
,
ν
ln
⁡
p
,
γ
ln
⁡
p
,where the components of the A
C
V are denoted by(28)
A
C
V
1
=
A
C
F
=
κ
ln
⁡
p
,
A
C
V
2
=
ν
ln
⁡
p
,
A
C
V
3
=
γ
ln
⁡
p
;see also [34].Theorem 5.
LetA
∈
R
n
×
n be a nonsingular matrix. For a given integer k
≥
1 if the method from the family (17) with convergence order p
=
4
k
+
3 is used to compute the approximate inverse of A with an initial approximation V
0 satisfying ρ
T
0
<
1 where T
0
=
I
-
A
V
0, then κ
=
k
+
4, ν
=
2, and γ
=
k
+
2
.Proof.
Consider the proposed family of methods (17). If k
=
1 then p
=
7 and from (16) gives Φ
^
7
(
T
m
)
=
Ω
(
T
m
)
I
+
Ψ
(
T
m
)
=
T
m
+
ℵ
(
T
m
)
I
+
ℵ
(
T
m
)
+
Γ
T
m
), m
=
0,1
,
…, where T
m
=
I
-
A
V
m
. Thus, iteration V
m
+
1
=
V
m
I
+
Φ
^
7
(
T
m
) requires five M
M
s, two M
a
s, and three additions by identity; therefore, κ
=
4
+
1
=
5, ν
=
2, and γ
=
3
=
1
+
2 hold true for k
=
1 per step. Assume that the proposition is true for k, that is, (29)
V
m
+
1
=
V
m
I
+
Φ
^
4
k
+
3
T
m
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
-
1
Γ
j
T
m
,
m
=
0,1
,
…
,
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
I
+
Γ
T
m
I
+
Γ
T
m
…with k
-
1 nested loops requiring κ
=
k
+
4, ν
=
2
,
γ
=
k
+
2 per step. Then for k
+
1 we have (30)
V
m
+
1
=
V
m
I
+
Φ
^
4
k
+
1
+
3
T
m
,
m
=
0,1
,
…
.
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
∑
j
=
0
k
Γ
j
T
m
=
V
m
I
+
Ω
T
m
I
+
Ψ
T
m
I
+
Γ
T
m
I
+
Γ
T
m
… with k nested loops which require κ
=
k
+
5
=
k
+
1
+
4
,
ν
=
2, and γ
=
k
+
3
=
k
+
1
+
2 and hold true for k
+
1.We presentA
C
V from (27) and the total multiplication count M
C and total addition count A
C per iteration required by the proposed methods (17) for integer k
≥
1 in Table 1 where the cost of addition with I or subtracting from I is taken as n. We will use the notations: M
1
(
k
=
1
) for the proposed 7
t
h order method in (17), M
2 for the 7
t
h order method in (4) [30], M
3
(
k
=
2
) for the 11
t
h order method in (17), M
4 for the 11
t
h order method in (5) [31], M
5
(
k
=
3
) for the 15
t
h order method in (17), M
6 for the 15
t
h order method in (6) [31], M
7
(
k
=
4
) for the 19
t
h order method in (17), and M
8 for the 19
t
h order method in (7) [31]. Table 2 shows the rounded A
C
V in (27) of the proposed methods M
1
,
M
3
,
M
5
,
M
7 for the orders p
=
7,11,15,19, respectively, the methods M
2, M
4
,
M
6
,
M
8, and the methods (8) for k
=
1, (9) for k
=
2 (see also [34]). This table demonstrates that the proposed iterative method (17) of convergence order p
=
11 has lower A
C
F, defined in (28), than the method (5) of the same order given in [31]. Moreover, the proposed methods of orders p
=
11,15,19 possess lower A
C
V
2 from (28) than methods (5)-(7) of the same orders, respectively, given in [31].Table 1
Computational complexity of the proposed methods for real nonsingular matrices of sizen
×
n.
Proposed Methods (17)
Order(
p
)
4
k
+
3
κ
k
+
4
A
C
V
k
+
4
l
n
(
4
k
+
3
)
,
2
l
n
(
4
k
+
3
)
,
k
+
2
l
n
(
4
k
+
3
)
M
C
k
+
4
n
3
A
C
k
+
4
n
2
n
-
1
+
2
n
2
+
k
+
2
nTable 2
ACV (27) for the proposed family of methods (17) and for the methods from the literature of the orders p
=
7,11,15,19.
Orderp, Method
κ
ν
γ
A
C
V
p
=
7
M
1
M
2
E
q
.
8
f
o
r
k
=
1
,
[
10
]
5
2
3
2.5695,1.0278,1.5417
5
2
3
2.5695,1.0278,1.5417
6
1
4
3.083,0.5139,2.0556
p
=
11
M
3
M
4
6
2
4
2.5022,0.8341,1.6681
7
3
4
2.9192,1.2511,1.6681
p
=
15
M
5
M
6
7
2
5
2.5849,0.7385,1.8463
7
3
4
2.5849,1.1078,1.4771
p
=
19
M
7
M
8
E
q
.
(
9
)
for
k
=
2
,
[
10
]
8
2
6
2.7170,0.6792,2.0377
8
3
4
2.7170,1.0189,1.3585
10
2
5
(
3.3962,0.6792,1.6981
)


## 5. Algorithms for Numerical Regularized Solution

Discrete ill-posed problems arise from the discretization of ill-posed problems such as Fredholm integral equation of the first kind with smooth kernel [37]. If the discretization leads to the linear system(31)
A
u
h
=
b
,
A
∈
R
l
×
n
,
l
≥
n
,
a
n
d
b
∈
R
lthen all the singular values of A as well as the SVD components of the solution on the average decay gradually to zero and we say that a discrete Picard condition is satisfied; see [37]. In this cases the solution is very sensitive to perturbations in the data such as measurement or approximation errors and thus regularization methods are essential for computing meaningful approximate solutions. In this study we consider the problem of computing stable approximations of discrete ill-posed problems (31) where the data obey the following assumptions: (a) problem (31) is consistent, i.e., b belongs to the column space of A
,
R
A. (b) Instead of b and A we consider a noisy vector b
~
=
b
+
Δ
b and a noisy matrix A
~
=
A
+
Δ
A with available noise levels δ
b
∗ and δ
A
∗, respectively, such that(32)
b
~
-
b
2
≤
δ
b
∗
,
A
~
-
A
2
≤
δ
A
∗
.Discrete ill-posed problems with data satisfying (32) appear in a number of applications such as inverse scattering [17] and potential theory [18]. If l
=
n and that A
~ is nonsingular the obtained perturbed algebraic linear system(33)
A
+
Δ
A
y
=
b
+
Δ
b
,may be solved using Algorithms 6 and 11 given in the following subsections. When l
>
n denote the set of all least squares solutions of (33) by(34)
S
=
y
∈
R
n
∣
A
~
y
-
b
~
2
=
m
i
n
,then y
∈
S if and only if the orthogonality condition(35)
A
~
T
A
~
y
=
A
~
T
b
~
,is satisfied. On the basis of Theorem 1.1.3 in [38] if R
a
n
k
(
A
~
)
=
n then the unique regularized least squares solution is(36)
y
=
A
~
T
A
~
-
1
A
~
T
b
~
.In this case we may apply Algorithms 6 and 11 for the normalized system (35).


### 5.1. Algorithm6 and Convergence Analysis

LetA
~
∈
R
n
×
n be nonsingular matrix. We denote the error of the approximation of the inverse matrix by E
~
m
=
A
~
-
1
-
V
~
m where V
~
m is the approximate inverse of A
~ obtained at the m
-
t
h iteration using the proposed methods (17) for integers k
≥
1. Thus,(37)
y
~
m
=
V
~
m
b
~
,
(38)
r
~
m
=
b
~
-
A
~
y
~
m
,are the corresponding numerical regularized solution and the regularized residual error respectively; also we have T
~
m
=
I
-
A
~
V
~
m
. In this section and the following sections the norm of a matrix . is considered to be a subordinate matrix norm. For a given integer 1
≤
k
≤
4 and the given predescribed accuracy ε
>
0, Algorithm 6 finds the approximate solution y
~
m
~ given in (37) for the perturbed system (33) using the methods (17) by performing m
~ iterations to reach the accuracy r
~
m
~
∞
/
b
~
∞
≤
ε.Algorithm 6.
This algorithm finds the numerical regularized solutiony
~
m
~ of the perturbed system A
~
y
=
b
~, by using the proposed methods (17) for 1
≤
k
≤
4 and is given on the basis of Algorithm 1 in [10].
Step 1. Choose an initial matrix V
~
0 such that T
~
0
<
1 and an integer 1
≤
k
≤
4
.
Step 2 (let m
=
0). Evaluate y
~
0
=
V
~
0
b
~ and r
~
0
=
b
~
-
A
~
y
~
0 using (37), (38), respectively. Then, calculate r
~
0
∞
/
b
~
∞
.
DoStep 3
-Step 5 until r
~
m
∞
/
b
~
∞
≤
ε.
Step 3. Evaluate T
~
m
=
I
-
A
~
V
~
m
.
Step 4. Apply the iteration V
~
m
+
1
=
V
~
m
I
+
Φ
^
p
(
T
~
m
), for the corresponding method in (17) to find V
~
m
+
1
.
Step 5. m
=
m
+
1, and evaluate y
~
m
=
V
~
m
b
~ and r
~
m
=
b
~
-
A
~
y
~
m using (37), (38), respectively. Then, calculate r
~
m
∞
/
b
~
∞
.
Step 6. If m
~ is the iteration number performed, then y
~
m
~ is the approximate regularized solution satisfying r
~
m
∞
/
b
~
∞
≤
ε
.Remark 7.
When the methods from (17) of orders p
=
4
k
+
3, for k
>
4, are to be applied in Algorithm 6 then it should be noted that as p increases then A
C
F value also increases. Therefore, to lower the total computational cost inStep 4 few iterations (m
^) of one of the proposed methods from (17) of order p
=
4
k
+
3, for 1
≤
k
≤
4, can be used first to find an approximate inverse V
~
m
^
. Next, this approximate inverse may be used as an initial approximate inverse for the considered method of order p
=
4
k
+
3, for k
>
4.Theorem 8.
LetA
~
∈
R
n
×
n be a nonsingular matrix and let Algorithm 6 be applied for the solution of the linear perturbed system (33). If the chosen initial approximation V
~
0 satisfies T
~
0
<
1 then the sequence of the norm of the errors T
~
m
+
1, r
~
m
+
1, and y
-
y
~
m
+
1 converges to zero as m
→
∞ and the following error estimates are obtained at the (
m
+
1
)
-
t
h iteration:(39)
T
~
m
+
1
≤
T
~
0
p
m
+
1
,
(40)
r
~
m
+
1
≤
T
~
0
p
m
+
1
b
~
,
(41)
y
-
y
~
m
+
1
≤
C
o
n
d
A
~
T
~
0
p
m
+
1
y
, where T
~
0
=
I
-
A
~
V
~
0, p
=
4
k
+
3, and C
o
n
d
A
~
=
A
~
-
1
A
~.Proof.
The proof is analogous to the proof of Theorem 3 in [10]. From Theorem 2 using (22) and (24) in a subordinate matrix norm we get(42)
T
~
m
+
1
≤
T
~
0
p
m
+
1
,
(43)
E
~
m
+
1
≤
A
~
-
1
T
~
0
p
m
+
1
,respectively, where E
~
m
=
A
~
-
1
-
V
~
m
. Using the residual error and the approximate solution at the (
m
+
1
)
-
t
h step, as r
~
m
+
1
=
b
~
-
A
~
y
~
m
+
1 and y
~
m
+
1
=
V
~
m
+
1
b
~, we get(44)
r
~
m
+
1
=
b
~
-
A
~
V
~
m
+
1
b
~
=
b
~
-
A
~
A
~
-
1
-
E
~
m
+
1
b
~
=
A
~
E
~
m
+
1
b
~
=
T
~
m
+
1
b
~
,and from (22) we obtain(45)
r
~
m
+
1
=
T
~
0
p
m
+
1
b
~
.Using norm inequalities it follows that r
~
m
+
1
≤
T
~
0
p
m
+
1
b
~ and because T
~
0
<
1 the sequence of the norm of the residual errors T
~
m
+
1 and r
~
m
+
1 converges to zero as m
→
∞
. Furthermore,(46)
y
~
m
+
1
=
V
~
m
+
1
b
~
=
A
~
-
1
-
E
~
m
+
1
b
~
=
y
-
E
~
m
+
1
b
~
,using (24), yields y
-
y
~
m
+
1
=
E
~
m
+
1
b
~
=
A
~
-
1
T
~
0
p
m
+
1
b
~, and hence we get(47)
y
-
y
~
m
+
1
≤
A
~
-
1
T
~
0
p
m
+
1
b
~
,and (41) results by using b
~
≤
A
~
y in (47). The sequence of the norm of the errors y
-
y
~
m
+
1 converges to zero as m
→
∞ since T
~
0
<
1
.Theorem 9.
For the linear systems (1) and (33) where A
∈
R
n
×
n is nonsingular matrix, assume that Δ
A
≤
ϵ
A and Δ
b
≤
ϵ
b and that ϵ
C
o
n
d
(
A
)
<
1. If Algorithm 6 is applied for the solution of the linear perturbed system (33) by performing m iterations to compute an approximate inverse V
~
m of A
~
=
A
+
Δ
A with a chosen initial approximation V
~
0 satisfying T
~
0
<
1, then the following normedwise error bounds are obtained:(48)
x
-
y
~
m
≤
ϵ
1
-
ϵ
C
o
n
d
A
A
-
1
b
+
C
o
n
d
A
x
+
A
~
-
1
T
~
0
p
m
b
~
,
(49)
x
-
y
~
m
x
≤
2
ϵ
C
o
n
d
A
1
-
ϵ
C
o
n
d
A
+
A
~
-
1
A
T
~
0
p
m
1
+
ϵ
.Here b
~
=
b
+
Δ
b, y
~
m
=
V
~
m
b
~ is the approximate solution obtained at the m
-
t
h iteration, x is the exact solution of (1), and p
=
4
k
+
3 is the order of the method (17) used in Algorithm 6.Proof.
Proof is obtained using Theorem8 equation (47) and is analogous to the proof of Theorem 5 in [10].


### 5.2. Algorithm11 and Convergence Analysis

WhenA
~
∈
R
n
×
n is nonsingular matrix and n is even, we consider a 2
×
2 block partitioning of A
~ as(50)
A
~
=
A
~
11
A
~
12
A
~
21
A
~
22
,where each block has size n
1
×
n
1 and n
1
=
n
/
2
. When A
~
11 and S
~ are nonsingular we propose a splitting of the perturbed matrix A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗, where L
m
1
∗ and U
m
1
∗ are 2
×
2 block lower and block upper triangular matrices, respectively, as follows:(51)
L
m
1
∗
=
I
n
1
O
n
1
A
~
21
V
~
1
,
m
1
∗
I
n
1
,
U
m
1
∗
=
A
~
11
A
~
12
O
n
1
S
~
.Here I
n
1 and O
n
1 are of size n
1
×
n
1 and present the identity and zero matrices, respectively. S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12 is approximation of the Schur complement S
=
A
~
22
-
A
~
21
A
~
11
-
1
A
~
12, of A
~
11 [39], and V
~
1
,
m
1
∗ is approximate inverse of A
~
11 obtained by using the proposed methods from (17) of order p
=
4
k
+
3, for 1
≤
k
≤
4 with an initial approximate inverse V
~
1,0 satisfying ρ
T
~
1,0
=
ρ
I
-
A
~
11
V
~
1,0
<
1 for an accuracy of T
~
1
,
m
1
∗
∞
<
η
<
1 by performing m
1
∗ iterations, where T
~
1
,
m
1
∗
=
I
-
A
~
11
V
~
1
,
m
1
∗
.Lemma 10.
LetA
~
∈
R
n
×
n be nonsingular matrix and n be even. Assume that A
~
11
,
S
~
∈
R
n
1
×
n
1 are nonsingular and A
~
11 is a leading block of A
~ as partitioned in (50), where n
1
=
n
/
2, S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12, and V
~
1
,
m
1
∗ is an approximate inverse of A
~
11 obtained by using the proposed method from (17) of order p
=
4
k
+
3, for a given integer 1
≤
k
≤
4 such that ρ
T
~
1,0
=
ρ
I
-
A
~
11
V
~
1,0
<
1 by performing m
1
∗ iterations for an accuracy of T
~
1
,
m
1
∗
∞
<
η
<
1. If A
~
11
V
~
1,0
=
V
~
1,0
A
~
11 then(52)
R
m
1
∗
=
O
n
1
O
n
1
-
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
,and that R
m
1
∗ converges to zero matrix O
∈
R
n
×
n as m
1
∗
→
∞, where A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ and L
m
1
∗,U
m
1
∗ are as given in (51). In addition, this incomplete decomposition of A
~ is a convergent splitting if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1.Proof.
Assume thatA
~
11 and S
~ are nonsingular. Hence L
m
1
∗
U
m
1
∗
-
1 is(53)
L
m
1
∗
U
m
1
∗
-
1
=
A
~
11
-
1
-
A
~
11
-
1
A
~
12
S
~
-
1
O
n
1
S
~
-
1
I
n
1
O
n
1
-
A
~
21
V
~
1
,
m
1
∗
I
n
1
.On the basis of Lemma 3, the residual matrix R
m
1
∗
=
L
m
1
∗
U
m
1
∗
-
A
~ is(54)
R
m
1
∗
=
O
n
1
O
n
1
-
A
~
21
I
-
V
~
1
,
m
1
∗
A
~
11
O
n
1
=
O
n
1
O
n
1
-
A
~
21
T
~
1
,
m
1
∗
O
n
1
.From Theorem 2 we have T
~
1
,
m
1
∗
=
T
~
1
,
m
1
∗
-
1
p
=
T
~
1,0
p
m
1
∗, where p
=
4
k
+
3, k
≥
1, and substituting in (54) gives (52). Since ρ
T
~
1,0
<
1 it follows that R
m
1
∗
→
O
∈
R
n
×
n as m
1
∗
→
∞
. Using (52) and (53) we get(55)
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
.The eigenvalues of L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗ are the roots of the equation d
e
t
λ
I
-
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
0 where det denotes the determinant of a square matrix. On the basis of determinant of 2
×
2 lower block matrices in [40] we get(56)
d
e
t
λ
I
-
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
d
e
t
λ
I
n
1
-
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
d
e
t
λ
I
n
1
=
λ
n
1
d
e
t
λ
I
n
1
-
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
.Hence,(57)
ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
.Therefore, on the basis of convergence of one step stationary iterative method (see Theorem 5.3 Chapter 5 in [4]) if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1 then ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
<
1 and the splitting A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ is a convergent splitting.Algorithm 11.
This algorithm constructs approximate2
×
2 block Schur-BILU decomposition of A
~ using the proposed methods (17) for 1
≤
k
≤
4 to approximate the Schur complement matrices and then finds the regularized solution y
~
l
∗ of the perturbed system (33) by using the constructed Schur-BILU decomposition of A
~ to precondition the one step stationary iterative method.
Let stage numberz
=
1, and the subblock B
=
A
~
11.
Step 1. Let m
=
0 and choose an initial matrix V
~
z
,
0 such that
ρ
T
~
z
,
0
=
ρ
I
-
B
V
~
z
,
0
<
1 and the method from (17) for the corresponding integer 1
≤
k
≤
4
. Find T
~
z
,
0
∞
.
DoSteps 2 and 3 until T
~
z
,
m
∞
<
η
<
1.
Step 2. Apply the iteration V
~
z
,
m
+
1
=
V
~
z
,
m
I
+
Φ
^
p
(
T
~
z
,
m
) for the corresponding method in (17) to find V
~
z
,
m
+
1
.
Step 3 (m
=
m
+
1). Evaluate T
~
z
,
m
=
I
-
B
V
~
z
,
m and T
~
z
,
m
∞
.
Step 4. Let m
1
∗ be the total iteration number performed inSteps 2 and 3 and denote the approximate inverse of A
~
11
=
B at m
1
∗ iterations by V
1
,
m
1
∗. Find S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12
.
(i)
An approximate inverse preconditioner forS
~ can be constructed usingSteps 1
-3 by taking z
=
2 and B
=
S
~ and repeatingSteps 1
-3. Let m
2
∗ be the total iteration number performed, then the obtained approximate inverse preconditioner of S
~ is denoted by V
2
,
m
2
∗.
Step 5. Construct the lower block triangular matrix L
m
1
∗ and the upper block triangular matrix U
m
1
∗ as defined in (51).
Step 6. Take l
=
0, and y
~
0
=
0,0
,
…
,
0
T
∈
R
n, and find r
~
0
=
b
~
-
A
~
y
~
0
.
Step 7. Solve the system L
m
1
∗
U
m
1
∗
d
1
=
r
~
0 by solving the block lower and block upper triangular systems. For the solution of block upper triangular system the obtained subsystems may be solved using Algorithm 6. Also preconditioned BICG method or preconditioned GMRES method may be used to solve these block subsystems of which the approximate inverse preconditioners V
~
1
,
m
1
∗ and V
2
,
m
2
∗ obtained inStep 4 can be used as preconditioners. Next denote the approximate solution by d
~
1 and find d
~
1
∞ and increase the value of l to 1.
DoSteps 8 and 9 until d
~
l
∞
<
ε
<
1.
Step 8. Calculate y
~
l
=
y
~
l
-
1
+
d
~
l and r
~
l
=
b
~
-
A
~
y
~
l
.
Step 9. Solve L
m
1
∗
U
m
1
∗
d
l
+
1
=
r
~
l and denote the obtained approximate solution by d
~
l
+
1 and find d
~
l
+
1
∞ and let l
=
l
+
1
.
Step 10. Ifl
∗ is the iteration number performed, inSteps 8 and 9, then y
~
l
∗ is the approximate solution satisfying y
~
l
∗
-
y
~
l
∗
-
1
∞
<
ε
.Theorem 12.
Lete
l
+
1
=
y
-
y
l
+
1 be the error between the exact solution y of (33) and the solution y
l
+
1 of the one step stationary iterative method where y
l
+
1
=
y
l
+
d
l
+
1
, and d
l
+
1 is the exact solution ofL
m
1
∗
U
m
1
∗
d
l
+
1
=
r
^
l and r
^
l
=
b
~
-
A
~
y
l inSteps 7
-9 of Algorithm 11. If the conditions of Lemma 10 are satisfied, then e
l
+
1 and r
^
l
+
1 converge to zero vector as l
→
∞
. Furthermore,(58)
e
l
+
1
∞
≤
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
l
+
1
e
0
∞
,
(59)
r
^
l
+
1
∞
≤
A
~
∞
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
l
+
1
e
0
∞
,and if(60)
max
⁡
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
∞
,
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
∞
<
1
,then e
l
+
1
∞ and r
^
l
+
1
∞ converge to zero for any initial solution y
0 as l
→
∞ where L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗ is as given in (55).Proof.
On the basis of Section 5.2.1 in [4] the proof is as follows: let e
l
=
y
-
y
l, and r
^
l
=
b
~
-
A
~
y
l
(61)
r
^
l
+
1
=
b
~
-
A
~
y
l
+
1
=
b
~
-
A
~
y
-
e
l
+
1
=
A
~
e
l
+
1
.FromSteps 6
-9 of Algorithm 11 and using (55) it follows that(62)
e
l
+
1
=
y
-
y
l
+
1
=
I
-
L
m
1
∗
U
m
1
∗
-
1
A
~
e
l
=
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
e
l
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
e
l
,and by recursion at the l
+
1
-
t
h iteration we obtain(63)
e
l
+
1
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
l
+
1
e
0
.Furthermore, from (61) and (63) we obtain(64)
r
^
l
+
1
=
A
~
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
l
+
1
e
0
;if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1 then ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
<
1 and the splitting A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ is a convergent splitting of A
~ and we get that e
l
+
1 and r
^
l
+
1 converge to the zero vector as l
→
∞. Using norm inequalities we get inequalities (58) and (59). Furthermore, if condition (60) is satisfied thenL
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
<
1 and e
l
+
1
∞ and r
^
l
+
1
∞ converge to zero for any initial solution y
0 as l
→
∞.Remark 13.
Lety
l
+
1 be as given in Theorem 12 and y
~
l
+
1 be the approximate solution obtained by Algorithm 11; the error e
~
l
+
1
=
y
l
+
1
-
y
~
l
+
1 occurs due to the floating points or due to the preconditioned iterative method used to solve the block subsystems in the backward block substitution inSteps 7
-9 of Algorithm 11.


## 5.1. Algorithm6 and Convergence Analysis

LetA
~
∈
R
n
×
n be nonsingular matrix. We denote the error of the approximation of the inverse matrix by E
~
m
=
A
~
-
1
-
V
~
m where V
~
m is the approximate inverse of A
~ obtained at the m
-
t
h iteration using the proposed methods (17) for integers k
≥
1. Thus,(37)
y
~
m
=
V
~
m
b
~
,
(38)
r
~
m
=
b
~
-
A
~
y
~
m
,are the corresponding numerical regularized solution and the regularized residual error respectively; also we have T
~
m
=
I
-
A
~
V
~
m
. In this section and the following sections the norm of a matrix . is considered to be a subordinate matrix norm. For a given integer 1
≤
k
≤
4 and the given predescribed accuracy ε
>
0, Algorithm 6 finds the approximate solution y
~
m
~ given in (37) for the perturbed system (33) using the methods (17) by performing m
~ iterations to reach the accuracy r
~
m
~
∞
/
b
~
∞
≤
ε.Algorithm 6.
This algorithm finds the numerical regularized solutiony
~
m
~ of the perturbed system A
~
y
=
b
~, by using the proposed methods (17) for 1
≤
k
≤
4 and is given on the basis of Algorithm 1 in [10].
Step 1. Choose an initial matrix V
~
0 such that T
~
0
<
1 and an integer 1
≤
k
≤
4
.
Step 2 (let m
=
0). Evaluate y
~
0
=
V
~
0
b
~ and r
~
0
=
b
~
-
A
~
y
~
0 using (37), (38), respectively. Then, calculate r
~
0
∞
/
b
~
∞
.
DoStep 3
-Step 5 until r
~
m
∞
/
b
~
∞
≤
ε.
Step 3. Evaluate T
~
m
=
I
-
A
~
V
~
m
.
Step 4. Apply the iteration V
~
m
+
1
=
V
~
m
I
+
Φ
^
p
(
T
~
m
), for the corresponding method in (17) to find V
~
m
+
1
.
Step 5. m
=
m
+
1, and evaluate y
~
m
=
V
~
m
b
~ and r
~
m
=
b
~
-
A
~
y
~
m using (37), (38), respectively. Then, calculate r
~
m
∞
/
b
~
∞
.
Step 6. If m
~ is the iteration number performed, then y
~
m
~ is the approximate regularized solution satisfying r
~
m
∞
/
b
~
∞
≤
ε
.Remark 7.
When the methods from (17) of orders p
=
4
k
+
3, for k
>
4, are to be applied in Algorithm 6 then it should be noted that as p increases then A
C
F value also increases. Therefore, to lower the total computational cost inStep 4 few iterations (m
^) of one of the proposed methods from (17) of order p
=
4
k
+
3, for 1
≤
k
≤
4, can be used first to find an approximate inverse V
~
m
^
. Next, this approximate inverse may be used as an initial approximate inverse for the considered method of order p
=
4
k
+
3, for k
>
4.Theorem 8.
LetA
~
∈
R
n
×
n be a nonsingular matrix and let Algorithm 6 be applied for the solution of the linear perturbed system (33). If the chosen initial approximation V
~
0 satisfies T
~
0
<
1 then the sequence of the norm of the errors T
~
m
+
1, r
~
m
+
1, and y
-
y
~
m
+
1 converges to zero as m
→
∞ and the following error estimates are obtained at the (
m
+
1
)
-
t
h iteration:(39)
T
~
m
+
1
≤
T
~
0
p
m
+
1
,
(40)
r
~
m
+
1
≤
T
~
0
p
m
+
1
b
~
,
(41)
y
-
y
~
m
+
1
≤
C
o
n
d
A
~
T
~
0
p
m
+
1
y
, where T
~
0
=
I
-
A
~
V
~
0, p
=
4
k
+
3, and C
o
n
d
A
~
=
A
~
-
1
A
~.Proof.
The proof is analogous to the proof of Theorem 3 in [10]. From Theorem 2 using (22) and (24) in a subordinate matrix norm we get(42)
T
~
m
+
1
≤
T
~
0
p
m
+
1
,
(43)
E
~
m
+
1
≤
A
~
-
1
T
~
0
p
m
+
1
,respectively, where E
~
m
=
A
~
-
1
-
V
~
m
. Using the residual error and the approximate solution at the (
m
+
1
)
-
t
h step, as r
~
m
+
1
=
b
~
-
A
~
y
~
m
+
1 and y
~
m
+
1
=
V
~
m
+
1
b
~, we get(44)
r
~
m
+
1
=
b
~
-
A
~
V
~
m
+
1
b
~
=
b
~
-
A
~
A
~
-
1
-
E
~
m
+
1
b
~
=
A
~
E
~
m
+
1
b
~
=
T
~
m
+
1
b
~
,and from (22) we obtain(45)
r
~
m
+
1
=
T
~
0
p
m
+
1
b
~
.Using norm inequalities it follows that r
~
m
+
1
≤
T
~
0
p
m
+
1
b
~ and because T
~
0
<
1 the sequence of the norm of the residual errors T
~
m
+
1 and r
~
m
+
1 converges to zero as m
→
∞
. Furthermore,(46)
y
~
m
+
1
=
V
~
m
+
1
b
~
=
A
~
-
1
-
E
~
m
+
1
b
~
=
y
-
E
~
m
+
1
b
~
,using (24), yields y
-
y
~
m
+
1
=
E
~
m
+
1
b
~
=
A
~
-
1
T
~
0
p
m
+
1
b
~, and hence we get(47)
y
-
y
~
m
+
1
≤
A
~
-
1
T
~
0
p
m
+
1
b
~
,and (41) results by using b
~
≤
A
~
y in (47). The sequence of the norm of the errors y
-
y
~
m
+
1 converges to zero as m
→
∞ since T
~
0
<
1
.Theorem 9.
For the linear systems (1) and (33) where A
∈
R
n
×
n is nonsingular matrix, assume that Δ
A
≤
ϵ
A and Δ
b
≤
ϵ
b and that ϵ
C
o
n
d
(
A
)
<
1. If Algorithm 6 is applied for the solution of the linear perturbed system (33) by performing m iterations to compute an approximate inverse V
~
m of A
~
=
A
+
Δ
A with a chosen initial approximation V
~
0 satisfying T
~
0
<
1, then the following normedwise error bounds are obtained:(48)
x
-
y
~
m
≤
ϵ
1
-
ϵ
C
o
n
d
A
A
-
1
b
+
C
o
n
d
A
x
+
A
~
-
1
T
~
0
p
m
b
~
,
(49)
x
-
y
~
m
x
≤
2
ϵ
C
o
n
d
A
1
-
ϵ
C
o
n
d
A
+
A
~
-
1
A
T
~
0
p
m
1
+
ϵ
.Here b
~
=
b
+
Δ
b, y
~
m
=
V
~
m
b
~ is the approximate solution obtained at the m
-
t
h iteration, x is the exact solution of (1), and p
=
4
k
+
3 is the order of the method (17) used in Algorithm 6.Proof.
Proof is obtained using Theorem8 equation (47) and is analogous to the proof of Theorem 5 in [10].


## 5.2. Algorithm11 and Convergence Analysis

WhenA
~
∈
R
n
×
n is nonsingular matrix and n is even, we consider a 2
×
2 block partitioning of A
~ as(50)
A
~
=
A
~
11
A
~
12
A
~
21
A
~
22
,where each block has size n
1
×
n
1 and n
1
=
n
/
2
. When A
~
11 and S
~ are nonsingular we propose a splitting of the perturbed matrix A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗, where L
m
1
∗ and U
m
1
∗ are 2
×
2 block lower and block upper triangular matrices, respectively, as follows:(51)
L
m
1
∗
=
I
n
1
O
n
1
A
~
21
V
~
1
,
m
1
∗
I
n
1
,
U
m
1
∗
=
A
~
11
A
~
12
O
n
1
S
~
.Here I
n
1 and O
n
1 are of size n
1
×
n
1 and present the identity and zero matrices, respectively. S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12 is approximation of the Schur complement S
=
A
~
22
-
A
~
21
A
~
11
-
1
A
~
12, of A
~
11 [39], and V
~
1
,
m
1
∗ is approximate inverse of A
~
11 obtained by using the proposed methods from (17) of order p
=
4
k
+
3, for 1
≤
k
≤
4 with an initial approximate inverse V
~
1,0 satisfying ρ
T
~
1,0
=
ρ
I
-
A
~
11
V
~
1,0
<
1 for an accuracy of T
~
1
,
m
1
∗
∞
<
η
<
1 by performing m
1
∗ iterations, where T
~
1
,
m
1
∗
=
I
-
A
~
11
V
~
1
,
m
1
∗
.Lemma 10.
LetA
~
∈
R
n
×
n be nonsingular matrix and n be even. Assume that A
~
11
,
S
~
∈
R
n
1
×
n
1 are nonsingular and A
~
11 is a leading block of A
~ as partitioned in (50), where n
1
=
n
/
2, S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12, and V
~
1
,
m
1
∗ is an approximate inverse of A
~
11 obtained by using the proposed method from (17) of order p
=
4
k
+
3, for a given integer 1
≤
k
≤
4 such that ρ
T
~
1,0
=
ρ
I
-
A
~
11
V
~
1,0
<
1 by performing m
1
∗ iterations for an accuracy of T
~
1
,
m
1
∗
∞
<
η
<
1. If A
~
11
V
~
1,0
=
V
~
1,0
A
~
11 then(52)
R
m
1
∗
=
O
n
1
O
n
1
-
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
,and that R
m
1
∗ converges to zero matrix O
∈
R
n
×
n as m
1
∗
→
∞, where A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ and L
m
1
∗,U
m
1
∗ are as given in (51). In addition, this incomplete decomposition of A
~ is a convergent splitting if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1.Proof.
Assume thatA
~
11 and S
~ are nonsingular. Hence L
m
1
∗
U
m
1
∗
-
1 is(53)
L
m
1
∗
U
m
1
∗
-
1
=
A
~
11
-
1
-
A
~
11
-
1
A
~
12
S
~
-
1
O
n
1
S
~
-
1
I
n
1
O
n
1
-
A
~
21
V
~
1
,
m
1
∗
I
n
1
.On the basis of Lemma 3, the residual matrix R
m
1
∗
=
L
m
1
∗
U
m
1
∗
-
A
~ is(54)
R
m
1
∗
=
O
n
1
O
n
1
-
A
~
21
I
-
V
~
1
,
m
1
∗
A
~
11
O
n
1
=
O
n
1
O
n
1
-
A
~
21
T
~
1
,
m
1
∗
O
n
1
.From Theorem 2 we have T
~
1
,
m
1
∗
=
T
~
1
,
m
1
∗
-
1
p
=
T
~
1,0
p
m
1
∗, where p
=
4
k
+
3, k
≥
1, and substituting in (54) gives (52). Since ρ
T
~
1,0
<
1 it follows that R
m
1
∗
→
O
∈
R
n
×
n as m
1
∗
→
∞
. Using (52) and (53) we get(55)
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
.The eigenvalues of L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗ are the roots of the equation d
e
t
λ
I
-
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
0 where det denotes the determinant of a square matrix. On the basis of determinant of 2
×
2 lower block matrices in [40] we get(56)
d
e
t
λ
I
-
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
d
e
t
λ
I
n
1
-
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
d
e
t
λ
I
n
1
=
λ
n
1
d
e
t
λ
I
n
1
-
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
.Hence,(57)
ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
=
ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
.Therefore, on the basis of convergence of one step stationary iterative method (see Theorem 5.3 Chapter 5 in [4]) if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1 then ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
<
1 and the splitting A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ is a convergent splitting.Algorithm 11.
This algorithm constructs approximate2
×
2 block Schur-BILU decomposition of A
~ using the proposed methods (17) for 1
≤
k
≤
4 to approximate the Schur complement matrices and then finds the regularized solution y
~
l
∗ of the perturbed system (33) by using the constructed Schur-BILU decomposition of A
~ to precondition the one step stationary iterative method.
Let stage numberz
=
1, and the subblock B
=
A
~
11.
Step 1. Let m
=
0 and choose an initial matrix V
~
z
,
0 such that
ρ
T
~
z
,
0
=
ρ
I
-
B
V
~
z
,
0
<
1 and the method from (17) for the corresponding integer 1
≤
k
≤
4
. Find T
~
z
,
0
∞
.
DoSteps 2 and 3 until T
~
z
,
m
∞
<
η
<
1.
Step 2. Apply the iteration V
~
z
,
m
+
1
=
V
~
z
,
m
I
+
Φ
^
p
(
T
~
z
,
m
) for the corresponding method in (17) to find V
~
z
,
m
+
1
.
Step 3 (m
=
m
+
1). Evaluate T
~
z
,
m
=
I
-
B
V
~
z
,
m and T
~
z
,
m
∞
.
Step 4. Let m
1
∗ be the total iteration number performed inSteps 2 and 3 and denote the approximate inverse of A
~
11
=
B at m
1
∗ iterations by V
1
,
m
1
∗. Find S
~
=
A
~
22
-
A
~
21
V
~
1
,
m
1
∗
A
~
12
.
(i)
An approximate inverse preconditioner forS
~ can be constructed usingSteps 1
-3 by taking z
=
2 and B
=
S
~ and repeatingSteps 1
-3. Let m
2
∗ be the total iteration number performed, then the obtained approximate inverse preconditioner of S
~ is denoted by V
2
,
m
2
∗.
Step 5. Construct the lower block triangular matrix L
m
1
∗ and the upper block triangular matrix U
m
1
∗ as defined in (51).
Step 6. Take l
=
0, and y
~
0
=
0,0
,
…
,
0
T
∈
R
n, and find r
~
0
=
b
~
-
A
~
y
~
0
.
Step 7. Solve the system L
m
1
∗
U
m
1
∗
d
1
=
r
~
0 by solving the block lower and block upper triangular systems. For the solution of block upper triangular system the obtained subsystems may be solved using Algorithm 6. Also preconditioned BICG method or preconditioned GMRES method may be used to solve these block subsystems of which the approximate inverse preconditioners V
~
1
,
m
1
∗ and V
2
,
m
2
∗ obtained inStep 4 can be used as preconditioners. Next denote the approximate solution by d
~
1 and find d
~
1
∞ and increase the value of l to 1.
DoSteps 8 and 9 until d
~
l
∞
<
ε
<
1.
Step 8. Calculate y
~
l
=
y
~
l
-
1
+
d
~
l and r
~
l
=
b
~
-
A
~
y
~
l
.
Step 9. Solve L
m
1
∗
U
m
1
∗
d
l
+
1
=
r
~
l and denote the obtained approximate solution by d
~
l
+
1 and find d
~
l
+
1
∞ and let l
=
l
+
1
.
Step 10. Ifl
∗ is the iteration number performed, inSteps 8 and 9, then y
~
l
∗ is the approximate solution satisfying y
~
l
∗
-
y
~
l
∗
-
1
∞
<
ε
.Theorem 12.
Lete
l
+
1
=
y
-
y
l
+
1 be the error between the exact solution y of (33) and the solution y
l
+
1 of the one step stationary iterative method where y
l
+
1
=
y
l
+
d
l
+
1
, and d
l
+
1 is the exact solution ofL
m
1
∗
U
m
1
∗
d
l
+
1
=
r
^
l and r
^
l
=
b
~
-
A
~
y
l inSteps 7
-9 of Algorithm 11. If the conditions of Lemma 10 are satisfied, then e
l
+
1 and r
^
l
+
1 converge to zero vector as l
→
∞
. Furthermore,(58)
e
l
+
1
∞
≤
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
l
+
1
e
0
∞
,
(59)
r
^
l
+
1
∞
≤
A
~
∞
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
l
+
1
e
0
∞
,and if(60)
max
⁡
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
∞
,
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
∞
<
1
,then e
l
+
1
∞ and r
^
l
+
1
∞ converge to zero for any initial solution y
0 as l
→
∞ where L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗ is as given in (55).Proof.
On the basis of Section 5.2.1 in [4] the proof is as follows: let e
l
=
y
-
y
l, and r
^
l
=
b
~
-
A
~
y
l
(61)
r
^
l
+
1
=
b
~
-
A
~
y
l
+
1
=
b
~
-
A
~
y
-
e
l
+
1
=
A
~
e
l
+
1
.FromSteps 6
-9 of Algorithm 11 and using (55) it follows that(62)
e
l
+
1
=
y
-
y
l
+
1
=
I
-
L
m
1
∗
U
m
1
∗
-
1
A
~
e
l
=
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
e
l
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
e
l
,and by recursion at the l
+
1
-
t
h iteration we obtain(63)
e
l
+
1
=
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
l
+
1
e
0
.Furthermore, from (61) and (63) we obtain(64)
r
^
l
+
1
=
A
~
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
-
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
O
n
1
l
+
1
e
0
;if ρ
A
~
11
-
1
A
~
12
S
~
-
1
A
~
21
T
~
1,0
p
m
1
∗
<
1 then ρ
L
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
<
1 and the splitting A
~
=
L
m
1
∗
U
m
1
∗
-
R
m
1
∗ is a convergent splitting of A
~ and we get that e
l
+
1 and r
^
l
+
1 converge to the zero vector as l
→
∞. Using norm inequalities we get inequalities (58) and (59). Furthermore, if condition (60) is satisfied thenL
m
1
∗
U
m
1
∗
-
1
R
m
1
∗
∞
<
1 and e
l
+
1
∞ and r
^
l
+
1
∞ converge to zero for any initial solution y
0 as l
→
∞.Remark 13.
Lety
l
+
1 be as given in Theorem 12 and y
~
l
+
1 be the approximate solution obtained by Algorithm 11; the error e
~
l
+
1
=
y
l
+
1
-
y
~
l
+
1 occurs due to the floating points or due to the preconditioned iterative method used to solve the block subsystems in the backward block substitution inSteps 7
-9 of Algorithm 11.


## 6. Numerical Results

The experimental investigation of the proposed algorithms is given on two examples of Fredholm integral equation of first kind and on some nonsymmetric linear systems with strong indefinite symmetric components sourced from simulation of computer systems and nuclear reactor core model. All the computations are performed using a personal computer with properties AMD Ryzen 7 1800X Eight Core Processor 3.60GHz. Calculations are carried by Fortran programs in double precision. In this section the figures and the tables adopt the following notations:T
C
S
A
L
1 is the total solution cost in seconds of Algorithm 6 forSteps 1
-6.
T
C
S
S
P
A
L
1 is the total solution cost in seconds of Algorithm 6 for successive perturbations.
P
C
S
A
L
2 is the cost in seconds for constructing the preconditioner L
m
1
∗
U
m
1
∗ inSteps 1
-5 of Algorithm 11.
I
C
S
A
L
2 is the cost in seconds of the total iterations performed by the preconditioned stationary one step iterative method inSteps 6-9 of Algorithm 11.
T
C
S
A
L
2 is the total solution cost in seconds of Algorithm 11, that is, T
C
S
A
L
2
=
P
C
S
A
L
2+I
C
S
A
L
2.
TBMMs denote the total block matrix by matrix multiplications.
TMMs denote the total matrix by matrix multiplications.Application 14 (harmonic continuation problem).
The first application is the harmonic continuation problem [2, 26]. Given a harmonic function u
(
r
,
θ
) in the unit circle with known values for some r
<
1, u
(
r
,
θ
)
=
f
(
θ
), find its values f
(
θ
) for r
=
1
. Now f
(
θ
) and h
(
θ
) are related by the Poisson integral:(65)
1
2
π
∫
0
2
π
1
-
r
2
1
-
2
r
c
o
s
θ
-
φ
+
r
2
h
φ
d
φ
=
f
θ
.We use the same data used by Franklin [26] by which the answer is known from the real part of the analytic function z
3
-
z
+
sin
⁡
z
. For z
=
1,(66)
h
φ
=
cos
⁡
3
φ
-
cos
⁡
φ
+
sin
⁡
cos
⁡
φ
cosh
⁡
sin
⁡
φ
,whereas for z
=
r
=
0.5,(67)
f
θ
=
1
8
cos
⁡
3
θ
-
1
2
cos
⁡
θ
+
s
i
n
1
2
cos
⁡
θ
cosh
⁡
1
2
sin
⁡
θ
.Using the quadrature nodes φ
j
=
2
π
j
/
n
,
j
=
1
,
…
,
n, evaluating f
(
θ
) at the points θ
i
=
2
π
i
/
n, i
=
1
,
…
,
n, and taking r
=
0.5, the discretization of (65) gives the algebraic linear system A
u
h
=
b where(68)
A
i
,
j
=
3
n
5
-
4
cos
⁡
2
π
i
-
j
/
n
,
b
j
=
f
2
π
j
n
.In all applications given in this section we take Δ
A
=
δ
A
I and Δ
b
=
δ
b
,
δ
b
,
…
,
δ
b
T
∈
R
n satisfying (32) where I is n
×
n identity matrix. The trace of the exact solution h
φ at the grid points φ
j
=
2
π
j
/
n
,
j
=
1
,
…
,
n is denoted by u. The condition number of matrix A
~
=
A
+
Δ
A is denoted by C
o
n
d
(
A
~
) and for Application 14,  C
o
n
d
(
A
~
) with respect to δ
A and the values of δ
b are given in Table 3 where δ
A
=
0.5
δ
b
1.5
. The shape of the coefficient matrix A
~ when δ
A
=
0.5
×
5
×
1
0
-
10
1.5 and n
=
800 of Application 14 is given in Figure 1, illustrating that the coefficient matrix is very dense. The T
C
S
A
L
1 to achieve the desired accuracy r
~
m
∞
/
b
~
∞
≤
5
×
1
0
-
11 when n
=
800 and V
~
0
=
A
~
T
/
A
~
1
A
~
∞ for Application 14 with respect to δ
b are given in Figure 2.
This figure shows that, for small values ofδ
b, the proposed method M
3 (17) requires the minimum T
C
S
A
L
1. Table 4 presents the iteration number m
~ performed, TMMs, T
C
S
A
L
1 to achieve the desired accuracy r
~
m
∞
/
b
~
∞
≤
5
×
1
0
-
11, and the relative L
2 norm of the errors using the methods M
i, i
=
1
,
…
,
8, for Application 14 when n
=
800 and δ
b
=
1
0
-
5
. Next we take δ
b
=
1
0
-
11 and use the methods M
i, i
=
1
,
…
,
8 in Step 4 of Algorithm 6 by applying successive perturbations for solving the perturbed systems (
A
+
δ
A
j
I
)
y
j
=
b
+
Δ
b
j, Δ
b
j
=
[
δ
b
j
,
…
,
δ
b
j
]
T for j
=
1,2
,
3, where δ
b
1
=
δ
b
=
1
0
-
11, δ
A
1
=
0.5
δ
b
1
1.5, and δ
b
j
=
0.999
δ
b
j
-
1, δ
A
j
=
0.999
δ
A
j
-
1
,
j
=
2,3 such that the obtained approximate inverse is used as an initial approximate inverse for the next perturbed system. Let y
~
m
~
j
j
j
=
1,2
,
3, be the approximate solution of y
j, obtained by Algorithm 6 by performing m
~
j iterations for an accuracy of r
~
m
~
j
∞
/
b
~
∞
≤
5
×
1
0
-
11 for the corresponding perturbed system. Table 5 shows the total iteration number m
~
=
m
~
1
+
m
~
2
+
m
~
3 performed, TMMs, T
C
S
S
P
A
L
1 for an accuracy of r
~
m
~
j
∞
/
b
~
∞
≤
5
×
1
0
-
11, and the relative L
2 norm of the errors using the methods M
i, i
=
1
,
…
,
8, by Algorithm 6 for Application 14 when n
=
800. Thus, both Tables 4 and 5 present that the proposed methods M
3, M
5, and M
7 with Algorithm 6 give solutions by performing less total solution cost in time compared with the other methods of same orders, M
4, M
6, and M
8 respectively. The T
C
S
A
L
2 with respect to δ
b by using the methods M
i
,
i
=
1
,
…
,
8, inStep 2 of Algorithm 11 are given in Figure 3 for Application 14 when n
=
800. This figure shows that, for the considered values of δ
b, the proposed method M
3 (17) requires the minimum total solution cost to achieve the desired accuracy T
~
1
,
m
∞
<
0.05 and d
~
l
∞
<
5
×
1
0
-
6 when n
=
800. P
C
S
A
L
2, I
C
S
A
L
2, and T
C
S
A
L
2 to achieve the desired accuracy T
~
1
,
m
∞
<
0.05 and d
~
l
∞
<
5
×
1
0
-
6, iteration numbers m
1
∗, l
∗ and TBMMs, maximum norm, and relative L
2 norm of the errors with respect to δ
b
=
1
0
-
5 of the methods M
i, i
=
1,2
,
…
,
8, by Algorithm 11, when n
=
800 for Application 14, are presented in Table 6.
Moreover, in [2] the best error of numerical solution occurring by singular value decomposition when n
=
50 for the harmonic continuation problem (Application 14) was obtained approximately 1
0
-
3 using single precision. However, relative maximum errors using the proposed methods M
1
,
M
3
,
M
5
,
M
7, for the corresponding discretization by Algorithm 6 when δ
b
=
1
0
-
11 are 2.05
×
1
0
-
9 and by Algorithm 11 when δ
b
=
1
0
-
6 are 1.19
×
1
0
-
6. Furthermore, we solved harmonic continuation problem by direct LU with pivoting and the accuracy of the solution is ≃
1
0
-
5 for n
=
800, and δ
b
=
1
0
-
11
. Also when this problem is solved for same value of n and δ
b by Algorithm 6 using the proposed methods M
1
,
M
3
,
M
5
,
M
7, relative L
2 norm of the errors is less than 5.5
×
1
0
-
11 for the third successive perturbation as shown in column seven of Table 5.Table 3
Condition number of matrixA
~ of Application 14 with respect to δ
A
=
0.5
(
δ
b
)
1.5 when n
=
800.
δ
b
C
o
n
d
(
A
~
)
1
0
-
1
1413.69
1
0
-
2
47633.8
1
0
-
3
1.56798
E
+
06
1
0
-
4
5.09537
E
+
07
1
0
-
5
1.64380
E
+
09
1
0
-
6
5.27663
E
+
10
1
0
-
7
1.68754
E
+
12
1
0
-
8
5.377040
E
+
13
1
0
-
9
1.759400
E
+
15
1
0
-
10
2.503850
E
+
17Table 4
T
C
S
A
L
1, iteration numbers and TMMs, and relative L
2 norm of the errors obtained by Algorithm 6 using the methods M
i, i
=
1
,
.
.
.
,
8, for Application 14, when n
=
800 and δ
b
=
1
0
-
5.
Method
T
C
S
A
L
1
m
~
TMMs
u
-
y
~
m
~
2
u
2
M
1
104.77
8
40
1.6969719
E
-
05
M
2
104.77
8
40
1.6969719
E
-
05
M
3
95.72
6
36
1.6969729
E
-
05
M
4
117.84
6
42
1.6969729
E
-
05
M
5
110.26
6
42
1.6969718
E
-
05
M
6
110.81
6
42
1.6969718
E
-
05
M
7
105.32
5
40
1.6969726
E
-
05
M
8
105.53
5
40
1.6969726
E
-
05Table 5
T
C
S
S
P
A
L
1, iteration numbers, and relative L
2 norm of the errors obtained by Algorithm 6 using the methods M
i, i
=
1
,
.
.
.
,
8, for the perturbed systems of Application 14, when n
=
800 and δ
b
1
=
1
0
-
11.
Method
T
C
S
S
P
A
L
1
m
~
TMMs
u
-
y
~
m
~
1
1
2
u
2
u
-
y
~
m
~
2
2
2
u
2
u
-
y
~
m
~
3
3
2
u
2
M
1
137.06
10
50
7.6070
E
-
09
1.0704
E
-
10
1.7237
E
-
11
M
2
137.06
10
50
7.6070
E
-
09
1.0704
E
-
10
1.7339
E
-
11
M
3
146.40
9
54
3.2302
E
-
10
1.8604
E
-
11
1.8129
E
-
11
M
4
156.88
9
63
3.2298
E
-
10
1.8711
E
-
11
1.7950
E
-
11
M
5
153.44
8
56
1.9954
E
-
09
2.2614
E
-
11
1.8071
E
-
11
M
6
153.51
8
56
1.9954
E
-
09
2.2705
E
-
11
1.7496
E
-
11
M
7
159.29
8
64
9.7145
E
-
11
1.7114
E
-
11
5.4697
E
-
11
M
8
159.86
8
64
9.7083
E
-
11
1.7198
E
-
11
4.1126
E
-
11Table 6
Computational costs, total iteration numbers and TBMMs, maximum norm, and relativeL
2 norm of the errors obtained by Algorithm 11 using the methods M
i, i
=
1
,
.
.
.
,
8, for the perturbed system of Application 14, when n
=
800 and δ
b
=
1
0
-
5.
Method
P
C
S
A
L
2
I
C
S
A
L
2
T
C
S
A
L
2
m
1
∗
TBMMs
l
∗
u
-
y
~
l
∗
∞
u
-
y
~
l
∗
2
u
2
M
1
42.74
0.05
42.79
22
110
4
1.2213
E
-
05
1.6992
E
-
05
M
2
42.74
0.05
42.79
22
110
4
1.2190
E
-
05
1.6992
E
-
05
M
3
40.26
0.05
40.31
18
108
4
1.2270
E
-
05
1.6992
E
-
05
M
4
48.33
0.06
48.39
18
126
4
1.2151
E
-
05
1.6992
E
-
05
M
5
43.45
0.06
43.51
16
112
4
1.2228
E
-
05
1.6993
E
-
05
M
6
43.51
0.06
43.57
16
112
4
1.2215
E
-
05
1.6993
E
-
05
M
7
43.62
0.08
43.70
15
120
5
1.2273
E
-
05
1.6991
E
-
05
M
8
44.42
0.11
44.53
15
120
7
1.2185
E
-
05
1.6991
E
-
05Figure 1
The shape of the coefficient matrixA
~ when δ
A
=
0.5
×
5
×
1
0
-
10
1.5 and n
=
800 for Application 14.Figure 2
Total computational cost comparisons of methodsM
i, i
=
1
,
…
,
8, with respect to δ
b for Application 14 when n
=
800 by Algorithm 6.Figure 3
T
C
S
A
L
2 comparisons of methods M
i, i
=
1
,
.
.
.
,
8, with respect to δ
b for Application 14 when n
=
800 by Algorithm 11.Application 15 (Phillip’s problem).
We take the following Fredholm integral equation of first kind discussed in [29]. Nowadays this problem is reconsidered by the authors of [41, 42] and solved using a new L-curve technique and a variant of L-curve technique, respectively, to estimate the Tikhonov regularization parameter for regularizing the obtained algebraic linear system.(69)
∫
-
6
6
k
s
,
t
h
t
d
t
=
f
s
.Its solution, kernel, and the right-hand side are given by(70)
h
t
=
1
+
c
o
s
π
3
t
,
if
t
<
3
,
0
,
o
t
h
e
r
w
i
s
e
k
s
,
t
=
h
s
-
t
,
f
s
=
6
-
s
1
+
1
2
cos
⁡
π
3
s
+
9
2
π
sin
⁡
π
3
s
,respectively. We take quadrature nodes t
j
=
-
6
+
12
j
/
n
,
j
=
1
,
…
,
n, evaluating f
(
s
) at the points s
i
=
-
6
+
12
i
/
n
i
=
1
,
…
,
n, and the discretization of (69), (70) gives the algebraic linear system A
u
h
=
b. The trace of the exact solution h
t at the grid points t
j
=
-
6
+
12
j
/
n
,
j
=
1
,
…
,
n, is denoted by u
. The condition number of matrix A
~ (C
o
n
d
(
A
~
)) of Application 15 with respect to δ
A where δ
A
=
0.5
(
δ
b
)
1.5 and the values of δ
b are given in Table 7. T
C
S
A
L
1 with respect to δ
b by using the methods M
i, i
=
1
,
…
,
8, for Phillip’s problem is given in Figure 4 when n
=
800. This figure demonstrates that the proposed method M
3 (17) requires the minimum total solution cost in seconds with respect to the perturbations δ
b by Algorithm 6 to achieve the desired accuracy r
~
m
∞
/
b
~
∞
≤
5
×
1
0
-
7 when n
=
800 for Application 15.
Table8 presents the iteration number m
~ performed, TMMs, T
C
S
A
L
1 for an accuracy of r
~
m
~
j
∞
/
b
~
∞
≤
5
×
1
0
-
7, and the relative L
2 norm of the errors using the methods M
i, i
=
1
,
…
,
8, by Algorithm 6 for Application 15 when n
=
800 and δ
b
=
1
0
-
7. It can be viewed from this table that the best relative L
2 norm error is approximately 1.85
×
1
0
-
4
. To get more smooth solution, we use the methods M
i, i
=
1
,
…
,
8, inStep 4 of Algorithm 6 by applying successive perturbations for solving the perturbed systems (
A
+
δ
A
j
I
)
y
j
=
b
+
Δ
b
j, Δ
b
j
=
[
δ
b
j
,
…
,
δ
b
j
]
T for j
=
1
,
…
,
5, where δ
b
1
=
δ
b
=
1
0
-
7, δ
A
1
=
0.5
δ
b
1
1.5, and δ
b
j
=
0.999
δ
b
j
-
1, δ
A
j
=
0.999
δ
A
j
-
1
,
j
=
2
,
…
,
5, such that the obtained approximate inverse is used as the initial approximate inverse for the next perturbed system. Let, y
~
m
~
j
j
j
=
1
,
…
,
5 be the approximate solutions obtained by Algorithm 6 by performing m
~
j iterations for an accuracy of r
~
m
~
j
∞
/
b
~
∞
≤
5
×
1
0
-
7 for the considered perturbed systems. Total iteration number m
~
=
∑
j
=
1
5
‍
m
~
j performed, the T
C
S
S
P
A
L
1 and TMMs, and the relative L
2 norm of the errors for the last three perturbed systems using the methods M
i, i
=
1
,
…
,
8, by Algorithm 6 for Application 15 when n
=
800 are presented in Table 9. This table shows that best relative L
2 norm error obtained after fifth successive perturbation is 1.33
×
1
0
-
6. Thus, the proposed methods M
3, M
5, and M
7 with Algorithm 6 give solutions by performing less total solution cost in seconds compared with the methods of same orders M
4, M
6, and M
8, respectively, for Application 15. T
C
S
A
L
2 with respect to δ
b by Algorithm 11 applied with the methods M
i, i
=
1
,
…
,
8, for Application 15 when n
=
800 are illustrated in Figure 5. Figure 5 shows that, for the considered values of δ
b, the proposed method M
3 (17) requires the minimum T
C
S
A
L
2 by Algorithm 11 to achieve the desired accuracy T
~
1
,
m
∞
<
0.05 and d
~
l
∞
<
5
×
1
0
-
7, when n
=
800 for Application 15. The P
C
S
A
L
2, I
C
S
A
L
2 and the T
C
S
A
L
2, iteration numbers m
1
∗
,
l
∗ and TBMMs, maximum norm errors, and relative L
2 norm errors with respect to δ
b
=
1
0
-
7 of the methods M
i, i
=
1,2
,
…
,
8, by Algorithm 11, when n
=
800 for Application 15 are shown in Table 10.
Moreover, to compare our results with the existing ones from the literature we taken
=
64 for Phillip’s problem. The best relative error of numerical solution occurring in [42] by Galerkin method with the code Phillips was 1.84
×
1
0
-
3. On the other hand relative maximum error using the proposed methods M
1
,
M
3
,
M
5
,
M
7 by Algorithm 6 when δ
b
=
1
0
-
7 is 1.52
×
1
0
-
4 and by Algorithm 11 when δ
b
=
1
0
-
6 is 6.45
×
1
0
-
4.Table 7
Condition number of matrixA
~ of Application 15, with respect to δ
A
=
0.5
δ
b
1.5 when n
=
800.
δ
b
C
o
n
d
(
A
~
)
1
0
-
1
66631.7
1
0
-
2
2.37544
E
+
07
1
0
-
3
1.00631
E
+
09
1
0
-
4
6.92523
E
+
10
1
0
-
5
2.04433
E
+
10
1
0
-
6
1.06393
E
+
10
1
0
-
7
1.03069
E
+
10Table 8
T
C
S
A
L
1, total iteration numbers and TMMs, and relative L
2 norm of the errors obtained by Algorithm 6 using the methods M
i, i
=
1,2
,
.
.
.
,
8, for the perturbed system of Application 15 when n
=
800 and δ
b
=
1
0
-
7.
Method
T
C
S
A
L
1
m
~
TMMs
u
-
y
~
m
~
2
u
2
M
1
111.52
8
40
4.4350585
E
-
04
M
2
111.52
8
40
4.4350585
E
-
04
M
3
113.78
7
42
2.6701195
E
-
04
M
4
132.03
7
49
2.6701194
E
-
04
M
5
113.53
6
42
3.3382269
E
-
04
M
6
113.86
6
42
3.3382269
E
-
04
M
7
126.18
6
48
1.8498884
E
-
04
M
8
126.19
6
48
1.8498884
E
-
04Table 9
T
C
S
S
P
A
L
1, total iteration numbers and TMMs, and relative L
2 norm of the errors obtained by Algorithm 6 using the methods M
i, i
=
1,2
,
.
.
.
,
8, for the perturbed systems of Application 15 when n
=
800, δ
b
1
=
1
0
-
7.
Method
T
C
S
S
P
A
L
1
m
~
TMMs
u
-
y
~
m
~
3
3
2
u
2
u
-
y
~
m
~
4
4
2
u
2
u
-
y
~
m
~
5
5
2
u
2
M
1
183.03
12
60
8.7677
E
-
05
3.8965
E
-
05
1.7299
E
-
05
M
2
183.03
12
60
8.7677
E
-
05
3.8965
E
-
05
1.7299
E
-
05
M
3
187.46
11
66
3.6207
E
-
05
1.3305
E
-
05
4.8529
E
-
06
M
4
196.85
11
77
3.6207
E
-
05
1.3305
E
-
05
4.8551
E
-
06
M
5
191.51
10
70
3.4970
E
-
05
1.1284
E
-
05
3.6010
E
-
06
M
6
191.52
10
70
3.4970
E
-
05
1.1284
E
-
05
3.5974
E
-
06
M
7
214.06
10
80
1.5882
E
-
05
4.6073
E
-
06
1.3287
E
-
06
M
8
215.97
10
80
1.5882
E
-
05
4.6067
E
-
06
1.3287
E
-
06Table 10
Computational costs, iteration numbers and TBMMs, maximum norm, and relativeL
2 norm of the errors obtained by Algorithm 11 using the methods M
i, i
=
1
,
.
.
.
,
8, for the perturbed system of Application 15, when n
=
800 and δ
b
=
1
0
-
7.
Method
P
C
S
A
L
2
I
C
S
A
L
2
T
C
S
A
L
2
m
1
∗
TBMMs
l
∗
u
-
y
~
l
∗
∞
u
-
y
~
l
∗
2
u
2
M
1
46.15
0.05
46.20
24
120
3
5.7171
E
-
06
7.9054
E
-
06
M
2
46.15
0.05
46.20
24
120
3
5.7565
E
-
06
7.9598
E
-
06
M
3
42.70
0.05
42.75
19
114
4
5.7841
E
-
06
8.0019
E
-
06
M
4
50.12
0.05
50.17
19
133
3
5.7256
E
-
06
7.9330
E
-
06
M
5
43.48
0.05
43.53
17
119
4
5.7543
E
-
06
7.9658
E
-
06
M
6
43.42
0.05
43.47
17
119
3
5.7695
E
-
06
7.9780
E
-
06
M
7
47.26
0.05
47.31
16
128
3
5.7174
E
-
06
7.8968
E
-
06
M
8
47.68
0.05
47.73
16
128
3
5.6986
E
-
06
7.9142
E
-
06Figure 4
T
C
S
A
L
1 comparisons of methods M
i, i
=
1
,
…
,
8, with respect to δ
b for Application 15 by Algorithm 6 when n
=
800.Figure 5
T
C
S
A
L
2 comparisons of methods M
i, i
=
1
,
…
,
8, with respect to δ
b for Application 15 by Algorithm 11.Application 16 (selected matrices from Harwell-Boeing collection).
Test matrices are selected from the sets SMTAPE, GRENOBLE, and NUCL originating from simplex method basis matrix, simulation of computer systems, and nuclear reactor core model respectively, in Harwell-Boeing collection which are available from “Matrix Market”, a repository organized by the National Institute of Standards and Technology. In all these problems right-hand side is generated from a known solution vectoru of all ones. We use the proposed method M
3 inStep 4 of Algorithm 6 by applying successive perturbations for solving the perturbed systems (
A
+
δ
A
j
I
)
y
j
=
b
+
Δ
b
j, Δ
b
j
=
[
δ
b
j
,
…
,
δ
b
j
]
T for j
=
1
,
…
,
6, where δ
b
1
=
δ
b, δ
A
1
=
0.5
δ
b
1
1.5, and δ
b
j
=
0.999
δ
b
j
-
1, δ
A
j
=
0.999
δ
A
j
-
1
,
j
=
2
,
…
,
6, such that the obtained approximate inverse is used as the initial approximate inverse for the next perturbed system. Let y
~
m
~
j
j
j
=
1
,
…
,
6 be the approximate solutions obtained by Algorithm 6 by performing m
~
j iterations for an accuracy of r
~
m
~
j
∞
/
b
~
∞
≤
5
×
δ
b for the considered perturbed systems. The minimal values of the relative L
2 norm of the errors with respect to δ
b using the method M
3 by Algorithm 6 for the test problems and the total iteration number m
~
=
∑
j
=
1
6
‍
m
~
j performed and T
C
S
S
P
A
L
1 are presented in Table 11. The numerical solution of the problems GRE216A, GRE216B, GRE343, GRE512, and GRE1107 is also studied in [43] and the authors provided results obtained by balance scheme in Table III of [43]. The experimental study of ILU preconditioners for the problems BP1000, GRE1107, and NNC1374 is given in [5], and the authors classify these problems hard problems to solve. We present the numerical results of these test problems from the studies [5, 43] in Table 12, of which the Pnc means that problem is not considered, F denotes the failure, and RNEng means relative norm error (RNE) obtained by the method cited in the corresponding reference which is not given. In Table 12 second column presents relative norm errors (RNE) and the iteration numbers (iter) of the problems from GRENOBLE set by balance scheme taken from the Table III of [43], in which the authors mentioned that the preconditioned GMRES failed to converge for any of these problems. This table also demonstrates the iteration numbers for the preconditioned GMRES method or possible reasons of the failure of this method for the solution of the problems BP1000, GRE1107, and NNC1374 with the preconditioners ILU0 taken from Table 3, ILUTP(
30,1.00
) taken from Table 5, and ILUTP(
30,0.01
) taken from Table 6 of [5], presented in third, fourth, and fifth columns of Table 12, respectively. Furthermore, we present the solution of some problems from GRENOBLE set in Application 16 for the minimal values of relative L
2 norm of the errors with respect to δ
b to achieve the desired accuracy T
~
1
,
m
∞
<
0.05 and d
~
l
∞
<
ε using the method M
3 by Algorithm 11 in Table 13. We conclude that the proposed algorithms give stable and highly accurate solution of the considered test problems when compared with the results in Table 12. Moreover, taking into consideration the fact that the authors of [5] classified the problem NNC1374 as very hard problem to solve, Algorithm 6 gave stable and almost accurate solution of this problem.Table 11
T
C
S
S
P
A
L
1, total iteration numbers, and relative L
2 norm errors obtained by Algorithm 6 with method M
3 for the test problems of Application 16.
Problem
m
~
T
C
S
S
P
A
L
1
δ
b
u
-
y
~
m
~
1
1
2
u
2
u
-
y
~
m
~
5
5
2
u
2
u
-
y
~
m
~
6
6
2
u
2
BP1000
20
364.23
1
0
-
10
4.0661
E
-
04
2.5258
E
-
07
2.5232
E
-
07
GRE216A
22
7.62
1
0
-
15
5.7460
E
-
15
4.1467
E
-
15
3.6724
E
-
15
GRE216B
7
2.70
1
0
-
6
1.3477
E
-
06
1.2275
E
-
06
9.0150
E
-
08
GRE343
41
51.59
1
0
-
15
3.7381
E
-
15
3.5990
E
-
15
4.2013
E
-
15
GRE512
11
50.06
1
0
-
12
1.0091
E
-
12
1.0053
E
-
12
1.0036
E
-
12
GRE1107
13
486.25
1
0
-
10
1.6864
E
-
06
1.1541
E
-
06
8.9042
E
-
07
NNC1374
12
1515
1
0
-
7
0.5529
0.5528
0.5524Table 12
Results of test problems from the literature.
[43]
[5]
Problem
Balance scheme
ILU0
ILUTP(
30,1.00
)
ILUTP(
30,0.01
)
(RNE, iter)
(RNE, iter)
(RNE, iter)
(RNE, iter)
BP1000
Pnc
F
R
N
E
n
g
,
13
R
N
E
n
g
,
32
GRE216A
0.14
E
-
04,57
Pnc
Pnc
Pnc
GRE216B
0.39
E
-
06,61
Pnc
Pnc
Pnc
GRE343
0.88
E
-
05,66
Pnc
Pnc
Pnc
GRE512
0.56
E
-
05,80
Pnc
Pnc
Pnc
GRE1107
0.72
E
-
04,200
F(Unstable solve)
F(Inaccuracy)
F(Unstable solve)
NNC1374
Pnc
F(Small pivot)
F(Unstable solve)
F(Unstable solve)Table 13
Computational costs in seconds, iteration numbers, and relativeL
2 norm errors obtained by Algorithm 11 with method M
3 for some problems in Application 16.
Problem
δ
b
ε
P
C
S
A
L
2
I
C
S
A
L
2
T
C
S
A
L
2
m
1
∗
l
∗
u
-
y
~
l
∗
2
u
2
GRE216A
1
0
-
15
5
×
1
0
-
15
0.297
negligible
0.297
6
4
1.1943
E
-
15
GRE216B
1
0
-
18
5
×
1
0
-
6
1.281
0.016
1.297
28
15
1.2494
E
-
06
GRE512
1
0
-
15
5
×
1
0
-
15
5.015
negligible
5.015
8
2
1.3643
E
-
15


## 7. Concluding Remarks

In this study for any integerk
≥
1, we propose a family of methods with recursive structure for computing approximate matrix inverse of a real nonsingular square matrix with convergence order p
=
4
k
+
3
. It is proven that these methods require κ
=
k
+
4, matrix by matrix multiplications per iteration, which are fewer than κ
=
p
=
4
k
+
3, for the standard hyperpower method of same order. The proposed family of methods perform υ
=
2, matrix by matrix additions other than addition with the identity matrix and γ
=
k
+
2, matrix addition by identity per iteration. An algorithm is proposed by constructing 2
×
2 block ILU decomposition based on approximate Schur complement (Schur-BILU) for the coefficient matrix of the algebraic linear system of equations arising from the ill-posed discrete problems with noisy data. From the proposed methods of approximate matrix inversion, the methods of orders p
=
7,11,15,19 are applied for approximating the Schur complement matrices by which the obtained approximate Schur-BILU preconditioner is used to precondition the one step stationary iterative method. This economic computational efficiency is useful in many practical problems such as the solution of first kind Fredholm integral equations and the numerical results justify the given theoretical results. Furthermore, the cost of construction of the approximate Schur-BILU preconditioner can be amortized over systems with same coefficient matrix and different right sides since the preconditioner is to be reused several times. Table 14 shows that Algorithm 11 is at least two times faster than Algorithm 6 used with the methods M
i, i
=
1,2
,
…
,
8, in both considered problems when n
=
800 and δ
b
=
1
0
-
5 for Application 14 and δ
b
=
1
0
-
7 for Application 15.Table 14
Ratios of total solution costs in seconds of Algorithms6 and 11 for both applications using the methods M
i, i
=
1
,
…
,
8, when n
=
800.
Method
T
C
S
A
L
1
T
C
S
A
L
2 for Application 14
T
C
S
A
L
1
T
C
S
A
L
2 for Application 15
M
1
2.45
2.41
M
2
2.45
2.41
M
3
2.37
2.66
M
4
2.44
2.63
M
5
2.53
2.61
M
6
2.54
2.62
M
7
2.41
2.67
M
8
2.37
2.64



---

*Source: 1912535-2019-04-24.xml*