(A, B=None, epsilon=4.0, k=2,
                                     proj_type='l2', block_flag=False,
                                     symmetrize_measure=True)
return Atilde
Evolution Strength Measure.

    Construct strength of connection matrix using an Evolution-based measure

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        Sparse NxN matrix
    B : string, array
        If B=None, then the near nullspace vector used is all ones.  If B is
        an (NxK) array, then B is taken to be the near nullspace vectors.
    epsilon : scalar
        Drop tolerance
    k : integer
        ODE num time steps, step size is assumed to be 1/rho(DinvA)
    proj_type : {'l2','D_A'}
        Define norm for constrained min prob, i.e. define projection
    block_flag : boolean
        If True, use a block D inverse as preconditioner for A during
        weighted-Jacobi

    Returns
    -------
    Atilde : csr_matrix
        Sparse matrix of strength values

    See [2008OlScTu]_ for more details.

    References
    ----------
    .. [2008OlScTu] Olson, L. N., Schroder, J., Tuminaro, R. S.,
       "A New Perspective on Strength Measures in Algebraic Multigrid",
       submitted, June, 2008.

    Examples
    --------
    >>> import numpy as np
    >>> from pyamg.gallery import stencil_grid
    >>> from pyamg.strength import evolution_strength_of_connection
    >>> n=3
    >>> stencil =  np.array([[-1.0,-1.0,-1.0],
    ...                        [-1.0, 8.0,-1.0],
    ...                        [-1.0,-1.0,-1.0]])
    >>> A = stencil_grid(stencil, (n,n), format='csr')
    >>> S = evolution_strength_of_connection(A,  np.ones((A.shape[0],1)))
