def rho_D_inv_A(A):
    """Return the (approx.) spectral radius of D^-1 * A.

    Parameters
    ----------
    A : sparse-matrix

    Returns
    -------
    approximate spectral radius of diag(A)^{-1} A

    Examples
    --------
    >>> from pyamg.gallery import poisson
    >>> from pyamg.relaxation.smoothing import rho_D_inv_A
    >>> from scipy.sparse import csr_matrix
    >>> import numpy as np
    >>> A = csr_matrix(np.array([[1.0,0,0],[0,2.0,0],[0,0,3.0]]))
    >>> print rho_D_inv_A(A)
    1.0

    """
    if not hasattr(A, 'rho_D_inv'):
        D_inv = get_diagonal(A, inv=True)
        D_inv_A = scale_rows(A, D_inv, copy=True)
        A.rho_D_inv = approximate_spectral_radius(D_inv_A)

    return A.rho_D_inv

def reference_classical_soc(A, theta, norm='abs'):
    """
    This complex extension of the classic Ruge-Stuben
    strength-of-connection has some theoretical justification in
    "AMG Solvers for Complex-Valued Matrices", Scott MacClachlan,
    Cornelis Oosterlee

    A connection is strong if,
      | a_ij| >= theta * max_{k != i} |a_ik|
    """
    S = coo_matrix(A)

    # remove diagonals
    mask = S.row != S.col
    S.row = S.row[mask]
    S.col = S.col[mask]
    S.data = S.data[mask]
    max_offdiag = np.empty(S.shape[0])
    max_offdiag[:] = np.finfo(S.data.dtype).min

    # Note abs(.) takes the complex modulus
    if norm == 'abs':
        for i, v in zip(S.row, S.data):
            max_offdiag[i] = max(max_offdiag[i], abs(v))
    if norm == 'min':
        for i, v in zip(S.row, S.data):
            max_offdiag[i] = max(max_offdiag[i], -v)

    # strong connections
    if norm == 'abs':
        mask = np.abs(S.data) >= (theta * max_offdiag[S.row])
    if norm == 'min':
        mask = -S.data >= (theta * max_offdiag[S.row])

    S.row = S.row[mask]
    S.col = S.col[mask]
    S.data = S.data[mask]

    # Add back diagonal
    D = scipy.sparse.eye(S.shape[0], S.shape[0], format="csr", dtype=A.dtype)
    D.data[:] = csr_matrix(A).diagonal()
    S = S.tocsr() + D

    # Strength represents "distance", so take the magnitude
    S.data = np.abs(S.data)

    # Scale S by the largest magnitude entry in each row
    largest_row_entry = np.zeros((S.shape[0],), dtype=S.dtype)
    for i in range(S.shape[0]):
        for j in range(S.indptr[i], S.indptr[i+1]):
            val = abs(S.data[j])
            if val > largest_row_entry[i]:
                largest_row_entry[i] = val

    largest_row_entry[largest_row_entry != 0] =\
        1.0 / largest_row_entry[largest_row_entry != 0]
    S = S.tocsr()
    S = scale_rows(S, largest_row_entry, copy=True)

    return S

