(A, T, B, BtBinv, pattern, maxiter, tol,
                              weighting='local', Cpt_params=None)
return T
Use CG to smooth T by solving A T = 0, subject to nullspace and sparsity constraints.

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        SPD sparse NxN matrix
    T : bsr_matrix
        Tentative prolongator, a NxM sparse matrix (M < N).
        This is initial guess for the equation A T = 0.
        Assumed that T B_c = B_f
    B : array
        Near-nullspace modes for coarse grid, i.e., B_c.
        Has shape (M,k) where k is the number of coarse candidate vectors.
    BtBinv : array
        3 dimensional array such that,
        BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted
        to the neighborhood (in the matrix graph) of dof of i.
    pattern : csr_matrix, bsr_matrix
        Sparse NxM matrix
        This is the sparsity pattern constraint to enforce on the
        eventual prolongator
    maxiter : int
        maximum number of iterations
    tol : float
        residual tolerance for A T = 0
    weighting : string
        'block', 'diagonal' or 'local' construction of the diagonal
        preconditioning
    Cpt_params : tuple
        Tuple of the form (bool, dict).  If the Cpt_params[0] = False, then
        the standard SA prolongation smoothing is carried out.  If True, then
        dict must be a dictionary of parameters containing, (1) P_I: P_I.T is
        the injection matrix for the Cpts, (2) I_F: an identity matrix
        for only the F-points (i.e. I, but with zero rows and columns for
        C-points) and I_C: the C-point analogue to I_F.

    Returns
    -------
    T : bsr_matrix
        Smoothed prolongator using conjugate gradients to solve A T = 0,
        subject to the constraints, T B_c = B_f, and T has no nonzero
        outside of the sparsity pattern in pattern.

    See Also
    --------
    The principal calling routine,
    pyamg.aggregation.smooth.energy_prolongation_smoother
