Coverage for kwant/solvers/sparse.py : 64%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# Copyright 2011-2013 Kwant authors. # # This file is part of Kwant. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution and at # http://kwant-project.org/license. A list of Kwant authors can be found in # the file AUTHORS.rst at the top-level directory of this distribution and at # http://kwant-project.org/authors.
# Note: previous code would have failed if UMFPACK was provided by scikit
# check if we are actually using UMFPACK or rather SuperLU # TODO: remove the try (only using the except clause) once we depend on # scipy >= 0.14.0.
umfpack = linsolve.umfpack
def factorized(A, piv_tol=1.0, sym_piv_tol=1.0): """ Return a fuction for solving a sparse linear system, with A pre-factorized.
Examples -------- >>> solve = factorized(A) # Makes LU decomposition. >>> x1 = solve(rhs1) # Uses the LU factors. >>> x2 = solve(rhs2) # Uses again the LU factors.
Parameters ---------- A : csc_matrix matrix to be factorized piv_tol : float, 0 <= piv_tol <= 1.0 sym_piv_tol : float, 0 <= piv_tol <= 1.0 thresholds used by UMFPACK for pivoting. 0 means no pivoting, 1.0 means full pivoting as in dense matrices (guaranteeing stability, but reducing possibly sparsity). Defaults of UMFPACK are 0.1 and 0.001 respectively. Whether piv_tol or sym_piv_tol are used is decided internally by UMFPACK, depending on whether the matrix is "symmetric" enough. """
if not sp.isspmatrix_csc(A): A = sp.csc_matrix(A)
A.sort_indices() A = A.asfptype() # upcast to a floating point format
if A.dtype.char not in 'dD': raise ValueError("convert matrix data to double, please, using" " .astype()")
family = {'d': 'di', 'D': 'zi'} umf = umfpack.UmfpackContext(family[A.dtype.char])
# adjust pivot thresholds umf.control[umfpack.UMFPACK_PIVOT_TOLERANCE] = piv_tol umf.control[umfpack.UMFPACK_SYM_PIVOT_TOLERANCE] = sym_piv_tol
# Make LU decomposition. umf.numeric(A)
def solve(b): return umf.solve(umfpack.UMFPACK_A, A, b, autoTranspose=True)
return solve else: # no UMFPACK found. SuperLU is being used, but usually abysmally slow # (SuperLu is not bad per se, somehow the SciPy version isn't good). # Since scipy doesn't include UMFPACK anymore due to software rot, # there is no warning here.
"Sparse Solver class based on the sparse direct solvers provided by SciPy."
|