8. The Symmetric Indefinite Method with Orthogonal Factorization¶

8.1. The symmlq Module¶

A Python implementation of SYMMLQ.

This is a line-by-line translation from Matlab code available at http://www.stanford.edu/group/SOL/software/symmlq.htm.

class symmlq.Symmlq(op, **kwargs)¶

Bases: pykrylov.generic.generic.KrylovMethod

SYMMLQ is designed to solve the system of linear equations A x = b where A is an n by n symmetric operator and b is a given vector.

If shift is nonzero, SYMMLQ solves (A - shift I) x = b.

SYMMLQ requires one operator-vector products with A, 2 dot products and 4 daxpys per iteration.

If a preconditioner is supplied, SYMMLQ needs to solve one preconditioning system per iteration. This is a Pythonized line-by-line translation of Michael Saunders’ original SYMMLQ implementation in Matlab

Parameters:	op: an operator describing the coefficient operator A. y = op * x must return the operator-vector product y = Ax for any given vector x.
Keywords:	precon: optional preconditioner. If not None, y = precon * x returns the vector y solution of the linear system M y = x. The preconditioner must be symmetric and positive definite.

Note: atol has no effect in this method.

References:

[PaiSau75]

C. C. Paige and M. A. Saunders, Solution of Sparse Indefinite Systems of Linear Equations, SIAM Journal on Numerical Analysis, 12 (4), pp. 617–629, 1975.

solve(rhs, **kwargs)¶

Solve a linear system with rhs as right-hand side by the SYMMLQ method.

Parameters:	rhs: right-hand side vector b. Should be a Numpy array.
Keywords:	matvec_max: Max. number of matrix-vector produts. Default: 2n+2. rtol: relative stopping tolerance. Default: 1.0e-9. shift: optional shift value. Default: 0. check: specify whether or not to check that matvec indeed describes a symmetric matrix and that precon indeed describes a symmetric positive-definite preconditioner.

8.2. Example¶

Here is an example using SYMMLQ on a linear system. The coefficient matrix is read from file in Matrix Market format:

import numpy as np
from pykrylov.symmlq import SYMMLQ as KSolver
from pykrylov.linop import PysparseLinearOperator
from pysparse import spmatrix
from pysparse.pysparseMatrix import PysparseMatrix as sp

A = sp(matrix=spmatrix.ll_mat_from_mtx('lund_a.mtx'))
n = A.shape[0]
e = np.ones(n)
rhs = A*e

ks = KSolver(PysparseLinearOperator(A))
ks.solve(rhs, matvec_max=2*n)

print 'Number of matvecs: ', ks.nMatvec
print 'Initial/final res: %8.2e/%8.2e' % (ks.residNorm0, ks.residNorm)
print 'Direct error: %8.2e' % (np.linalg.norm(ks.bestSolution-e)/sqrt(n))

Running this script produces the following output:

Number of matvecs:  308
Initial/final res: 1.98e+09/2.83e+01
Direct error: 2.02e-04