linSolDefCg¶

Solves a real symmetric definite linear system using a conjugate gradient method. Using optional arguments, a preconditioner can be supplied.

Synopsis¶

linSolDefCg (amultp, b)

Required Arguments¶

amultp (p, z): User-supplied function which computes \(z = Ap\).
float b (Input): Vector of length n containing the right-hand side.

Return Value¶

The solution x of the linear system \(Ax = b\). If no solution was computed, then None is returned.

Optional Arguments¶

maxIter (Input/Output)

An integer, initially set to the maximum number of iterations allowed. On exit, the number of iterations used is returned.

relErr, float (Input)

The relative error desired.

Default: relErr = sqrt(machine(4))

precond (r, z) (Input)

User supplied function which sets \(z = M^{-1}r\), where M is the preconditioning matrix.

jacobi, float[] (Input)

Use the Jacobi preconditioner, i.e. M = diag (A). The user-supplied vector diagonal should be set so that jacobi [i] = \(A_{ii}\).

Description¶

The function linSolDefCg solves the symmetric definite linear system \(Ax = b\) using the conjugate gradient method with optional preconditioning. This method is described in detail by Golub and Van Loan (1983, Chapter 10), and in Hageman and Young (1981, Chapter 7).

The preconditioning matrix M is a matrix that approximates A, and for which the linear system \(Mz = r\) is easy to solve. These two properties are in conflict; balancing them is a topic of much current research. In the default use of linSolDefCg, \(M = I\). If the option jacobi is selected, M is set to the diagonal of A.

The number of iterations needed depends on the matrix and the error tolerance. As a rough guide,

\[\mathtt{maxIter} = \sqrt{n} \text{ for } n >> 1\]

See the references mentioned above for details.

Let M be the preconditioning matrix, let b, p, r, x, and z be vectors and let τ be the desired relative error. Then the algorithm used is as follows:

\[\begin{split}\begin{array}{c} \lambda = -1 \hfill\\ p_0 = x_0 \hfill \\ r_1 = b - Ap \hfill \\ \text{for } k = 1, \ldots, \text{maxIter} \hfill \\ \phantom{....} z_k = M^{-1}r_k \hfill \\ \phantom{....} \text{if } k = 1, \text{then} \hfill \\ \phantom{........} \beta_k = 1 \hfill \\ \phantom{........} p_k = z_k \hfill \\ \phantom{....} \text{else} \hfill \\ \phantom{........} \beta_k = \left(z_k^Tr_k\right) / \left(z_{k-1}^Tr_{k-1}\right) \hfill \\ \phantom{........} p_k = z_k + \beta_kp_k \hfill \\ \phantom{....} \text{endif} \hfill \\ \phantom{....} z_k = Ap \hfill \\ \phantom{....} \alpha_k = \left(z_{k-1}^Tz_{k-1}\right) / \left(z_k^Tp_k\right) \hfill \\ \phantom{....} x_k = x_k + \alpha_kp_k \hfill \\ \phantom{....} r_k = r_k - \alpha_kz_k \hfill \\ \phantom{....} \text{if } \left(\|z_k\|_2 \leq \tau\left(1-\lambda\right)\|x_k\|_2\right) \text{then} \hfill \\ \phantom{........} \text{recompute } \lambda \hfill \\ \phantom{........} \text{if } \left(\|z_k\|_2 \leq \tau\left(1-\lambda\right)\|x_k\|_2\right) \text{exit} \hfill \\ \phantom{....} \text{endif} \hfill \\ \text{endfor} \hfill \\ \end{array}\end{split}\]

Here λ is an estimate of \(λ_{max} G\), the largest eigenvalue of the iteration matrix \(G = I - M^{-1}A\). The stopping criterion is based on the result (Hageman and Young 1981, pp. 148-151)

\[\frac{\|x_k-x\|_M}{\|x\|_M} \leq \left(\frac{1}{1-\lambda_{\max}(G)}\right) \left(\frac{\|z_k\|_M}{\|x_k\|_M}\right)\]

where

\[\|x\|_M^2 = x^TMx\]

It is also known that

\[\lambda_{\max}\left(T_1\right) \leq \lambda_{\max}\left(T_2\right) \leq \ldots \leq \lambda_{\max}(G) < 1\]

where the \(T_n\) are the symmetric, tridiagonal matrices

\[\begin{split}T_n = \begin{bmatrix} \mu_1 & \omega_2 & & \\ \omega_2 & \mu_2 & \omega_3 & \\ & \omega_3 & \mu_3 & \ddots \\ & & \ddots & \ddots \\ \end{bmatrix}\end{split}\]