Use CG to smooth T by solving A T = 0, subject to nullspace and sparsity constraints.
Use CG to smooth T by solving A T = 0 subject to nullspace and sparsity constraints .
def cg_prolongation_smoothing(A, T, B, BtBinv, pattern, maxiter, tol,
                              weighting='local', Cpt_params=None):
    """Use CG to smooth T by solving A T = 0, subject to nullspace and sparsity constraints.

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        SPD sparse NxN matrix
    T : bsr_matrix
        Tentative prolongator, a NxM sparse matrix (M < N).
        This is initial guess for the equation A T = 0.
        Assumed that T B_c = B_f
    B : array
        Near-nullspace modes for coarse grid, i.e., B_c.
        Has shape (M,k) where k is the number of coarse candidate vectors.
    BtBinv : array
        3 dimensional array such that,
        BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted
        to the neighborhood (in the matrix graph) of dof of i.
    pattern : csr_matrix, bsr_matrix
        Sparse NxM matrix
        This is the sparsity pattern constraint to enforce on the
        eventual prolongator
    maxiter : int
        maximum number of iterations
    tol : float
        residual tolerance for A T = 0
    weighting : string
        'block', 'diagonal' or 'local' construction of the diagonal
        preconditioning
    Cpt_params : tuple
        Tuple of the form (bool, dict).  If the Cpt_params[0] = False, then
        the standard SA prolongation smoothing is carried out.  If True, then
        dict must be a dictionary of parameters containing, (1) P_I: P_I.T is
        the injection matrix for the Cpts, (2) I_F: an identity matrix
        for only the F-points (i.e. I, but with zero rows and columns for
        C-points) and I_C: the C-point analogue to I_F.

    Returns
    -------
    T : bsr_matrix
        Smoothed prolongator using conjugate gradients to solve A T = 0,
        subject to the constraints, T B_c = B_f, and T has no nonzero
        outside of the sparsity pattern in pattern.

    See Also
    --------
    The principal calling routine,
    pyamg.aggregation.smooth.energy_prolongation_smoother

    """
    # Preallocate
    AP = sparse.bsr_matrix((np.zeros(pattern.data.shape, dtype=T.dtype),
                            pattern.indices, pattern.indptr),
                           shape=(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 = sparse.bsr_matrix((Dinv, np.arange(Dinv.shape[0]),
                                  np.arange(Dinv.shape[0]+1)),
                                 shape=A.shape)
    elif weighting == 'local':
        # Based on Gershgorin estimate
        D = np.abs(A)*np.ones((A.shape[0], 1), dtype=A.dtype)
        Dinv = np.zeros_like(D)
        Dinv[D != 0] = 1.0 / np.abs(D[D != 0])
    else:
        raise ValueError('weighting value is invalid')

    # Calculate initial residual
    #   Equivalent to R = -A*T;    R = R.multiply(pattern)
    #   with the added constraint that R has an explicit 0 wherever
    #   R is 0 and pattern is not
    uones = np.zeros(pattern.data.shape, dtype=T.dtype)
    R = sparse.bsr_matrix((uones, pattern.indices,
                           pattern.indptr),
                          shape=(pattern.shape))
    amg_core.incomplete_mat_mult_bsr(A.indptr, A.indices,
                                     np.ravel(A.data),
                                     T.indptr, T.indices,
                                     np.ravel(T.data),
                                     R.indptr, R.indices,
                                     np.ravel(R.data),
                                     int(T.shape[0]/T.blocksize[0]),
                                     int(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  # np.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 in ('local', 'diagonal'):
            Z = scale_rows(R, Dinv)
        else:
            Z = Dinv*R

        # Frobenius inner-product of (R,Z) = sum( np.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
            oldsum = newsum
        else:
            beta = newsum / oldsum
            P = Z + beta*P
        oldsum = newsum

        # Calculate new direction and enforce constraints
        #   Equivalent to:  AP = A*P;    AP = AP.multiply(pattern)
        #   with the added constraint that explicit zeros are in AP wherever
        #   AP = 0 and pattern does not  !!!!
        AP.data[:] = 0.0
        amg_core.incomplete_mat_mult_bsr(A.indptr, A.indices,
                                         np.ravel(A.data),
                                         P.indptr, P.indices,
                                         np.ravel(P.data),
                                         AP.indptr, AP.indices,
                                         np.ravel(AP.data),
                                         int(T.shape[0]/T.blocksize[0]),
                                         int(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  # np.sqrt((R.data.conjugate()*R.data).sum())
        # print "Energy Minimization of Prolongator \
        # --- Iteration " + str(i) + " --- r = " + str(resid)

    return T
def cg_prolongation_smoothing ( A , T , B , BtBinv , pattern , maxiter , tol , weighting = 'local' , Cpt_params = None ) : # Preallocate AP = sparse . bsr_matrix ( ( np . zeros ( pattern . data . shape , dtype = T . dtype ) , pattern . indices , pattern . indptr ) , shape = ( 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 = sparse . bsr_matrix ( ( Dinv , np . arange ( Dinv . shape [ 0 ] ) , np . arange ( Dinv . shape [ 0 ] + 1 ) ) , shape = A . shape ) elif weighting == 'local' : # Based on Gershgorin estimate D = np . abs ( A ) * np . ones ( ( A . shape [ 0 ] , 1 ) , dtype = A . dtype ) Dinv = np . zeros_like ( D ) Dinv [ D != 0 ] = 1.0 / np . abs ( D [ D != 0 ] ) else : raise ValueError ( 'weighting value is invalid' ) # Calculate initial residual #   Equivalent to R = -A*T;    R = R.multiply(pattern) #   with the added constraint that R has an explicit 0 wherever #   R is 0 and pattern is not uones = np . zeros ( pattern . data . shape , dtype = T . dtype ) R = sparse . bsr_matrix ( ( uones , pattern . indices , pattern . indptr ) , shape = ( pattern . shape ) ) amg_core . incomplete_mat_mult_bsr ( A . indptr , A . indices , np . ravel ( A . data ) , T . indptr , T . indices , np . ravel ( T . data ) , R . indptr , R . indices , np . ravel ( R . data ) , int ( T . shape [ 0 ] / T . blocksize [ 0 ] ) , int ( 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 # np.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 in ( 'local' , 'diagonal' ) : Z = scale_rows ( R , Dinv ) else : Z = Dinv * R # Frobenius inner-product of (R,Z) = sum( np.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 oldsum = newsum else : beta = newsum / oldsum P = Z + beta * P oldsum = newsum # Calculate new direction and enforce constraints #   Equivalent to:  AP = A*P;    AP = AP.multiply(pattern) #   with the added constraint that explicit zeros are in AP wherever #   AP = 0 and pattern does not  !!!! AP . data [ : ] = 0.0 amg_core . incomplete_mat_mult_bsr ( A . indptr , A . indices , np . ravel ( A . data ) , P . indptr , P . indices , np . ravel ( P . data ) , AP . indptr , AP . indices , np . ravel ( AP . data ) , int ( T . shape [ 0 ] / T . blocksize [ 0 ] ) , int ( 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 # np.sqrt((R.data.conjugate()*R.data).sum()) # print "Energy Minimization of Prolongator \ # --- Iteration " + str(i) + " --- r = " + str(resid) return T