def reference_symmetric_soc(A, theta):
    # This is just a direct complex extension of the classic
    # SA strength-of-connection measure.  The extension continues
    # to compare magnitudes. This should reduce to the classic
    # measure if A is all real.

    # if theta == 0:
    #    return A

    D = np.abs(A.diagonal())

    S = coo_matrix(A)

    mask = S.row != S.col
    DD = np.array(D[S.row] * D[S.col]).reshape(-1,)
    # Note that abs takes the complex modulus element-wise
    # Note that using the square of the measure is the technique used
    # in the C++ routine, so we use it here.  Doing otherwise causes errors.
    mask &= ((real(S.data)**2 + imag(S.data)**2) >= theta*theta*DD)

    S.row = S.row[mask]
    S.col = S.col[mask]
    S.data = S.data[mask]

    # Add back diagonal
    D = scipy.sparse.eye(S.shape[0], S.shape[0], format="csr", dtype=A.dtype)
    D.data[:] = csr_matrix(A).diagonal()
    S = S.tocsr() + D

    # Strength represents "distance", so take the magnitude
    S.data = np.abs(S.data)

    # Scale S by the largest magnitude entry in each row
    largest_row_entry = np.zeros((S.shape[0],), dtype=S.dtype)
    for i in range(S.shape[0]):
        for j in range(S.indptr[i], S.indptr[i+1]):
            val = abs(S.data[j])
            if val > largest_row_entry[i]:
                largest_row_entry[i] = val

    largest_row_entry[largest_row_entry != 0] =\
        1.0 / largest_row_entry[largest_row_entry != 0]
    S = S.tocsr()
    S = scale_rows(S, largest_row_entry, copy=True)

    return S

def my_vis(ml, V, error=None, fname="", E2V=None, Pcols=None):
    """Coarse grid visualization for 2-D problems, for use with Paraview
       For all levels, outputs meshes, aggregates, near nullspace modes B, and selected
       prolongator basis functions.  Coarse level meshes are constructed by doing a
       Delaunay triangulation of interpolated fine grid vertices.

    Parameters
    ----------
    ml : {multilevel hiearchy}
        defines the multilevel hierarchy to visualize
    V : {array}
        coordinate array (N x D)
    Error : {array}
        Fine grid error to plot (N x D)
    fname : {string}
        string to be appended to all output files, e.g. 'diffusion1'
    E2V : {array}
        Element index array (Nel x Nelnodes) for the finest level.  If None,
        then a Delaunay triangulation is done for the finest level.  All coarse
        levels use an internally calculated Delaunay triangulation
    P_cols : {list of tuples}
        Optional input list of tuples of the form [(lvl, [ints]), ...]
        where lvl is an integer defining the level on which to output
        the list of columns in [ints].

    Returns
    -------
        - Writes data to .vtk files for use in paraview (xml 0.1 format)
    
    Notes
    -----


    Examples
    --------

     """
    system('rm -f *.vtu')

    ##
    # For the purposes of clearer plotting, perturb vertices slightly
    V += rand(V.shape[0], V.shape[1])*1e-6

    ## 
    # Create a list of vertices and meshes for all levels
    levels = ml.levels
    Vlist = [V]
    if E2V is None:
        [circ_cent,edges,E2V,tri_nbs]=delaunay.delaunay(V[:,0], V[:,1])
    E2Vlist = [E2V]

    mesh_type_list = []
    mesh_num_list = []
    if E2V.shape[1] == 1:
        mesh_type_list.append('vertex')
        mesh_num_list.append(1)
    if E2V.shape[1] == 3:
        mesh_type_list.append('tri')
        mesh_num_list.append(5)
    if E2V.shape[1] == 4:
        if vertices.shape[1] == 2:
            mesh_type_list.append('quad')
            mesh_num_list.append(9)
    
    if sparse.isspmatrix_bsr(levels[0].A):
        nPDEs = levels[0].A.blocksize[0]
    else:
        nPDEs = 1
    
    Agglist = []
    Agg = sparse.eye(levels[0].A.shape[0]/nPDEs, levels[0].A.shape[1]/nPDEs, format='csr') 
    for i in range(1,len(levels)):
        ##
        # Interpolate the vertices to the next level by taking each
        # aggregate's center of gravity (i.e. average x and y value).
        Agg = Agg.tocsr()*levels[i-1].AggOp.tocsr()
        Agg.data[:] = 1.0
        Agglist.append(Agg)
            
        AggX = scale_rows(Agg, Vlist[0][:,0], copy=True) 
        AggY = scale_rows(Agg, Vlist[0][:,1], copy=True) 
        AggX = ones((1, AggX.shape[0]))*AggX
        AggY = ones((1, AggY.shape[0]))*AggY
        Agg = Agg.tocsc()
        count = Agg.indptr[1:]-Agg.indptr[:-1]
        AggX = (ravel(AggX)/count).reshape(-1,1)
        AggY = (ravel(AggY)/count).reshape(-1,1)
        Vlist.append(hstack((AggX, AggY)))

        [circ_cent,edges,E2Vnew,tri_nbs]=delaunay.delaunay(Vlist[i][:,0], Vlist[i][:,1])
        E2Vlist.append(E2Vnew)
        mesh_type_list.append('tri')
        mesh_num_list.append(5)

        
    ##
    # On each level, output aggregates, B, the mesh
    for i in range(len(levels)):
        mesh_num = mesh_num_list[i]
        mesh_type = mesh_type_list[i]