with \(μ_k = 1 - β_k / α_{k-1} - 1/α_k\), \(μ_1 = 1 - 1/α_1\) and

\[\omega_k = \sqrt{B_k} / \alpha_{k-1}\]

Usually the eigenvalue computation is needed for only a few of the iterations.

Examples¶

Example 1¶

In this example, the solution to a linear system is found. The coefficient matrix is stored as a full matrix.

from numpy import *
from pyimsl.util.imslConvert import toCtypes
from pyimsl.math.linSolDefCg import linSolDefCg
from pyimsl.math.matMulRect import matMulRect
from pyimsl.math.writeMatrix import writeMatrix
from pyimsl.util.imslConvert import toNdarray

n = 3
a = [[1.0, -3.0, 2.0],
     [-3.0, 10.0, -5.0],
     [2.0, -5.0, 6.0]]


def amultp(p, z):
    ptmp = toNdarray(p, (n))
    ztmp = matMulRect("A*x",
                      aMatrix=a,
                      xVector=ptmp)
    toCtypes(ztmp, z)


b = [27.0, -78.0, 64.0]
x = linSolDefCg(amultp, b)
writeMatrix("x", x)

Output¶

 
                  x
          1            2            3
          1           -4            7

Example 2¶

In this example, two different preconditioners are used to find the solution of a linear system which occurs in a finite difference solution of Laplace’s equation on a regular c × c grid, c = 100. The matrix is \(A = E\) (\(c^2\), c). For the first solution, select Jacobi preconditioning and supply the diagonal, so M = diag (A). The number of iterations performed and the maximum absolute error are printed. Next, use a more complicated preconditioning matrix, M, consisting of the symmetric tridiagonal part of A.

Notice that the symmetric positive definite band solver is used to factor M once, and subsequently just perform forward and back solves. Again, the number of iterations performed and the maximum absolute error are printed. Note the substantial reduction in iterations.

from __future__ import print_function
from numpy import *
from pyimsl.util.imslConvert import toCtypes
from pyimsl.math.generateTestBand import generateTestBand
from pyimsl.math.generateTestCoordinate import generateTestCoordinate
from pyimsl.math.linSolDefCg import linSolDefCg
from pyimsl.math.linSolPosdefBand import linSolPosdefBand
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
from pyimsl.math.vectorNorm import vectorNorm
from pyimsl.util.imslConvert import toNdarray

c = 50
n = 2500
b = empty((n), dtype=double)

# Generate coefficient matrix
a = generateTestCoordinate(n, c)
aMatrix = {'nrowa': n, 'ncola': n, 'a': a}


def amultp(p, z):
    ptmp = toNdarray(p, (n))
    ztmp = matMulRectCoordinate("A*x",
                                aMatrix=aMatrix,
                                xVector=ptmp)
    toCtypes(ztmp, z)


def precond1(r, z):
    # Factor the first time through
    factor = []
    m = generateTestBand(n, 1, symmetricStorage=True)
    res = linSolPosdefBand(m, 1, None, factor=factor,
                           factorOnly=True)

    tempr = toNdarray(r, n)
    # Perform the forward and back solves
    ztmp = linSolPosdefBand(m, 1, tempr,
                            factor=factor,
                            solveOnly=True)
    toCtypes(ztmp, z)


index = []
maxit = [1000]

# Set a predetermined answer and diagonal
mod_five = empty((n), dtype=double)
diagonal = empty((n), dtype=double)
for i in range(0, n):
    mod_five[i] = (i % 5)
    diagonal[i] = 4.0

# Get right hand side
b = matMulRectCoordinate("A*x",
                         aMatrix=aMatrix,
                         xVector=mod_five)

# Solve with jacobi preconditioning
x = linSolDefCg(amultp, b,
                maxIter=maxit,
                jacobi=diagonal)

# Find max absolute error, print results
norm = vectorNorm(x, secondVector=mod_five,
                  infNorm=index)
print("iterations = %d, norm = %e" % (maxit[0], norm))

# Solve same system, with different preconditioner
maxit = [1000]
x = linSolDefCg(amultp, b,
                maxIter=maxit,
                precond=precond1)
norm = vectorNorm(x, secondVector=mod_five, infNorm=index)
print("iterations = %d, norm = %e" % (maxit[0], norm))

Output¶

iterations = 187, norm = 4.463445e-10
iterations = 127, norm = 5.134553e-10

Fatal Errors¶

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.