Evolution Strength Measure.
Evolution Strength Measure .
def evolution_strength_of_connection(A, B=None, epsilon=4.0, k=2,
                                     proj_type='l2', block_flag=False,
                                     symmetrize_measure=True):
    """Evolution Strength Measure.

    Construct strength of connection matrix using an Evolution-based measure

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        Sparse NxN matrix
    B : string, array
        If B=None, then the near nullspace vector used is all ones.  If B is
        an (NxK) array, then B is taken to be the near nullspace vectors.
    epsilon : scalar
        Drop tolerance
    k : integer
        ODE num time steps, step size is assumed to be 1/rho(DinvA)
    proj_type : {'l2','D_A'}
        Define norm for constrained min prob, i.e. define projection
    block_flag : boolean
        If True, use a block D inverse as preconditioner for A during
        weighted-Jacobi

    Returns
    -------
    Atilde : csr_matrix
        Sparse matrix of strength values

    See [2008OlScTu]_ for more details.

    References
    ----------
    .. [2008OlScTu] Olson, L. N., Schroder, J., Tuminaro, R. S.,
       "A New Perspective on Strength Measures in Algebraic Multigrid",
       submitted, June, 2008.

    Examples
    --------
    >>> import numpy as np
    >>> from pyamg.gallery import stencil_grid
    >>> from pyamg.strength import evolution_strength_of_connection
    >>> n=3
    >>> stencil =  np.array([[-1.0,-1.0,-1.0],
    ...                        [-1.0, 8.0,-1.0],
    ...                        [-1.0,-1.0,-1.0]])
    >>> A = stencil_grid(stencil, (n,n), format='csr')
    >>> S = evolution_strength_of_connection(A,  np.ones((A.shape[0],1)))

    """
    # ====================================================================
    # Check inputs
    if epsilon < 1.0:
        raise ValueError('expected epsilon > 1.0')
    if k <= 0:
        raise ValueError('number of time steps must be > 0')
    if proj_type not in ['l2', 'D_A']:
        raise ValueError('proj_type must be "l2" or "D_A"')
    if (not sparse.isspmatrix_csr(A)) and (not sparse.isspmatrix_bsr(A)):
        raise TypeError('expected csr_matrix or bsr_matrix')

    # ====================================================================
    # Format A and B correctly.
    # B must be in mat format, this isn't a deep copy
    if B is None:
        Bmat = np.ones((A.shape[0], 1), dtype=A.dtype)
    else:
        Bmat = np.asarray(B)

    # Pre-process A.  We need A in CSR, to be devoid of explicit 0's and have
    # sorted indices
    if not sparse.isspmatrix_csr(A):
        csrflag = False
        numPDEs = A.blocksize[0]
        D = A.diagonal()
        # Calculate Dinv*A
        if block_flag:
            Dinv = get_block_diag(A, blocksize=numPDEs, inv_flag=True)
            Dinv = sparse.bsr_matrix((Dinv, np.arange(Dinv.shape[0]),
                                      np.arange(Dinv.shape[0] + 1)),
                                     shape=A.shape)
            Dinv_A = (Dinv * A).tocsr()
        else:
            Dinv = np.zeros_like(D)
            mask = (D != 0.0)
            Dinv[mask] = 1.0 / D[mask]
            Dinv[D == 0] = 1.0
            Dinv_A = scale_rows(A, Dinv, copy=True)
        A = A.tocsr()
    else:
        csrflag = True
        numPDEs = 1
        D = A.diagonal()
        Dinv = np.zeros_like(D)
        mask = (D != 0.0)
        Dinv[mask] = 1.0 / D[mask]
        Dinv[D == 0] = 1.0
        Dinv_A = scale_rows(A, Dinv, copy=True)

    A.eliminate_zeros()
    A.sort_indices()

    # Handle preliminaries for the algorithm
    dimen = A.shape[1]
    NullDim = Bmat.shape[1]

    # Get spectral radius of Dinv*A, this will be used to scale the time step
    # size for the ODE
    rho_DinvA = approximate_spectral_radius(Dinv_A)

    # Calculate D_A for later use in the minimization problem
    if proj_type == 'D_A':
        D_A = sparse.spdiags([D], [0], dimen, dimen, format='csr')
    else:
        D_A = sparse.eye(dimen, dimen, format='csr', dtype=A.dtype)

    # Calculate (I - delta_t Dinv A)^k
    #      In order to later access columns, we calculate the transpose in
    #      CSR format so that columns will be accessed efficiently
    # Calculate the number of time steps that can be done by squaring, and
    # the number of time steps that must be done incrementally
    nsquare = int(np.log2(k))
    ninc = k - 2**nsquare

    # Calculate one time step
    Id = sparse.eye(dimen, dimen, format='csr', dtype=A.dtype)
    Atilde = (Id - (1.0 / rho_DinvA) * Dinv_A)
    Atilde = Atilde.T.tocsr()

    # Construct a sparsity mask for Atilde that will restrict Atilde^T to the
    # nonzero pattern of A, with the added constraint that row i of Atilde^T
    # retains only the nonzeros that are also in the same PDE as i.
    mask = A.copy()

    # Restrict to same PDE
    if numPDEs > 1:
        row_length = np.diff(mask.indptr)
        my_pde = np.mod(np.arange(dimen), numPDEs)
        my_pde = np.repeat(my_pde, row_length)
        mask.data[np.mod(mask.indices, numPDEs) != my_pde] = 0.0
        del row_length, my_pde
        mask.eliminate_zeros()

    # If the total number of time steps is a power of two, then there is
    # a very efficient computational short-cut.  Otherwise, we support
    # other numbers of time steps, through an inefficient algorithm.
    if ninc > 0:
        warn('The most efficient time stepping for the Evolution Strength '
             f'Method is done in powers of two.\nYou have chosen {k} time steps.')

        # Calculate (Atilde^nsquare)^T = (Atilde^T)^nsquare
        for _i in range(nsquare):
            Atilde = Atilde * Atilde

        JacobiStep = (Id - (1.0 / rho_DinvA) * Dinv_A).T.tocsr()
        for _i in range(ninc):
            Atilde = Atilde * JacobiStep
        del JacobiStep

        # Apply mask to Atilde, zeros in mask have already been eliminated at
        # start of routine.
        mask.data[:] = 1.0
        Atilde = Atilde.multiply(mask)
        Atilde.eliminate_zeros()
        Atilde.sort_indices()

    elif nsquare == 0:
        if numPDEs > 1:
            # Apply mask to Atilde, zeros in mask have already been eliminated
            # at start of routine.
            mask.data[:] = 1.0
            Atilde = Atilde.multiply(mask)
            Atilde.eliminate_zeros()
            Atilde.sort_indices()

    else:
        # Use computational short-cut for case (ninc == 0) and (nsquare > 0)
        # Calculate Atilde^k only at the sparsity pattern of mask.
        for _i in range(nsquare - 1):
            Atilde = Atilde * Atilde

        # Call incomplete mat-mat mult
        AtildeCSC = Atilde.tocsc()
        AtildeCSC.sort_indices()
        mask.sort_indices()
        Atilde.sort_indices()
        amg_core.incomplete_mat_mult_csr(Atilde.indptr, Atilde.indices,
                                         Atilde.data, AtildeCSC.indptr,
                                         AtildeCSC.indices, AtildeCSC.data,
                                         mask.indptr, mask.indices, mask.data,
                                         dimen)

        del AtildeCSC, Atilde
        Atilde = mask
        Atilde.eliminate_zeros()
        Atilde.sort_indices()

    del Dinv, Dinv_A, mask

    # Calculate strength based on constrained min problem of
    # min( z - B*x ), such that
    # (B*x)|_i = z|_i, i.e. they are equal at point i
    # z = (I - (t/k) Dinv A)^k delta_i
    #
    # Strength is defined as the relative point-wise approx. error between
    # B*x and z.  We don't use the full z in this problem, only that part of
    # z that is in the sparsity pattern of A.
    #
    # Can use either the D-norm, and inner product, or l2-norm and inner-prod
    # to solve the constrained min problem.  Using D gives scale invariance.
    #
    # This is a quadratic minimization problem with a linear constraint, so
    # we can build a linear system and solve it to find the critical point,
    # i.e. minimum.
    #
    # We exploit a known shortcut for the case of NullDim = 1.  The shortcut is
    # mathematically equivalent to the longer constrained min. problem

    if NullDim == 1:
        # Use shortcut to solve constrained min problem if B is only a vector
        # Strength(i,j) = | 1 - (z(i)/b(j))/(z(j)/b(i)) |
        # These ratios can be calculated by diagonal row and column scalings

        # Create necessary vectors for scaling Atilde
        #   Its not clear what to do where B == 0.  This is an
        #   an easy programming solution, that may make sense.
        Bmat_forscaling = np.ravel(Bmat)
        Bmat_forscaling[Bmat_forscaling == 0] = 1.0
        DAtilde = Atilde.diagonal()
        DAtildeDivB = np.ravel(DAtilde) / Bmat_forscaling

        # Calculate best approximation, z_tilde, in span(B)
        #   Importantly, scale_rows and scale_columns leave zero entries
        #   in the matrix.  For previous implementations this was useful
        #   because we assume data and Atilde.data are the same length below
        data = Atilde.data.copy()
        Atilde.data[:] = 1.0
        Atilde = scale_rows(Atilde, DAtildeDivB)
        Atilde = scale_columns(Atilde, np.ravel(Bmat_forscaling))

        # If angle in the complex plane between z and z_tilde is
        # greater than 90 degrees, then weak.  We can just look at the
        # dot product to determine if angle is greater than 90 degrees.
        angle = np.multiply(np.real(Atilde.data), np.real(data)) +\
            np.multiply(np.imag(Atilde.data), np.imag(data))
        angle = angle < 0.0
        angle = np.array(angle, dtype=bool)

        # Calculate Approximation ratio
        Atilde.data = Atilde.data / data

        # If approximation ratio is less than tol, then weak connection
        weak_ratio = (np.abs(Atilde.data) < 1e-4)

        # Calculate Approximation error
        Atilde.data = abs(1.0 - Atilde.data)

        # Set small ratios and large angles to weak
        Atilde.data[weak_ratio] = 0.0
        Atilde.data[angle] = 0.0

        # Set near perfect connections to 1e-4
        Atilde.eliminate_zeros()
        Atilde.data[Atilde.data < np.sqrt(np.finfo(float).eps)] = 1e-4

        del data, weak_ratio, angle

    else:
        # For use in computing local B_i^H*B, precompute the element-wise
        # multiply of each column of B with each other column.  We also scale
        # by 2.0 to account for BDB's eventual use in a constrained
        # minimization problem
        BDBCols = int(np.sum(np.arange(NullDim + 1)))
        BDB = np.zeros((dimen, BDBCols), dtype=A.dtype)
        counter = 0
        for i in range(NullDim):
            for j in range(i, NullDim):
                BDB[:, counter] = 2.0 *\
                    (np.conjugate(np.ravel(Bmat[:, i])) * np.ravel(D_A * Bmat[:, j]))
                counter = counter + 1

        # Choose tolerance for dropping "numerically zero" values later
        tol = set_tol(Atilde.dtype)

        # Use constrained min problem to define strength
        amg_core.evolution_strength_helper(Atilde.data,
                                           Atilde.indptr,
                                           Atilde.indices,
                                           Atilde.shape[0],
                                           np.ravel(Bmat),
                                           np.ravel((D_A * B.conj()).T),
                                           np.ravel(BDB),
                                           BDBCols, NullDim, tol)

        Atilde.eliminate_zeros()

    # All of the strength values are real by this point, so ditch the complex
    # part
    Atilde.data = np.array(np.real(Atilde.data), dtype=float)

    # Apply drop tolerance
    if epsilon != np.inf:
        amg_core.apply_distance_filter(dimen, epsilon, Atilde.indptr,
                                       Atilde.indices, Atilde.data)
        Atilde.eliminate_zeros()

    # Symmetrize
    if symmetrize_measure:
        Atilde = 0.5 * (Atilde + Atilde.T)

    # Set diagonal to 1.0, as each point is strongly connected to itself.
    Id = sparse.eye(dimen, dimen, format='csr')
    Id.data -= Atilde.diagonal()
    Atilde = Atilde + Id

    # If converted BSR to CSR, convert back and return amalgamated matrix,
    #   i.e. the sparsity structure of the blocks of Atilde
    if not csrflag:
        Atilde = Atilde.tobsr(blocksize=(numPDEs, numPDEs))

        n_blocks = Atilde.indices.shape[0]
        blocksize = Atilde.blocksize[0] * Atilde.blocksize[1]
        CSRdata = np.zeros((n_blocks,))
        amg_core.min_blocks(n_blocks, blocksize,
                            np.ravel(np.asarray(Atilde.data)), CSRdata)
        # Atilde = sparse.csr_matrix((data, row, col), shape=(*,*))
        Atilde = sparse.csr_matrix((CSRdata, Atilde.indices, Atilde.indptr),
                                   shape=(int(Atilde.shape[0] / numPDEs),
                                          int(Atilde.shape[1] / numPDEs)))

    # Standardized strength values require small values be weak and large
    # values be strong.  So, we invert the algebraic distances computed here
    Atilde.data = 1.0 / Atilde.data

    # Scale C by the largest magnitude entry in each row
    Atilde = scale_rows_by_largest_entry(Atilde)

    return Atilde
