C/C++模板：线性系统解法与迭代方法

需积分: 9 109 浏览量更新于2024-07-31 收藏 598KB PDF 举报

"Templates for the Solution of Linear Systems 是一本关于使用C/C++模板解决线性系统的书籍，重点介绍迭代法在矩阵计算中的应用。作者团队包括多位知名专家，该书由SIAM出版，并受到多个机构的支持。" 在矩阵计算领域，求解线性系统是一个核心问题，通常涉及到大量的数值计算和算法实现。这本书"Templates for the Solution of Linear Systems"提供了一种利用C/C++模板来高效处理这类问题的方法。模板是一种强大的C++特性，允许在编译时进行类型检查和代码生成，这使得代码更灵活且易于复用。书中强调了迭代法在求解线性系统中的应用。迭代法是一种非直接求解线性系统Ax=b的方法，它通过不断地近似来逐渐提高解的精度，而不是一次性求出精确解。常见的迭代方法包括：雅可比迭代、高斯-塞德尔迭代、共轭梯度法（CG）、最小二乘 QR 迭代等。这些方法在处理大规模稀疏矩阵时特别有效，因为它们避免了大矩阵的直接乘法，降低了计算复杂性和内存需求。分块技术是另一种优化线性系统求解的重要策略。在处理大型矩阵时，可以将矩阵分成较小的子块，分别处理后再组合，这样可以减少计算量并提高计算效率。此外，分块技术还可以与迭代法结合，如分块雅可比或分块共轭梯度法，以进一步提升性能。线性系统在工程、科学和许多其他领域有广泛应用，例如在物理模拟、图像处理、机器学习和数据科学中。理解如何有效地使用迭代法和分块技术是任何从事相关工作的人必须掌握的关键技能。这本书不仅提供了理论基础，还包含实际的C/C++模板代码，使得读者可以直接在自己的项目中应用这些方法。 "Templates for the Solution of Linear Systems"是一本宝贵的资源，它不仅教导读者如何用现代编程工具解决复杂的数学问题，还通过实例展示了如何将理论转化为实际代码，这对于数值计算和高性能计算领域的研究人员和工程师来说极具价值。

6 CHAPTER 2. ITERATIVE METHODS

– SSOR.

Symmetric Successive Overrelaxation (SSOR) has no advantage over SOR as a stand-

alone iterative method; however, it is useful as a preconditioner for nonstationary meth-

ods.

• Nonstationary Methods

– Conjugate Gradient (CG).

The conjugate gradient method derives its name from the fact that it generates a se-

quence of conjugate (or orthogonal) vectors. These vectors are the residuals of the iter-

ates. They are also the gradients of a quadratic functional, the minimization of which

is equivalent to solving the linear system. CG is an extremely effective method when

the coefﬁcient matrix is symmetric positive deﬁnite, since storage for only a limited

number of vectors is required.

– Minimum Residual (MINRES) and Symmetric LQ (SYMMLQ).

These methods are computational alternatives for CG for coefﬁcient matrices that are

symmetric but possibly indeﬁnite. SYMMLQ will generate the same solution iterates

as CG if the coefﬁcient matrix is symmetric positive deﬁnite.

– Conjugate Gradient on the Normal Equations: CGNE and CGNR.

These methods are based on the application of the CG method to one of two forms of

the normal equations for Ax = b. CGNE solves the system (AA

)y = b for y and then

computes the solution x = A

y. CGNR solves (A

A)x =

b for the solution vector

x where

b = A

b. When the coefﬁcient matrix A is nonsymmetric and nonsingular,

the normal equations matrices AA

and A

A will be symmetric and positive deﬁnite,

and hence CG can be applied. The convergence may be slow, since the spectrum of the

normal equations matrices will be less favorable than the spectrum of A.

– Generalized Minimal Residual (GMRES).

The Generalized Minimal Residual method computes a sequence of orthogonal vectors

(like MINRES), and combines these through a least-squares solve and update. How-

ever, unlike MINRES (and CG) it requires storing the whole sequence, so that a large

