Direct Solution of Symmetric Systems of Equations

The pyma27 Module

Ma27: Direct multifrontal solution of symmetric systems.

class hsl.solvers.pyma27.PyMa27Solver(A, **kwargs)

Bases: hsl.solvers.sils.Sils

fetch_perm()

Return the permutation vector p.

The permutation vector is used to compute the factorization of A. Rows and columns were permuted so that

Pᵀ A P = L B Lᵀ

where i-th row of P is the p(i)-th row of the identity matrix, L is unit upper triangular and B is block diagonal with 1x1 and 2x2 blocks.

inertia

Return inertia of matrix A.

Inertia of a matrix is given by (# of λ>0, # of λ<0, # λ=0).

refine(b, nitref=3, tol=1e-08, **kwargs)

Perform iterative refinement if necessary.

Iterative refinement is performed until the scaled residual norm ‖b-Ax‖/(1+‖b‖) falls below the threshold ‘tol’ or until nitref steps are taken. Make sure you have called solve() with the same right-hand side b before calling refine(). The residual vector self.residual will be updated to reflect the updated approximate solution.

By default, tol = 1.0e-8 and nitref = 3.

solve(b, get_resid=True)

Solve the linear system of equations Ax = b.

The solution will be found in self.x and residual in self.residual.

The pyma57 Module

Ma57: Direct multifrontal solution of symmetric systems.

class hsl.solvers.pyma57.PyMa57Solver(A, factorize=True, **kwargs)

Bases: hsl.solvers.sils.Sils

factorize(A)

Perform numerical factorization.

Before this can be done, symbolic factorization (the analyze phase) must have been performed.

The values of the elements of the matrix may have been altered since the analyze phase but the sparsity pattern must not have changed. Use the optional argument newA to specify the updated matrix if applicable.

fetch_perm()

Return the permutation vector p.

Permutation vector is used to compute the factorization of A. Rows and columns were permuted so that

Pᵀ A P = L B Lᵀ

where i-th row of P is the p(i)-th row of the identity matrix, L is unit upper triangular and B is block diagonal with 1x1 and 2x2 blocks.

inertia

Return inertia of matrix A.

Inertia of a matrix is given by (# of λ>0, # of λ<0, # λ=0).

refine(b, nitref=3, **kwargs)

Perform iterative refinement if necessary.

Iterative refinement is performed until the scaled residual norm ‖b-Ax‖/(1+‖b‖) falls below the threshold ‘tol’ or until nitref steps are taken. Make sure you have called solve() with the same right-hand side b before calling refine(). The residual vector self.residual will be updated to reflect the updated approximate solution.

By default, nitref = 3.

solve(b, get_resid=True)

Solve the linear system of equations Ax = b.

The solution will be found in self.x and residual in self.residual.

Example

"""Exemple from MA57 spec sheet: http://www.hsl.rl.ac.uk/specs/ma57.pdf."""

import sys
import numpy as np
from hsl.solvers.src._cyma57_numpy_INT32_FLOAT64 import NumpyMA57Solver_INT32_FLOAT64

n = 5
nnz = 7
A = np.array([[2.0, 3.0, 0, 0, 0], [0, 0, 4.0, 0, 6.0], [0, 0, 1, 5, 0],
              [0, 0, 0, 0, 0], [0, 0, 0, 0, 1]], dtype=np.float32)
arow = np.array([0, 0, 1, 1, 2, 2, 4], dtype=np.int32)
acol = np.array([0, 1, 2, 4, 2, 3, 4], dtype=np.int32)
aval = np.array([2.0, 3.0, 4.0, 6.0, 1.0, 5.0, 1.0], dtype=np.float64)

rhs = np.array([8, 45, 31, 15, 17], dtype=np.float64)


context = NumpyMA57Solver_INT32_FLOAT64(n, n, nnz)
context.get_matrix_data(arow, acol, aval)

context.analyze()

context.factorize()

print 'Solve:'
x, residual = context.solve(rhs, True)
# x = context.solve(rhs, False)
print '  x:'
print x
print '  residual:'
print residual

print 'Fetch_perm:'
perm = context.fetch_perm()
print '  perm:'
print perm

print 'Refine:'
(new_x, new_res) = context.refine(x, rhs, residual, 5)
print '  new_x: '
print new_x
print '  new_res: '
print new_res

from cysparse.sparse.ll_mat import *
import cysparse.types.cysparse_types as types
import numpy as np

from hsl.solvers.src._cyma57_cysparse_INT32_FLOAT64 import CySparseMA57Solver_INT32_FLOAT64
import sys


A = LLSparseMatrix(mm_filename=sys.argv[1], itype=types.INT32_T,
                   dtype=types.FLOAT64_T)

(n, m) = A.shape
e = np.ones(n, 'd')
# rhs = np.zeros(n, 'd')
rhs = A * e

context = CySparseMA57Solver_INT32_FLOAT64(A.nrow, A.ncol, A.nnz)
context.get_matrix_data(A)

context.analyze()

context.factorize()

print 'Fetch_perm:'
perm = context.fetch_perm()
print '  perm:'
print perm


print 'Solve:'
x, residual = context.solve(rhs, True)
print '  x:'
print x
print '  residual:'
print residual

print 'Refine:'
(new_x, new_res) = context.refine(x, rhs, residual, 5)
print '  cond: ', context.cond
print '  new_x: '
print new_x
print '  new_res: '
print new_res

# Demonstrate usage of PyMa27 abstract class for the solution of symmetric
# systems of linear equations.
# Example usage: python demo_ma27.py file1.mtx ... fileN.mtx
# where each fileK.mtx is in MatrixMarket format.

# from hsl.solvers.pyma27 import PyMa27Solver as LBLContext
from hsl.solvers.pyma57 import PyMa57Solver as LBLContext
from pysparse import spmatrix
from numpy.linalg import norm
import timeit
import numpy as np


def Hilbert(n):
    """The cream of ill conditioning: the Hilbert matrix.

    See Higham, "Accuracy and Stability of Numerical Algoriths", section 28.1.
    The matrix has elements H(i,j) = 1/(i+j-1) when indexed i,j=1..n.  However,
    here we index as i,j=0..n-1, so the elements are H(i,j) = 1/(i+j+1).
    """
    if n <= 0:
        return None
    if n == 1:
        return 1.0
    nnz = n * (n - 1) / 2
    H = spmatrix.ll_mat_sym(n, nnz)
    for i in range(n):
        for j in range(i + 1):
            H[i, j] = 1.0 / (i + j + 1)
    return H


def Ma27SpecSheet():
    """This is the example from the MA27 spec sheet.

    Solution should be [1,2,3,4,5]"""
    A = spmatrix.ll_mat_sym(5, 7)
    A[0, 0] = 2
    A[1, 0] = 3
    A[2, 1] = 4
    A[2, 2] = 1
    A[3, 2] = 5
    A[4, 1] = 6
    A[4, 4] = 1

    rhs = np.ones(5, 'd')
    rhs[0] = 8
    rhs[1] = 45
    rhs[2] = 31
    rhs[3] = 15
    rhs[4] = 17

    return (A, rhs)


def solve_system(A, rhs, itref_threshold=1.0e-6, nitrefmax=5, **kwargs):

    # Obtain Sils context object
    t = timeit.default_timer()
    LBL = LBLContext(A, **kwargs)
    t_analyze = timeit.default_timer() - t

    # Solve system and compute residual
    t = timeit.default_timer()
    LBL.solve(rhs)
    t_solve = timeit.default_timer() - t

    # Compute residual norm
    nrhsp1 = norm(rhs, ord=np.infty) + 1
    nr = norm(LBL.residual) / nrhsp1

    # If residual not small, perform iterative refinement
    LBL.refine(rhs, tol=itref_threshold, nitref=nitrefmax)
    nr1 = norm(LBL.residual, ord=np.infty) / nrhsp1

    return (LBL.x, LBL.residual, nr, nr1, t_analyze, t_solve, LBL.neig)

if __name__ == '__main__':
    import sys
    import os
    matrices = sys.argv[1:]

    hdr_fmt = '%-13s  %-11s  %-11s  %-11s  %-7s  %-7s  %-5s\n'
    res_fmt = '%-13s  %-11.5e  %-11.5e  %-11.5e  %-7.3f  %-7.3f  %-5d\n'
    hdrs = ('Name', 'Rel. resid.', 'Residual', 'Resid itref', 'Analyze',
            'Solve', 'neig')
    header = hdr_fmt % hdrs
    lhead = len(header)
    sys.stderr.write('-' * lhead + '\n')
    sys.stderr.write(header)
    sys.stderr.write('-' * lhead + '\n')

    # Solve example from the spec sheet
    (A, rhs) = Ma27SpecSheet()
    (x, r, nr, nr1, t_an, t_sl, neig) = solve_system(A, rhs)
    exact = np.arange(5, dtype='d') + 1
    relres = norm(x - exact) / norm(exact)
    sys.stdout.write(res_fmt %
                     ('Spec sheet', relres, nr, nr1, t_an, t_sl, neig))

    # Solve example with Hilbert matrix
    n = 10
    H = Hilbert(n)
    e = np.ones(n, 'd')
    rhs = np.empty(n, 'd')
    H.matvec(e, rhs)
    (x, r, nr, nr1, t_an, t_sl, neig) = solve_system(H, rhs)
    relres = norm(x - e) / norm(e)
    sys.stdout.write(res_fmt % ('Hilbert', relres, nr, nr1, t_an, t_sl, neig))

    # Process matrices given on the command line
    for matrix in matrices:
        M = spmatrix.ll_mat_from_mtx(matrix)
        (m, n) = M.shape
        if m != n:
            break
        e = np.ones(n, 'd')
        rhs = np.empty(n, 'd')
        M.matvec(e, rhs)
        (x, r, nr, nr1, t_an, t_sl, neig) = solve_system(M, rhs)
        relres = norm(x - e) / norm(e)
        probname = os.path.basename(matrix)
        if probname[-4:] == '.mtx':
            probname = probname[:-4]
        sys.stdout.write(res_fmt %
                         (probname, relres, nr, nr1, t_an, t_sl, neig))
    sys.stderr.write('-' * lhead + '\n')