(A, T, B, BtBinv, pattern, maxiter,
                                tol, weighting='diagonal', Cpt_params=None)
return T
Smooth T with CGNR by solving A T = 0, subject to nullspace and sparsity constraints.

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        SPD sparse NxN matrix
        Should be at least nonsymmetric or indefinite
    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
        IGNORED here, only 'diagonal' preconditioning is used.
    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 CGNR 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
Smooth T with CGNR by solving A T = 0, subject to nullspace and sparsity constraints.
Smooth T with CGNR by solving A T = 0 subject to nullspace and sparsity constraints .
def cgnr_prolongation_smoothing(A, T, B, BtBinv, pattern, maxiter,
                                tol, weighting='diagonal', Cpt_params=None):
    """Smooth T with CGNR by solving A T = 0, subject to nullspace and sparsity constraints.

    Parameters
    ----------
    A : csr_matrix, bsr_matrix
        SPD sparse NxN matrix
        Should be at least nonsymmetric or indefinite
    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
        IGNORED here, only 'diagonal' preconditioning is used.
    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 CGNR 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

    """
    if weighting != 'diagonal':
        warn(f'Weighting of {weighting} unused.', stacklevel=2)

    # For non-SPD system, apply CG on Normal Equations with Diagonal
    # Preconditioning (requires transpose)
    Ah = A.H
    Ah.sort_indices()

    # Preallocate
    uones = np.zeros(pattern.data.shape, dtype=T.dtype)
    AP = sparse.bsr_matrix((uones, pattern.indices, pattern.indptr),
                           shape=(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(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))
    AT = -1.0*A*T
    R.data[:] = 0.0
    amg_core.incomplete_mat_mult_bsr(Ah.indptr, Ah.indices,
                                     np.ravel(Ah.data),
                                     AT.indptr, AT.indices,
                                     np.ravel(AT.data),
                                     R.indptr, R.indices,
                                     np.ravel(R.data),
                                     int(T.shape[0]/T.blocksize[0]),
                                     int(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')
        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:

        # vect = np.ravel((A*T).data)
        # print "Iteration " + str(i) + "   \
        # Energy = %1.3e"%np.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
            oldsum = newsum
        else:
            beta = newsum/oldsum
            P = Z + beta*P
        oldsum = newsum

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

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

    return T
def cgnr_prolongation_smoothing ( A , T , B , BtBinv , pattern , maxiter , tol , weighting = 'diagonal' , Cpt_params = None ) : if weighting != 'diagonal' : warn ( f'Weighting of {weighting} unused.' , stacklevel = 2 ) # For non-SPD system, apply CG on Normal Equations with Diagonal # Preconditioning (requires transpose) Ah = A . H Ah . sort_indices ( ) # Preallocate uones = np . zeros ( pattern . data . shape , dtype = T . dtype ) AP = sparse . bsr_matrix ( ( uones , pattern . indices , pattern . indptr ) , shape = ( 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(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 ) ) AT = - 1.0 * A * T R . data [ : ] = 0.0 amg_core . incomplete_mat_mult_bsr ( Ah . indptr , Ah . indices , np . ravel ( Ah . data ) , AT . indptr , AT . indices , np . ravel ( AT . data ) , R . indptr , R . indices , np . ravel ( R . data ) , int ( T . shape [ 0 ] / T . blocksize [ 0 ] ) , int ( 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' ) 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 : # vect = np.ravel((A*T).data) # print "Iteration " + str(i) + "   \ # Energy = %1.3e"%np.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 oldsum = newsum else : beta = newsum / oldsum P = Z + beta * P oldsum = newsum # Calculate new direction #  Equivalent to:  AP = Ah*(A*P);    AP = AP.multiply(pattern) #  with the added constraint that explicit zeros are in AP wherever #  AP = 0 and pattern does not AP_temp = A * P AP . data [ : ] = 0.0 amg_core . incomplete_mat_mult_bsr ( Ah . indptr , Ah . indices , np . ravel ( Ah . data ) , AP_temp . indptr , AP_temp . indices , np . ravel ( AP_temp . data ) , AP . indptr , AP . indices , np . ravel ( AP . data ) , int ( T . shape [ 0 ] / T . blocksize [ 0 ] ) , int ( 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 # np.sqrt((R.data.conjugate()*R.data).sum()) # print "Energy Minimization of Prolongator \ # --- Iteration " + str(i) + " --- r = " + str(resid) # vect = np.ravel((A*T).data) # print "Final Iteration " + str(i) + "   \ # Energy = %1.3e"%np.sqrt( (vect.conjugate()*vect).sum() ) return T