amount of storage is needed. For this reason, restarted versions of this method are used.

In restarted versions, computation and storage costs are limited by specifying a ﬁxed

number of vectors to be generated. This method is useful for general nonsymmetric

matrices.

– BiConjugate Gradient (BiCG).

The Biconjugate Gradient method generates two CG-like sequences of vectors, one

based on a system with the original coefﬁcient matrix A, and one on A

. Instead of

orthogonalizing each sequence, they are made mutually orthogonal, or “bi-orthogonal”.

This method, like CG, uses limited storage. It is useful when the matrix is nonsymmet-

ric and nonsingular; however, convergence may be irregular, and there is a possibility

that the method will break down. BiCG requires a multiplication with the coefﬁcient

matrix and with its transpose at each iteration.

– Quasi-Minimal Residual (QMR).

The Quasi-Minimal Residual method applies a least-squares solve and update to the

BiCG residuals, thereby smoothing out the irregular convergence behavior of BiCG.

Also, QMR largely avoids the breakdown that can occur in BiCG. On the other hand,

it does not effect a true minimization of either the error or the residual, and while it

converges smoothly, it does not essentially improve on the BiCG.

– Conjugate Gradient Squared (CGS).

2.2. STATIONARY ITERATIVE METHODS 7

The Conjugate Gradient Squared method is a variant of BiCG that applies the updating

operations for the A-sequence and the A

-sequences both to the same vectors. Ideally,

this would double the convergence rate, but in practice convergence may be much more

irregular than for BiCG. A practical advantage is that the method does not need the

multiplications with the transpose of the coefﬁcient matrix.

– Biconjugate Gradient Stabilized (Bi-CGSTAB).

The Biconjugate Gradient Stabilized method is a variant of BiCG, like CGS, but using

different updates for the A

-sequence in order to obtain smoother convergence than

CGS.

– Chebyshev Iteration.

The Chebyshev Iteration recursively determines polynomials with coefﬁcients chosen

to minimize the norm of the residual in a min-max sense. The coefﬁcient matrix must

be positive deﬁnite and knowledge of the extremal eigenvalues is required. This method

has the advantage of requiring no inner products.

2.2 Stationary Iterative Methods

Iterative methods that can be expressed in the simple form

(k)

= Bx

(k−1)

+ c (2.1)

(where neither B nor c depend upon the iteration count k) are called stationary iterative methods.

In this section, we present the four main stationary iterative methods: the Jacobi method, the

Gauss-Seidel method, the Successive Overrelaxation (SOR) method and the Symmetric Successive

Overrelaxation (SSOR) method. In each case, we summarize their convergence behavior and their

effectiveness, and discuss how and when they should be used. Finally, in §2.2.5, we give some

historical background and further notes and references.

2.2.1 The Jacobi Method

The Jacobi method is easily derived by examining each of the n equations in the linear system

Ax = b in isolation. If in the ith equation

j=1

i,j

= b

we solve for the value of x

while assuming the other entries of x remain ﬁxed, we obtain

= (b

−

j6=i

i,j

)/a

i,i

. (2.2)

This suggests an iterative method deﬁned by

(k)

= (b

−

j6=i

i,j

(k−1)

)/a

i,i

, (2.3)

which is the Jacobi method. Note that the order in which the equations are examined is irrelevant,

since the Jacobi method treats them independently. For this reason, the Jacobi method is also

known as the method of simultaneous displacements, since the updates could in principle be done

simultaneously.

In matrix terms, the deﬁnition of the Jacobi method in (2.3) can be expressed as

(k)

= D

−1

(L + U )x

(k−1)

+ D

−1

b, (2.4)

8 CHAPTER 2. ITERATIVE METHODS

Choose an initial guess x

(0)

to the solution x.

for k = 1, 2, . . .

for i = 1, 2, . . . , n

¯x

= 0

for j = 1, 2, . . . , i − 1, i + 1, . . . , n

¯x

= ¯x

+ a

i,j

(k−1)

end

¯x

= (b

− ¯x

)/a

