Source code for pyamg.classical.interpolate

"""Classical AMG Interpolation methods"""


__docformat__ = "restructuredtext en"

import numpy
from scipy.sparse import csr_matrix, isspmatrix_csr
from pyamg import amg_core

__all__ = ['direct_interpolation']


[docs]def direct_interpolation(A, C, splitting): """Create prolongator using direct interpolation Parameters ---------- A : {csr_matrix} NxN matrix in CSR format C : {csr_matrix} Strength-of-Connection matrix Must have zero diagonal splitting : array C/F splitting stored in an array of length N Returns ------- P : {csr_matrix} Prolongator using direct interpolation Examples -------- >>> from pyamg.gallery import poisson >>> from pyamg.classical import direct_interpolation >>> import numpy >>> A = poisson((5,),format='csr') >>> splitting = numpy.array([1,0,1,0,1], dtype='intc') >>> P = direct_interpolation(A, A, splitting) >>> print P.todense() [[ 1. 0. 0. ] [ 0.5 0.5 0. ] [ 0. 1. 0. ] [ 0. 0.5 0.5] [ 0. 0. 1. ]] """ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix for A') if not isspmatrix_csr(C): raise TypeError('expected csr_matrix for C') # Interpolation weights are computed based on entries in A, but subject to # the sparsity pattern of C. So, copy the entries of A into the # sparsity pattern of C. C = C.copy() C.data[:] = 1.0 C = C.multiply(A) Pp = numpy.empty_like(A.indptr) amg_core.rs_direct_interpolation_pass1(A.shape[0], C.indptr, C.indices, splitting, Pp) nnz = Pp[-1] Pj = numpy.empty(nnz, dtype=Pp.dtype) Px = numpy.empty(nnz, dtype=A.dtype) amg_core.rs_direct_interpolation_pass2(A.shape[0], A.indptr, A.indices, A.data, C.indptr, C.indices, C.data, splitting, Pp, Pj, Px) return csr_matrix((Px, Pj, Pp))