(A, tol=1e-5, maxiter=10)
return (1/scale)*C
Binormalize matrix A.  Attempt to create unit l_1 norm rows.

    Parameters
    ----------
    A : csr_matrix
        sparse matrix (n x n)
    tol : float
        tolerance
    x : array
        guess at the diagonal
    maxiter : int
        maximum number of iterations to try

    Returns
    -------
    C : csr_matrix
        diagonally scaled A, C=DAD

    Notes
    -----
        - Goal: Scale A so that l_1 norm of the rows are equal to 1:
        - B = DAD
        - want row sum of B = 1
        - easily done with tol=0 if B=DA, but this is not symmetric
        - algorithm is O(N log (1.0/tol))

    Examples
    --------
    >>> from pyamg.gallery import poisson
    >>> from pyamg.classical.cr import binormalize
    >>> A = poisson((10,),format='csr')
    >>> C = binormalize(A)

    References
    ----------
    .. [1] Livne, Golub, "Scaling by Binormalization"
       Tech Report SCCM-03-12, SCCM, Stanford, 2003
       http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
Binormalize matrix A.  Attempt to create unit l_1 norm rows.
Binormalize matrix A . Attempt to create unit l_1 norm rows .
def binormalize(A, tol=1e-5, maxiter=10):
    """Binormalize matrix A.  Attempt to create unit l_1 norm rows.

    Parameters
    ----------
    A : csr_matrix
        sparse matrix (n x n)
    tol : float
        tolerance
    x : array
        guess at the diagonal
    maxiter : int
        maximum number of iterations to try

    Returns
    -------
    C : csr_matrix
        diagonally scaled A, C=DAD

    Notes
    -----
        - Goal: Scale A so that l_1 norm of the rows are equal to 1:
        - B = DAD
        - want row sum of B = 1
        - easily done with tol=0 if B=DA, but this is not symmetric
        - algorithm is O(N log (1.0/tol))

    Examples
    --------
    >>> from pyamg.gallery import poisson
    >>> from pyamg.classical.cr import binormalize
    >>> A = poisson((10,),format='csr')
    >>> C = binormalize(A)

    References
    ----------
    .. [1] Livne, Golub, "Scaling by Binormalization"
       Tech Report SCCM-03-12, SCCM, Stanford, 2003
       http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679

    """
    if not isspmatrix(A):
        raise TypeError('expecting sparse matrix A')

    if A.dtype == complex:
        raise NotImplementedError('complex A not implemented')

    n = A.shape[0]
    it = 0
    x = np.ones((n, 1)).ravel()

    # 1.
    B = A.multiply(A).tocsc()  # power(A,2) inconsistent in numpy, scipy.sparse
    d = B.diagonal().ravel()

    # 2.
    beta = B * x
    betabar = (1.0/n) * np.dot(x, beta)
    stdev = rowsum_stdev(x, beta)

    # 3
    while stdev > tol and it < maxiter:
        for i in range(0, n):
            # solve equation x_i, keeping x_j's fixed
            # see equation (12)
            c2 = (n-1)*d[i]
            c1 = (n-2)*(beta[i] - d[i]*x[i])
            c0 = -d[i]*x[i]*x[i] + 2*beta[i]*x[i] - n*betabar

            if -c0 < 1e-14:
                print('warning: A nearly un-binormalizable...')
                return A

            # see equation (12)
            xnew = (2*c0)/(-c1 - np.sqrt(c1*c1 - 4*c0*c2))
            dx = xnew - x[i]

            # here we assume input matrix is symmetric since we grab a row of B
            # instead of a column
            ii = B.indptr[i]
            iii = B.indptr[i+1]
            dot_Bcol = np.dot(x[B.indices[ii:iii]], B.data[ii:iii])

            betabar = betabar + (1.0/n)*dx*(dot_Bcol + beta[i] + d[i]*dx)
            beta[B.indices[ii:iii]] += dx*B.data[ii:iii]

            x[i] = xnew

        stdev = rowsum_stdev(x, beta)
        it += 1

    # rescale for unit 2-norm
    d = np.sqrt(x)
    D = spdiags(d.ravel(), [0], n, n)
    C = D * A * D
    C = C.tocsr()
    beta = C.multiply(C).sum(axis=1)
    scale = np.sqrt((1.0/n) * np.sum(beta))
    return (1/scale)*C
def binormalize ( A , tol = 1e-5 , maxiter = 10 ) : if not isspmatrix ( A ) : raise TypeError ( 'expecting sparse matrix A' ) if A . dtype == complex : raise NotImplementedError ( 'complex A not implemented' ) n = A . shape [ 0 ] it = 0 x = np . ones ( ( n , 1 ) ) . ravel ( ) # 1. B = A . multiply ( A ) . tocsc ( ) # power(A,2) inconsistent in numpy, scipy.sparse d = B . diagonal ( ) . ravel ( ) # 2. beta = B * x betabar = ( 1.0 / n ) * np . dot ( x , beta ) stdev = rowsum_stdev ( x , beta ) # 3 while stdev > tol and it < maxiter : for i in range ( 0 , n ) : # solve equation x_i, keeping x_j's fixed # see equation (12) c2 = ( n - 1 ) * d [ i ] c1 = ( n - 2 ) * ( beta [ i ] - d [ i ] * x [ i ] ) c0 = - d [ i ] * x [ i ] * x [ i ] + 2 * beta [ i ] * x [ i ] - n * betabar if - c0 < 1e-14 : print ( 'warning: A nearly un-binormalizable...' ) return A # see equation (12) xnew = ( 2 * c0 ) / ( - c1 - np . sqrt ( c1 * c1 - 4 * c0 * c2 ) ) dx = xnew - x [ i ] # here we assume input matrix is symmetric since we grab a row of B # instead of a column ii = B . indptr [ i ] iii = B . indptr [ i + 1 ] dot_Bcol = np . dot ( x [ B . indices [ ii : iii ] ] , B . data [ ii : iii ] ) betabar = betabar + ( 1.0 / n ) * dx * ( dot_Bcol + beta [ i ] + d [ i ] * dx ) beta [ B . indices [ ii : iii ] ] += dx * B . data [ ii : iii ] x [ i ] = xnew stdev = rowsum_stdev ( x , beta ) it += 1 # rescale for unit 2-norm d = np . sqrt ( x ) D = spdiags ( d . ravel ( ) , [ 0 ] , n , n ) C = D * A * D C = C . tocsr ( ) beta = C . multiply ( C ) . sum ( axis = 1 ) scale = np . sqrt ( ( 1.0 / n ) * np . sum ( beta ) ) return ( 1 / scale ) * C