i,i

end

(k)

= ¯x

check convergence; continue if necessary

end

Figure 2.1: The Jacobi Method

where the matrices D, −L and −U represent the diagonal, the strictly lower-triangular, and the

strictly upper-triangular parts of A, respectively.

The pseudocode for the Jacobi method is given in Figure 2.1. Note that an auxiliary storage

vector, ¯x is used in the algorithm. It is not possible to update the vector x in place, since values

from x

(k−1)

are needed throughout the computation of x

(k)

Convergence of the Jacobi method

Iterative methods are often used for solving discretized partial differential equations. In that context

a rigorous analysis of the convergence of simple methods such as the Jacobi method can be given.

As an example, consider the boundary value problem

Lu = −u

= f on (0, 1), u(0) = u

, u(1) = u

discretized by

Lu(x

) = 2u(x

) − u(x

i−1

) − u(x

i+1

) = f(x

)/N

for x

= i/N, i = 1 . . . N − 1.

The eigenfunctions of the Land L operator are the same: for n = 1 . . . N − 1 the function u

(x) =

sin nπx is an eigenfunction corresponding to λ = 4 sin

nπ/(2N). The eigenvalues of the Jacobi

iteration matrix B are then λ(B) = 1 − 1/2λ(L) = 1 −2 sin

nπ/(2N).

From this it is easy to see that the high frequency modes (i.e., eigenfunction u

with n large)

are damped quickly, whereas the damping factor for modes with n small is close to 1. The spectral

radius of the Jacobi iteration matrix is ≈ 1 − 10/N

, and it is attained for the eigenfunction

u(x) = sin πx.

The type of analysis applied to this example can be generalized to higher dimensions and other

stationary iterative methods. For both the Jacobi and Gauss-Seidel method (below) the spectral

radius is found to be 1 − O(h

) where h is the discretization mesh width, i.e., h = N

−d

where N

is the number of variables and d is the number of space dimensions.

2.2.2 The Gauss-Seidel Method

Consider again the linear equations in (2.2). If we proceed as with the Jacobi method, but now

assume that the equations are examined one at a time in sequence, and that previously computed

2.2. STATIONARY ITERATIVE METHODS 9

Choose an initial guess x

(0)

to the solution x.

for k = 1, 2, . . .

for i = 1, 2, . . . , n

σ = 0

for j = 1, 2, . . . , i − 1

σ = σ + a

i,j

(k)

end

for j = i + 1, . . . , n

σ = σ + a

i,j

(k−1)

end

(k)

= (b

− σ)/a

i,i

end

check convergence; continue if necessary

end

Figure 2.2: The Gauss-Seidel Method

results are used as soon as they are available, we obtain the Gauss-Seidel method:

(k)

= (b

−

j<i

i,j

(k)

−

j>i

i,j

(k−1)

)/a

i,i

. (2.5)

Two important facts about the Gauss-Seidel method should be noted. First, the computations

in (2.5) appear to be serial. Since each component of the new iterate depends upon all previ-

ously computed components, the updates cannot be done simultaneously as in the Jacobi method.

Second, the new iterate x

(k)

depends upon the order in which the equations are examined. The

Gauss-Seidel method is sometimes called the method of successive displacements to indicate the

dependence of the iterates on the ordering. If this ordering is changed, the components of the new

iterate (and not just their order) will also change.

These two points are important because if A is sparse, the dependency of each component

of the new iterate on previous components is not absolute. The presence of zeros in the matrix

may remove the inﬂuence of some of the previous components. Using a judicious ordering of

the equations, it may be possible to reduce such dependence, thus restoring the ability to make

updates to groups of components in parallel. However, reordering the equations can affect the rate

at which the Gauss-Seidel method converges. A poor choice of ordering can degrade the rate of

convergence; a good choice can enhance the rate of convergence. For a practical discussion of this

tradeoff (parallelism versus convergence rate) and some standard reorderings, the reader is referred

to Chapter 3 and §4.4.

In matrix terms, the deﬁnition of the Gauss-Seidel method in (2.5) can be expressed as