#.........这里部分代码省略.........

#.........这里部分代码省略.........
                                                           #   converted to upper tri with Givens Rots 


    # GMRES will be run with diagonal preconditioning
    if weighting == 'diagonal':
        Dinv = get_diagonal(A, norm_eq=False, inv=True)
    elif weighting == 'block':
        Dinv = get_block_diag(A, blocksize=A.blocksize[0], inv_flag=True)
        Dinv = bsr_matrix( (Dinv, numpy.arange(Dinv.shape[0]), numpy.arange(Dinv.shape[0]+1)), shape = A.shape)
    elif weighting == 'local':
        # Based on Gershgorin estimate
        D = numpy.abs(A)*numpy.ones((A.shape[0],1), dtype=A.dtype)
        Dinv = numpy.zeros_like(D)
        Dinv[D != 0] = 1.0 / numpy.abs(D[D != 0])
    else:
        raise ValueError('weighting value is invalid')

    # Calculate initial residual
    #   Equivalent to R = -A*T;    R = R.multiply(Sparsity_Pattern)
    #   with the added constraint that R has an explicit 0 wherever 
    #   R is 0 and Sparsity_Pattern is not
    R = bsr_matrix((numpy.zeros(Sparsity_Pattern.data.shape, dtype=T.dtype), 
                    Sparsity_Pattern.indices, Sparsity_Pattern.indptr), 
                    shape=(Sparsity_Pattern.shape) )
    pyamg.amg_core.incomplete_mat_mult_bsr(A.indptr,  A.indices,  numpy.ravel(A.data), 
                                           T.indptr,  T.indices,  numpy.ravel(T.data),
                                           R.indptr,  R.indices,  numpy.ravel(R.data),
                                           T.shape[0]/T.blocksize[0], T.shape[1]/T.blocksize[1],
                                           A.blocksize[0], A.blocksize[1], T.blocksize[1])
    R.data *= -1.0
 
    #Apply diagonal preconditioner
    if weighting == 'local' or weighting == 'diagonal':
        R = scale_rows(R, Dinv)
    else:
        R = Dinv*R
    
    # Enforce R*B = 0
    Satisfy_Constraints(R, B, BtBinv)

    if R.nnz == 0:
        print "Error in sa_energy_min(..).  Initial R no nonzeros on a level.  Returning tentative prolongator\n"
        return T
    
    # This is the RHS vector for the problem in the Krylov Space
    normr = numpy.sqrt((R.data.conjugate()*R.data).sum())
    g = numpy.zeros((maxiter+1,), dtype=xtype) 
    g[0] = normr
    
    # First Krylov vector
    # V[0] = r/normr
    if normr > 0.0:
        V.append((1.0/normr)*R)

    #print "Energy Minimization of Prolongator --- Iteration 0 --- r = " + str(normr)
    i = -1
    #vect = numpy.ravel((A*T).data)
    #print "Iteration " + str(i+1) + "   Energy = %1.3e"%numpy.sqrt( (vect.conjugate()*vect).sum() )
    #print "Iteration " + str(i+1) + "   Normr  %1.3e"%normr
    while i < maxiter-1 and normr > tol:
        i = i+1

        # Calculate new search direction
        #   Equivalent to:  AV = A*V;    AV = AV.multiply(Sparsity_Pattern)
        #   with the added constraint that explicit zeros are in AP wherever 
        #   AP = 0 and Sparsity_Pattern does not

