Source code for pyamg.krylov._fgmres

from warnings import warn
from numpy import array, zeros, sqrt, ravel, abs, max, dot, arange, conjugate, hstack, ceil, isnan, isinf 
from scipy.sparse.linalg.isolve.utils import make_system
from scipy.sparse.sputils import upcast
from pyamg.util.linalg import norm
from pyamg import amg_core
from scipy.linalg import get_blas_funcs 
import scipy

__docformat__ = "restructuredtext en"

__all__ = ['fgmres']

def mysign(x):
    if x == 0.0:
        return 1.0
    else:
        # return the complex "sign"
        return x/abs(x)

[docs]def fgmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None, M=None, callback=None, residuals=None): '''Flexible Generalized Minimum Residual Method (fGMRES) fGMRES iteratively refines the initial solution guess to the system Ax = b. fGMRES is flexible in the sense that the right preconditioner (M) can vary from iteration to iteration. Parameters ---------- A : {array, matrix, sparse matrix, LinearOperator} n x n, linear system to solve b : {array, matrix} right hand side, shape is (n,) or (n,1) x0 : {array, matrix} initial guess, default is a vector of zeros tol : float relative convergence tolerance, i.e. tol is scaled by ||r_0||_2 restrt : {None, int} - if int, restrt is max number of inner iterations and maxiter is the max number of outer iterations - if None, do not restart GMRES, and max number of inner iterations is maxiter maxiter : {None, int} - if restrt is None, maxiter is the max number of inner iterations and GMRES does not restart - if restrt is int, maxiter is the max number of outer iterations, and restrt is the max number of inner iterations xtype : type dtype for the solution, default is automatic type detection M : {array, matrix, sparse matrix, LinearOperator} n x n, inverted preconditioner, i.e. solve A M x = b. M need not be stationary for fgmres callback : function User-supplied function is called after each iteration as callback( ||rk||_2 ), where rk is the current residual vector residuals : list residuals has the residual norm history, including the initial residual, appended to it Returns ------- (xNew, info) xNew : an updated guess to the solution of Ax = b info : halting status of gmres == ============================================= 0 successful exit >0 convergence to tolerance not achieved, return iteration count instead. This value is precisely the order of the Krylov space. <0 numerical breakdown, or illegal input == ============================================= Notes ----- - The LinearOperator class is in scipy.sparse.linalg.interface. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix. A.psolve(..) is still supported as a legacy. - fGMRES allows for non-stationary preconditioners, as opposed to GMRES - For robustness, Householder reflections are used to orthonormalize the Krylov Space Givens Rotations are used to provide the residual norm each iteration Flexibility implies that the right preconditioner, M or A.psolve, can vary from iteration to iteration Examples -------- >>> from pyamg.krylov.fgmres import fgmres >>> from pyamg.util.linalg import norm >>> import numpy >>> from pyamg.gallery import poisson >>> A = poisson((10,10)) >>> b = numpy.ones((A.shape[0],)) >>> (x,flag) = fgmres(A,b, maxiter=2, tol=1e-8) >>> print norm(b - A*x) 6.5428213057 References ---------- .. [1] Yousef Saad, "Iterative Methods for Sparse Linear Systems, Second Edition", SIAM, pp. 151-172, pp. 272-275, 2003 http://www-users.cs.umn.edu/~saad/books.html ''' # Convert inputs to linear system, with error checking A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) dimen = A.shape[0] ## # Ensure that warnings are always reissued from this function import warnings warnings.filterwarnings('always', module='pyamg\.krylov\._fgmres') # Choose type if not hasattr(A, 'dtype'): Atype = upcast(x.dtype, b.dtype) else: Atype = A.dtype if not hasattr(M, 'dtype'): Mtype = upcast(x.dtype, b.dtype) else: Mtype = M.dtype xtype = upcast(Atype, x.dtype, b.dtype, Mtype) # Should norm(r) be kept if residuals == []: keep_r = True else: keep_r = False # Set number of outer and inner iterations if restrt: if maxiter: max_outer = maxiter else: max_outer = 1 if restrt > dimen: warn('Setting number of inner iterations (restrt) to maximum allowed, which is A.shape[0] ') restrt = dimen max_inner = restrt else: max_outer = 1 if maxiter > dimen: warn('Setting number of inner iterations (maxiter) to maximum allowed, which is A.shape[0] ') maxiter = dimen elif maxiter == None: maxiter = min(dimen, 40) max_inner = maxiter # Get fast access to underlying BLAS routines [rotg] = get_blas_funcs(['rotg'], [x]) # Is this a one dimensional matrix? if dimen == 1: entry = ravel(A*array([1.0], dtype=xtype)) return (postprocess(b/entry), 0) # Prep for method r = b - ravel(A*x) normr = norm(r) if keep_r: residuals.append(normr) # Check initial guess ( scaling by b, if b != 0, # must account for case when norm(b) is very small) normb = norm(b) if normb == 0.0: normb = 1.0 if normr < tol*normb: if callback != None: callback(norm(r)) return (postprocess(x), 0) # Scale tol by ||r_0||_2, we don't use the preconditioned # residual because this is right preconditioned GMRES. if normr != 0.0: tol = tol*normr # Use separate variable to track iterations. If convergence fails, we cannot # simply report niter = (outer-1)*max_outer + inner. Numerical error could cause # the inner loop to halt while the actual ||r|| > tol. niter = 0 # Begin fGMRES for outer in range(max_outer): # Calculate vector w, which defines the Householder reflector # Take shortcut in calculating, # w = r + sign(r[1])*||r||_2*e_1 w = r beta = mysign(w[0])*normr w[0] += beta w /= norm(w) # Preallocate for Krylov vectors, Householder reflectors and Hessenberg matrix # Space required is O(dimen*max_inner) Q = zeros( (4*max_inner,), dtype=xtype) # Givens Rotations H = zeros( (max_inner, max_inner), dtype=xtype) # upper Hessenberg matrix (actually made upper tri with Givens Rotations) W = zeros( (max_inner, dimen), dtype=xtype) # Householder reflectors Z = zeros( (dimen, max_inner), dtype=xtype) # For fGMRES, preconditioned vectors must be stored # No Horner-like scheme exists that allow us to avoid this W[0,:] = w # Multiply r with (I - 2*w*w.T), i.e. apply the Householder reflector # This is the RHS vector for the problem in the Krylov Space g = zeros((dimen,), dtype=xtype) g[0] = -beta for inner in range(max_inner): # Calculate Krylov vector in two steps # (1) Calculate v = P_j = (I - 2*w*w.T)v, where k = inner v = -2.0*conjugate(w[inner])*w v[inner] += 1.0 # (2) Calculate the rest, v = P_1*P_2*P_3...P_{j-1}*ej. #for j in range(inner-1,-1,-1): # v = v - 2.0*dot(conjugate(W[j,:]), v)*W[j,:] amg_core.apply_householders(v, ravel(W), dimen, inner-1, -1, -1) #Apply preconditioner v = ravel(M*v) ## Check for nan, inf #if isnan(v).any() or isinf(v).any(): # warn('inf or nan after application of preconditioner') # return(postprocess(x), -1) Z[:,inner] = v # Calculate new search direction v = ravel(A*v) # Factor in all Householder orthogonal reflections on new search direction #for j in range(inner+1): # v = v - 2.0*dot(conjugate(W[j,:]), v)*W[j,:] amg_core.apply_householders(v, ravel(W), dimen, 0, inner+1, 1) # Calculate next Householder reflector, w # w = v[inner+1:] + sign(v[inner+1])*||v[inner+1:]||_2*e_{inner+1) # Note that if max_inner = dimen, then this is unnecessary for the last inner # iteration, when inner = dimen-1. Here we do not need to calculate a # Householder reflector or Givens rotation because nnz(v) is already the # desired length, i.e. we do not need to zero anything out. if inner != dimen-1: if inner < (max_inner-1): w = W[inner+1,:] vslice = v[inner+1:] alpha = norm(vslice) if alpha != 0: alpha = mysign(vslice[0])*alpha # We do not need the final reflector for future calculations if inner < (max_inner-1): w[inner+1:] = vslice w[inner+1] += alpha w /= norm(w) # Apply new reflector to v # v = v - 2.0*w*(w.T*v) v[inner+1] = -alpha v[inner+2:] = 0.0 if inner > 0: # Apply all previous Givens Rotations to v amg_core.apply_givens(Q, v, dimen, inner) # Calculate the next Givens rotation, where j = inner # Note that if max_inner = dimen, then this is unnecessary for the last inner # iteration, when inner = dimen-1. Here we do not need to calculate a # Householder reflector or Givens rotation because nnz(v) is already the # desired length, i.e. we do not need to zero anything out. if inner != dimen-1: if v[inner+1] != 0: [c,s] = rotg(v[inner], v[inner+1]) Qblock = array([[c, s], [-conjugate(s), c]], dtype=xtype) Q[(inner*4) : ((inner+1)*4)] = ravel(Qblock).copy() # Apply Givens Rotation to g, # the RHS for the linear system in the Krylov Subspace. # Note that this dot does a matrix multiply, not an actual # dot product where a conjugate transpose is taken g[inner:inner+2] = dot(Qblock, g[inner:inner+2]) # Apply effect of Givens Rotation to v v[inner] = dot(Qblock[0,:], v[inner:inner+2]) v[inner+1] = 0.0 # Write to upper Hessenberg Matrix, # the LHS for the linear system in the Krylov Subspace H[:,inner] = v[0:max_inner] # Don't update normr if last inner iteration, because # normr is calculated directly after this loop ends. if inner < max_inner-1: normr = abs(g[inner+1]) if normr < tol: break # Allow user access to residual if callback != None: callback( normr ) if keep_r: residuals.append(normr) niter += 1 # end inner loop, back to outer loop # Find best update to x in Krylov Space, V. Solve inner+1 x inner+1 system. # Apparently this is the best way to solve a triangular # system in the magical world of scipy #piv = arange(inner+1) #y = lu_solve((H[0:(inner+1),0:(inner+1)], piv), g[0:(inner+1)], trans=0) y = scipy.linalg.solve(H[0:(inner+1),0:(inner+1)], g[0:(inner+1)]) # No Horner like scheme exists because the preconditioner can change each iteration # Hence, we must store each preconditioned vector update = dot(Z[:,0:inner+1], y) x = x + update r = b - ravel(A*x) normr = norm(r) # Allow user access to residual if callback != None: callback( normr ) if keep_r: residuals.append(normr) # Has fGMRES stagnated? indices = (x != 0) if indices.any(): change = max(abs( update[indices] / x[indices] )) if change < 1e-12: # No change, halt return (postprocess(x), -1) # test for convergence if normr < tol: return (postprocess(x),0) # end outer loop return (postprocess(x), niter)