(k)

= (D − L)

−1

(Ux

(k−1)

+ b). (2.6)

As before, D, −L and −U represent the diagonal, lower-triangular, and upper-triangular parts of A,

respectively.

The pseudocode for the Gauss-Seidel algorithm is given in Figure 2.2.

2.2.3 The Successive Overrelaxation Method

The Successive Overrelaxation Method, or SOR, is devised by applying extrapolation to the Gauss-

Seidel method. This extrapolation takes the form of a weighted average between the previous iterate

10 CHAPTER 2. ITERATIVE METHODS

Choose an initial guess x

(0)

to the solution x.

for k = 1, 2, . . .

for i = 1, 2, . . . , n

σ = 0

for j = 1, 2, . . . , i − 1

σ = σ + a

i,j

(k)

end

for j = i + 1, . . . , n

σ = σ + a

i,j

(k−1)

end

σ = (b

− σ)/a

i,i

(k)

= x

(k−1)

+ ω(σ − x

(k−1)

)

end

check convergence; continue if necessary

end

Figure 2.3: The SOR Method

and the computed Gauss-Seidel iterate successively for each component:

(k)

= ω¯x

(k)

+ (1 − ω)x

(k−1)

(where ¯x denotes a Gauss-Seidel iterate, and ω is the extrapolation factor). The idea is to choose a

value for ω that will accelerate the rate of convergence of the iterates to the solution.

In matrix terms, the SOR algorithm can be written as follows:

(k)

= (D − ωL)

−1

(ωU + (1 − ω)D)x

(k−1)

+ ω(D −ωL)

−1

b. (2.7)

The pseudocode for the SOR algorithm is given in Figure 2.3.

Choosing the Value of ω

If ω = 1, the SOR method simpliﬁes to the Gauss-Seidel method. A theorem due to Kahan [129]

shows that SOR fails to converge if ω is outside the interval (0, 2). Though technically the term

underrelaxation should be used when 0 < ω < 1, for convenience the term overrelaxation is now

used for any value of ω ∈ (0, 2).

In general, it is not possible to compute in advance the value of ω that is optimal with respect

to the rate of convergence of SOR. Even when it is possible to compute the optimal value for ω, the

expense of such computation is usually prohibitive. Frequently, some heuristic estimate is used,

such as ω = 2 −O(h) where h is the mesh spacing of the discretization of the underlying physical

domain.

If the coefﬁcient matrix A is symmetric and positive deﬁnite, the SOR iteration is guaranteed to

converge for any value of ω between 0 and 2, though the choice of ω can signiﬁcantly affect the rate

at which the SOR iteration converges. Sophisticated implementations of the SOR algorithm (such

as that found in ITPACK [139]) employ adaptive parameter estimation schemes to try to home in

on the appropriate value of ω by estimating the rate at which the iteration is converging.

For coefﬁcient matrices of a special class called consistently ordered with property A (see

Young [215]), which includes certain orderings of matrices arising from the discretization of el-

liptic PDEs, there is a direct relationship between the spectra of the Jacobi and SOR iteration

剩余116页未读，继续阅读

tuzbing

粉丝: 0
资源: 6

C/C++模板：线性系统解法与迭代方法

Templates for the Solution of Linear Systems - Building Blocks for Iterative Methods.ps

Solution of Linear Systems

稳定的双共轭梯度法算法设计与matlab程序

CUDA Templates for Linear Algebra Subroutines and Solvers - headers only

spring.thymeleaf.prefix=classpath:/templates/ spring.thymeleaf.suffix=.html spring.thymeleaf.cache=false spring.thymeleaf.mode=LEGACYHTML5

c++ templates - the complete guide, 2nd edition pdf github

use C++ to write a program to enable input and output of a vector, respectively，Define a class template called Vector(a single-column- Matrix). The templates can instantiate a vector of any element type. Overloaded >> and << operators: to enable input and output of a vector, respectively

masm for windows

c++ templates the complete guide 2nd pdf

dart Live Templates

最新资源