def jacobi_prolongation_smoother(S, T, C, B, omega=4.0/3.0, degree=1, filter=False, weighting='diagonal'):
    """Jacobi prolongation smoother
   
    Parameters
    ----------
    S : {csr_matrix, bsr_matrix}
        Sparse NxN matrix used for smoothing.  Typically, A.
    T : {csr_matrix, bsr_matrix}
        Tentative prolongator
    C : {csr_matrix, bsr_matrix}
        Strength-of-connection matrix
    B : {array}
        Near nullspace modes for the coarse grid such that T*B 
        exactly reproduces the fine grid near nullspace modes
    omega : {scalar}
        Damping parameter
    filter : {boolean}
        If true, filter S before smoothing T.  This option can greatly control
        complexity.
    weighting : {string}
        'block', 'diagonal' or 'local' weighting for constructing the Jacobi D
        'local': Uses a local row-wise weight based on the Gershgorin estimate.
          Avoids any potential under-damping due to inaccurate spectral radius
          estimates.
        'block': If A is a BSR matrix, use a block diagonal inverse of A  
        'diagonal': Classic Jacobi D = diagonal(A)

    Returns
    -------
    P : {csr_matrix, bsr_matrix}
        Smoothed (final) prolongator defined by P = (I - omega/rho(K) K) * T
        where K = diag(S)^-1 * S and rho(K) is an approximation to the 
        spectral radius of K.

    Notes
    -----
    If weighting is not 'local', then results using Jacobi prolongation
    smoother are not precisely reproducible due to a random initial guess used
    for the spectral radius approximation.  For precise reproducibility, 
    set numpy.random.seed(..) to the same value before each test. 
    
    Examples
    --------
    >>> from pyamg.aggregation import jacobi_prolongation_smoother
    >>> from pyamg.gallery import poisson
    >>> from scipy.sparse import coo_matrix
    >>> import numpy
    >>> data = numpy.ones((6,))
    >>> row = numpy.arange(0,6)
    >>> col = numpy.kron([0,1],numpy.ones((3,)))
    >>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
    >>> T.todense()
    matrix([[ 1.,  0.],
            [ 1.,  0.],
            [ 1.,  0.],
            [ 0.,  1.],
            [ 0.,  1.],
            [ 0.,  1.]])
    >>> A = poisson((6,),format='csr')
    >>> P = jacobi_prolongation_smoother(A,T,A,numpy.ones((2,1)))
    >>> P.todense()
    matrix([[ 0.64930164,  0.        ],
            [ 1.        ,  0.        ],
            [ 0.64930164,  0.35069836],
            [ 0.35069836,  0.64930164],
            [ 0.        ,  1.        ],
            [ 0.        ,  0.64930164]])

    """

    # preprocess weighting
    if weighting == 'block':
        if isspmatrix_csr(S):
            weighting = 'diagonal'
        elif isspmatrix_bsr(S):
            if S.blocksize[0] == 1:
                weighting = 'diagonal'
    
    if filter:
        ##
        # Implement filtered prolongation smoothing for the general case by
        # utilizing satisfy constraints

        if isspmatrix_bsr(S):
            numPDEs = S.blocksize[0]
        else:
            numPDEs = 1

        # Create a filtered S with entries dropped that aren't in C
        C = UnAmal(C, numPDEs, numPDEs)
        S = S.multiply(C)
        S.eliminate_zeros()

    if weighting == 'diagonal':
        # Use diagonal of S
        D_inv = get_diagonal(S, inv=True)
        D_inv_S = scale_rows(S, D_inv, copy=True)
        D_inv_S = (omega/approximate_spectral_radius(D_inv_S))*D_inv_S
    elif weighting == 'block':
        # Use block diagonal of S