def evolution_strength_of_connection ( A , B = None , epsilon = 4.0 , k = 2 , proj_type = 'l2' , block_flag = False , symmetrize_measure = True ) : # ==================================================================== # Check inputs if epsilon < 1.0 : raise ValueError ( 'expected epsilon > 1.0' ) if k <= 0 : raise ValueError ( 'number of time steps must be > 0' ) if proj_type not in [ 'l2' , 'D_A' ] : raise ValueError ( 'proj_type must be "l2" or "D_A"' ) if ( not sparse . isspmatrix_csr ( A ) ) and ( not sparse . isspmatrix_bsr ( A ) ) : raise TypeError ( 'expected csr_matrix or bsr_matrix' ) # ==================================================================== # Format A and B correctly. # B must be in mat format, this isn't a deep copy if B is None : Bmat = np . ones ( ( A . shape [ 0 ] , 1 ) , dtype = A . dtype ) else : Bmat = np . asarray ( B ) # Pre-process A.  We need A in CSR, to be devoid of explicit 0's and have # sorted indices if not sparse . isspmatrix_csr ( A ) : csrflag = False numPDEs = A . blocksize [ 0 ] D = A . diagonal ( ) # Calculate Dinv*A if block_flag : Dinv = get_block_diag ( A , blocksize = numPDEs , inv_flag = True ) Dinv = sparse . bsr_matrix ( ( Dinv , np . arange ( Dinv . shape [ 0 ] ) , np . arange ( Dinv . shape [ 0 ] + 1 ) ) , shape = A . shape ) Dinv_A = ( Dinv * A ) . tocsr ( ) else : Dinv = np . zeros_like ( D ) mask = ( D != 0.0 ) Dinv [ mask ] = 1.0 / D [ mask ] Dinv [ D == 0 ] = 1.0 Dinv_A = scale_rows ( A , Dinv , copy = True ) A = A . tocsr ( ) else : csrflag = True numPDEs = 1 D = A . diagonal ( ) Dinv = np . zeros_like ( D ) mask = ( D != 0.0 ) Dinv [ mask ] = 1.0 / D [ mask ] Dinv [ D == 0 ] = 1.0 Dinv_A = scale_rows ( A , Dinv , copy = True ) A . eliminate_zeros ( ) A . sort_indices ( ) # Handle preliminaries for the algorithm dimen = A . shape [ 1 ] NullDim = Bmat . shape [ 1 ] # Get spectral radius of Dinv*A, this will be used to scale the time step # size for the ODE rho_DinvA = approximate_spectral_radius ( Dinv_A ) # Calculate D_A for later use in the minimization problem if proj_type == 'D_A' : D_A = sparse . spdiags ( [ D ] , [ 0 ] , dimen , dimen , format = 'csr' ) else : D_A = sparse . eye ( dimen , dimen , format = 'csr' , dtype = A . dtype ) # Calculate (I - delta_t Dinv A)^k #      In order to later access columns, we calculate the transpose in #      CSR format so that columns will be accessed efficiently # Calculate the number of time steps that can be done by squaring, and # the number of time steps that must be done incrementally nsquare = int ( np . log2 ( k ) ) ninc = k - 2 ** nsquare # Calculate one time step Id = sparse . eye ( dimen , dimen , format = 'csr' , dtype = A . dtype ) Atilde = ( Id - ( 1.0 / rho_DinvA ) * Dinv_A ) Atilde = Atilde . T . tocsr ( ) # Construct a sparsity mask for Atilde that will restrict Atilde^T to the # nonzero pattern of A, with the added constraint that row i of Atilde^T # retains only the nonzeros that are also in the same PDE as i. mask = A . copy ( ) # Restrict to same PDE if numPDEs > 1 : row_length = np . diff ( mask . indptr ) my_pde = np . mod ( np . arange ( dimen ) , numPDEs ) my_pde = np . repeat ( my_pde , row_length ) mask . data [ np . mod ( mask . indices , numPDEs ) != my_pde ] = 0.0 del row_length , my_pde mask . eliminate_zeros ( ) # If the total number of time steps is a power of two, then there is # a very efficient computational short-cut.  Otherwise, we support # other numbers of time steps, through an inefficient algorithm. if ninc > 0 : warn ( 'The most efficient time stepping for the Evolution Strength ' f'Method is done in powers of two.\nYou have chosen {k} time steps.' ) # Calculate (Atilde^nsquare)^T = (Atilde^T)^nsquare for _i in range ( nsquare ) : Atilde = Atilde * Atilde JacobiStep = ( Id - ( 1.0 / rho_DinvA ) * Dinv_A ) . T . tocsr ( ) for _i in range ( ninc ) : Atilde = Atilde * JacobiStep del JacobiStep # Apply mask to Atilde, zeros in mask have already been eliminated at # start of routine. mask . data [ : ] = 1.0 Atilde = Atilde . multiply ( mask ) Atilde . eliminate_zeros ( ) Atilde . sort_indices ( ) elif nsquare == 0 : if numPDEs > 1 : # Apply mask to Atilde, zeros in mask have already been eliminated # at start of routine. mask . data [ : ] = 1.0 Atilde = Atilde . multiply ( mask ) Atilde . eliminate_zeros ( ) Atilde . sort_indices ( ) else : # Use computational short-cut for case (ninc == 0) and (nsquare > 0) # Calculate Atilde^k only at the sparsity pattern of mask. for _i in range ( nsquare - 1 ) : Atilde = Atilde * Atilde # Call incomplete mat-mat mult AtildeCSC = Atilde . tocsc ( ) AtildeCSC . sort_indices ( ) mask . sort_indices ( ) Atilde . sort_indices ( ) amg_core . incomplete_mat_mult_csr ( Atilde . indptr , Atilde . indices , Atilde . data , AtildeCSC . indptr , AtildeCSC . indices , AtildeCSC . data , mask . indptr , mask . indices , mask . data , dimen ) del AtildeCSC , Atilde Atilde = mask Atilde . eliminate_zeros ( ) Atilde . sort_indices ( ) del Dinv , Dinv_A , mask # Calculate strength based on constrained min problem of # min( z - B*x ), such that # (B*x)|_i = z|_i, i.e. they are equal at point i # z = (I - (t/k) Dinv A)^k delta_i # # Strength is defined as the relative point-wise approx. error between # B*x and z.  We don't use the full z in this problem, only that part of # z that is in the sparsity pattern of A. # # Can use either the D-norm, and inner product, or l2-norm and inner-prod # to solve the constrained min problem.  Using D gives scale invariance. # # This is a quadratic minimization problem with a linear constraint, so # we can build a linear system and solve it to find the critical point, # i.e. minimum. # # We exploit a known shortcut for the case of NullDim = 1.  The shortcut is # mathematically equivalent to the longer constrained min. problem if NullDim == 1 : # Use shortcut to solve constrained min problem if B is only a vector # Strength(i,j) = | 1 - (z(i)/b(j))/(z(j)/b(i)) | # These ratios can be calculated by diagonal row and column scalings # Create necessary vectors for scaling Atilde #   Its not clear what to do where B == 0.  This is an #   an easy programming solution, that may make sense. Bmat_forscaling = np . ravel ( Bmat ) Bmat_forscaling [ Bmat_forscaling == 0 ] = 1.0 DAtilde = Atilde . diagonal ( ) DAtildeDivB = np . ravel ( DAtilde ) / Bmat_forscaling # Calculate best approximation, z_tilde, in span(B) #   Importantly, scale_rows and scale_columns leave zero entries #   in the matrix.  For previous implementations this was useful #   because we assume data and Atilde.data are the same length below data = Atilde . data . copy ( ) Atilde . data [ : ] = 1.0 Atilde = scale_rows ( Atilde , DAtildeDivB ) Atilde = scale_columns ( Atilde , np . ravel ( Bmat_forscaling ) ) # If angle in the complex plane between z and z_tilde is # greater than 90 degrees, then weak.  We can just look at the # dot product to determine if angle is greater than 90 degrees. angle = np . multiply ( np . real ( Atilde . data ) , np . real ( data ) ) + np . multiply ( np . imag ( Atilde . data ) , np . imag ( data ) ) angle = angle < 0.0 angle = np . array ( angle , dtype = bool ) # Calculate Approximation ratio Atilde . data = Atilde . data / data # If approximation ratio is less than tol, then weak connection weak_ratio = ( np . abs ( Atilde . data ) < 1e-4 ) # Calculate Approximation error Atilde . data = abs ( 1.0 - Atilde . data ) # Set small ratios and large angles to weak Atilde . data [ weak_ratio ] = 0.0 Atilde . data [ angle ] = 0.0 # Set near perfect connections to 1e-4 Atilde . eliminate_zeros ( ) Atilde . data [ Atilde . data < np . sqrt ( np . finfo ( float ) . eps ) ] = 1e-4 del data , weak_ratio , angle else : # For use in computing local B_i^H*B, precompute the element-wise # multiply of each column of B with each other column.  We also scale # by 2.0 to account for BDB's eventual use in a constrained # minimization problem BDBCols = int ( np . sum ( np . arange ( NullDim + 1 ) ) ) BDB = np . zeros ( ( dimen , BDBCols ) , dtype = A . dtype ) counter = 0 for i in range ( NullDim ) : for j in range ( i , NullDim ) : BDB [ : , counter ] = 2.0 * ( np . conjugate ( np . ravel ( Bmat [ : , i ] ) ) * np . ravel ( D_A * Bmat [ : , j ] ) ) counter = counter + 1 # Choose tolerance for dropping "numerically zero" values later tol = set_tol ( Atilde . dtype ) # Use constrained min problem to define strength amg_core . evolution_strength_helper ( Atilde . data , Atilde . indptr , Atilde . indices , Atilde . shape [ 0 ] , np . ravel ( Bmat ) , np . ravel ( ( D_A * B . conj ( ) ) . T ) , np . ravel ( BDB ) , BDBCols , NullDim , tol ) Atilde . eliminate_zeros ( ) # All of the strength values are real by this point, so ditch the complex # part Atilde . data = np . array ( np . real ( Atilde . data ) , dtype = float ) # Apply drop tolerance if epsilon != np . inf : amg_core . apply_distance_filter ( dimen , epsilon , Atilde . indptr , Atilde . indices , Atilde . data ) Atilde . eliminate_zeros ( ) # Symmetrize if symmetrize_measure : Atilde = 0.5 * ( Atilde + Atilde . T ) # Set diagonal to 1.0, as each point is strongly connected to itself. Id = sparse . eye ( dimen , dimen , format = 'csr' ) Id . data -= Atilde . diagonal ( ) Atilde = Atilde + Id # If converted BSR to CSR, convert back and return amalgamated matrix, #   i.e. the sparsity structure of the blocks of Atilde if not csrflag : Atilde = Atilde . tobsr ( blocksize = ( numPDEs , numPDEs ) ) n_blocks = Atilde . indices . shape [ 0 ] blocksize = Atilde . blocksize [ 0 ] * Atilde . blocksize [ 1 ] CSRdata = np . zeros ( ( n_blocks , ) ) amg_core . min_blocks ( n_blocks , blocksize , np . ravel ( np . asarray ( Atilde . data ) ) , CSRdata ) # Atilde = sparse.csr_matrix((data, row, col), shape=(*,*)) Atilde = sparse . csr_matrix ( ( CSRdata , Atilde . indices , Atilde . indptr ) , shape = ( int ( Atilde . shape [ 0 ] / numPDEs ) , int ( Atilde . shape [ 1 ] / numPDEs ) ) ) # Standardized strength values require small values be weak and large # values be strong.  So, we invert the algebraic distances computed here Atilde . data = 1.0 / Atilde . data # Scale C by the largest magnitude entry in each row Atilde = scale_rows_by_largest_entry ( Atilde ) return Atilde