#.........这里部分代码省略.........

#.........这里部分代码省略.........
    Ah = A.H
    Ah.sort_indices()
    
    # Preallocate
    AP = bsr_matrix((numpy.zeros(Sparsity_Pattern.data.shape, dtype=T.dtype), 
                     Sparsity_Pattern.indices, Sparsity_Pattern.indptr), 
                     shape=(Sparsity_Pattern.shape) )
 
    # D for A.H*A
    Dinv = get_diagonal(A, norm_eq=1, inv=True)
 
    # Calculate initial residual
    #   Equivalent to R = -Ah*(A*T);    R = R.multiply(Sparsity_Pattern)
    #   with the added constraint that R has an explicit 0 wherever 
    #   R is 0 and Sparsity_Pattern is not
    R = bsr_matrix((numpy.zeros(Sparsity_Pattern.data.shape, dtype=T.dtype), 
                    Sparsity_Pattern.indices, Sparsity_Pattern.indptr), 
                    shape=(Sparsity_Pattern.shape) )
    AT = -1.0*A*T
    R.data[:] = 0.0
    pyamg.amg_core.incomplete_mat_mult_bsr(Ah.indptr, Ah.indices, numpy.ravel(Ah.data), 
                                           AT.indptr, AT.indices, numpy.ravel(AT.data),
                                           R.indptr,  R.indices,  numpy.ravel(R.data),
                                           T.shape[0]/T.blocksize[0], T.shape[1]/T.blocksize[1],
                                           Ah.blocksize[0], Ah.blocksize[1], T.blocksize[1])

    # Enforce R*B = 0
    Satisfy_Constraints(R, B, BtBinv)
 
    if R.nnz == 0:
        print "Error in sa_energy_min(..).  Initial R no nonzeros on a level.  Returning tentative prolongator\n"
        return T
    
    #Calculate Frobenius norm of the residual
    resid = R.nnz #numpy.sqrt((R.data.conjugate()*R.data).sum())
    #print "Energy Minimization of Prolongator --- Iteration 0 --- r = " + str(resid)

    i = 0
    while i < maxiter and resid > tol:
        
        vect = numpy.ravel((A*T).data)
        #print "Iteration " + str(i) + "   Energy = %1.3e"%numpy.sqrt( (vect.conjugate()*vect).sum() )

        #Apply diagonal preconditioner
        Z = scale_rows(R, Dinv)
 
        #Frobenius innerproduct of (R,Z) = sum(rk.*zk)
        newsum = (R.conjugate().multiply(Z)).sum()
        if newsum < tol:
            # met tolerance, so halt
            break
            
        #P is the search direction, not the prolongator, which is T.    
        if(i == 0):
            P = Z
        else:
            beta = newsum/oldsum
            P = Z + beta*P
        oldsum = newsum
 
        #Calculate new direction
        #  Equivalent to:  AP = Ah*(A*P);    AP = AP.multiply(Sparsity_Pattern)
        #  with the added constraint that explicit zeros are in AP wherever 
        #  AP = 0 and Sparsity_Pattern does not
        AP_temp = A*P
        AP.data[:] = 0.0
        pyamg.amg_core.incomplete_mat_mult_bsr(Ah.indptr,      Ah.indices,      numpy.ravel(Ah.data), 
                                               AP_temp.indptr, AP_temp.indices, numpy.ravel(AP_temp.data),
                                               AP.indptr,      AP.indices,      numpy.ravel(AP.data),
                                               T.shape[0]/T.blocksize[0], T.shape[1]/T.blocksize[1],
                                               Ah.blocksize[0], Ah.blocksize[1], T.blocksize[1])
        del AP_temp 


        # Enforce AP*B = 0
        Satisfy_Constraints(AP, B, BtBinv)
        
        #Frobenius inner-product of (P, AP)
        alpha = newsum/(P.conjugate().multiply(AP)).sum()
 
        #Update the prolongator, T
        T = T + alpha*P 
 
        # Ensure identity at C-pts 
        if Cpt_params[0]: 
            T = Cpt_params[1]['I_F']*T + Cpt_params[1]['P_I']

        #Update residual
        R = R - alpha*AP
        
        i += 1

        #Calculate Frobenius norm of the residual
        resid = R.nnz #numpy.sqrt((R.data.conjugate()*R.data).sum())
        #print "Energy Minimization of Prolongator --- Iteration " + str(i) + " --- r = " + str(resid)

    vect = numpy.ravel((A*T).data)
    #print "Final Iteration " + str(i) + "   Energy = %1.3e"%numpy.sqrt( (vect.conjugate()*vect).sum() )

    return T

#.........这里部分代码省略.........
    '''

    # Preallocate
    AP = bsr_matrix((numpy.zeros(Sparsity_Pattern.data.shape, dtype=T.dtype), Sparsity_Pattern.indices, Sparsity_Pattern.indptr), 
                     shape=(Sparsity_Pattern.shape) )

    # CG will be run with diagonal preconditioning
    if weighting == 'diagonal':
        Dinv = get_diagonal(A, norm_eq=False, inv=True)
    elif weighting == 'block':
        Dinv = get_block_diag(A, blocksize=A.blocksize[0], inv_flag=True)
        Dinv = bsr_matrix( (Dinv, numpy.arange(Dinv.shape[0]), numpy.arange(Dinv.shape[0]+1)), shape = A.shape)
    elif weighting == 'local':
        # Based on Gershgorin estimate
        D = numpy.abs(A)*numpy.ones((A.shape[0],1), dtype=A.dtype)
        Dinv = numpy.zeros_like(D)
        Dinv[D != 0] = 1.0 / numpy.abs(D[D != 0])
    else:
        raise ValueError('weighting value is invalid')

    # Calculate initial residual
    #   Equivalent to R = -A*T;    R = R.multiply(Sparsity_Pattern)
    #   with the added constraint that R has an explicit 0 wherever 
    #   R is 0 and Sparsity_Pattern is not
    R = bsr_matrix((numpy.zeros(Sparsity_Pattern.data.shape, dtype=T.dtype), Sparsity_Pattern.indices, 
                   Sparsity_Pattern.indptr), shape=(Sparsity_Pattern.shape) )
    pyamg.amg_core.incomplete_mat_mult_bsr(A.indptr,  A.indices,  numpy.ravel(A.data), 
                                           T.indptr,  T.indices,  numpy.ravel(T.data),
                                           R.indptr, R.indices,   numpy.ravel(R.data),
                                           T.shape[0]/T.blocksize[0], T.shape[1]/T.blocksize[1],
                                           A.blocksize[0], A.blocksize[1], T.blocksize[1])
    R.data *= -1.0

    # Enforce R*B = 0
    Satisfy_Constraints(R, B, BtBinv)

    if R.nnz == 0:
        print "Error in sa_energy_min(..).  Initial R no nonzeros on a level.  Returning tentative prolongator\n"
        return T
    
    #Calculate Frobenius norm of the residual
    resid = R.nnz ##numpy.sqrt((R.data.conjugate()*R.data).sum())
    #print "Energy Minimization of Prolongator --- Iteration 0 --- r = " + str(resid)
    
    i = 0
    while i < maxiter and resid > tol:
        #Apply diagonal preconditioner
        if weighting == 'local' or weighting == 'diagonal':
            Z = scale_rows(R, Dinv)
        else:
            Z = Dinv*R

        #Frobenius inner-product of (R,Z) = sum( numpy.conjugate(rk).*zk)
        newsum = (R.conjugate().multiply(Z)).sum()
        if newsum < tol:
            # met tolerance, so halt
            break

        #P is the search direction, not the prolongator, which is T.    
        if(i == 0):
            P = Z
        else:
            beta = newsum/oldsum
            P = Z + beta*P
        oldsum = newsum
        
        # Calculate new direction and enforce constraints
        #   Equivalent to:  AP = A*P;    AP = AP.multiply(Sparsity_Pattern)
        #   with the added constraint that explicit zeros are in AP wherever 
        #   AP = 0 and Sparsity_Pattern does not  !!!!
        AP.data[:] = 0.0
        pyamg.amg_core.incomplete_mat_mult_bsr(A.indptr,  A.indices,  numpy.ravel(A.data), 
                                               P.indptr,  P.indices,  numpy.ravel(P.data),
                                               AP.indptr, AP.indices, numpy.ravel(AP.data),
                                               T.shape[0]/T.blocksize[0], T.shape[1]/T.blocksize[1],
                                               A.blocksize[0], A.blocksize[1], P.blocksize[1])

        # Enforce AP*B = 0
        Satisfy_Constraints(AP, B, BtBinv)
        
        #Frobenius inner-product of (P, AP)
        alpha = newsum/(P.conjugate().multiply(AP)).sum()

        #Update the prolongator, T
        T = T + alpha*P 

        # Ensure identity at C-pts 
        if Cpt_params[0]:
            T = Cpt_params[1]['I_F']*T + Cpt_params[1]['P_I']

        #Update residual
        R = R - alpha*AP
        
        i += 1

        #Calculate Frobenius norm of the residual
        resid = R.nnz #numpy.sqrt((R.data.conjugate()*R.data).sum())
        #print "Energy Minimization of Prolongator --- Iteration " + str(i) + " --- r = " + str(resid)
   
    return T

def evolution_strength_of_connection(A, B='ones', epsilon=4.0, k=2,
                                     proj_type="l2", block_flag=False,
                                     symmetrize_measure=True):
    """
    Construct strength of connection matrix using an Evolution-based measure

    Parameters
    ----------
    A : {csr_matrix, bsr_matrix}
        Sparse NxN matrix
    B : {string, array}
        If B='ones', 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

    References
    ----------
    .. [1] 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)))
    """
    # local imports for evolution_strength_of_connection
    from pyamg.util.utils import scale_rows, get_block_diag, scale_columns
    from pyamg.util.linalg import approximate_spectral_radius

    # ====================================================================
    # 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 == 'ones':
        Bmat = np.mat(np.ones((A.shape[0], 1), dtype=A.dtype))
    else:
        Bmat = np.mat(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()
#.........这里部分代码省略.........

def reference_distance_soc(A, V, theta=2.0, relative_drop=True):
    '''
    Reference routine for distance based strength of connection
    '''

    # deal with the supernode case
    if isspmatrix_bsr(A):
        dimen = int(A.shape[0]/A.blocksize[0])
        C = csr_matrix((np.ones((A.data.shape[0],)), A.indices, A.indptr),
                       shape=(dimen, dimen))
    else:
        A = A.tocsr()
        dimen = A.shape[0]
        C = A.copy()
        C.data = np.real(C.data)

    if V.shape[1] == 2:
        three_d = False
    elif V.shape[1] == 3:
        three_d = True

    for i in range(dimen):
        rowstart = C.indptr[i]
        rowend = C.indptr[i+1]
        pt_i = V[i, :]
        for j in range(rowstart, rowend):
            if C.indices[j] == i:
                # ignore the diagonal entry by making it large
                C.data[j] = np.finfo(np.float).max
            else:
                # distance between entry j and i
                pt_j = V[C.indices[j], :]
                dist = (pt_i[0] - pt_j[0])**2
                dist += (pt_i[1] - pt_j[1])**2
                if three_d:
                    dist += (pt_i[2] - pt_j[2])**2
                C.data[j] = np.sqrt(dist)

        # apply drop tolerance
        this_row = C.data[rowstart:rowend]
        if relative_drop:
            tol_i = theta*this_row.min()
            this_row[this_row > tol_i] = 0.0
        else:
            this_row[this_row > theta] = 0.0

        C.data[rowstart:rowend] = this_row

    C.eliminate_zeros()
    C = C + 2.0*scipy.sparse.eye(C.shape[0], C.shape[1], format='csr')

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

    # Scale C by the largest magnitude entry in each row
    largest_row_entry = np.zeros((C.shape[0],), dtype=C.dtype)
    for i in range(C.shape[0]):
        for j in range(C.indptr[i], C.indptr[i+1]):
            val = abs(C.data[j])
            if val > largest_row_entry[i]:
                largest_row_entry[i] = val

    largest_row_entry[largest_row_entry != 0] =\
        1.0 / largest_row_entry[largest_row_entry != 0]
    C = C.tocsr()
    C = scale_rows(C, largest_row_entry, copy=True)

    return C

def reference_evolution_soc(A, B, epsilon=4.0, k=2, proj_type="l2"):
    """
    All python reference implementation for Evolution Strength of Connection

    --> If doing imaginary test, both A and B should be imaginary type upon
    entry

    --> This does the "unsymmetrized" version of the ode measure
    """

    # number of PDEs per point is defined implicitly by block size
    csrflag = isspmatrix_csr(A)
    if csrflag:
        numPDEs = 1
    else:
        numPDEs = A.blocksize[0]
        A = A.tocsr()

    # Preliminaries
    near_zero = np.finfo(float).eps
    sqrt_near_zero = np.sqrt(np.sqrt(near_zero))
    Bmat = np.mat(B)
    A.eliminate_zeros()
    A.sort_indices()
    dimen = A.shape[1]
    NullDim = Bmat.shape[1]

    # Get spectral radius of Dinv*A, this is the time step size for the ODE
    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)
    rho_DinvA = approximate_spectral_radius(Dinv_A)

    # Calculate (Atilde^k) naively
    S = (scipy.sparse.eye(dimen, dimen, format="csr") - (1.0/rho_DinvA)*Dinv_A)
    Atilde = scipy.sparse.eye(dimen, dimen, format="csr")
    for i in range(k):
        Atilde = S*Atilde

    # Strength Info should be row-based, so transpose Atilde
    Atilde = Atilde.T.tocsr()

    # Construct and apply a sparsity mask for Atilde that restricts 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()

    # Only consider strength at dofs from your PDE.  Use mask to enforce this
    # by zeroing out all entries in Atilde that aren't from your 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()

    # 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()
    del mask

    # Calculate strength based on constrained min problem of
    LHS = np.mat(np.zeros((NullDim+1, NullDim+1)), dtype=A.dtype)
    RHS = np.mat(np.zeros((NullDim+1, 1)), dtype=A.dtype)

    # Choose tolerance for dropping "numerically zero" values later
    t = Atilde.dtype.char
    eps = np.finfo(np.float).eps
    feps = np.finfo(np.single).eps
    geps = np.finfo(np.longfloat).eps
    _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G': 2}
    tol = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]]

    for i in range(dimen):

        # Get rowptrs and col indices from Atilde
        rowstart = Atilde.indptr[i]
        rowend = Atilde.indptr[i+1]
        length = rowend - rowstart
        colindx = Atilde.indices[rowstart:rowend]

        # Local diagonal of A is used for scale invariant min problem
        D_A = np.mat(np.eye(length, dtype=A.dtype))
        if proj_type == "D_A":
            for j in range(length):
                D_A[j, j] = D[colindx[j]]

        # Find row i's position in colindx, matrix must have sorted column
        # indices.
        iInRow = colindx.searchsorted(i)

        if length <= NullDim:
#.........这里部分代码省略.........