(A, T, B, BtBinv, pattern, maxiter,
                                 tol, weighting='local', Cpt_params=None)
return T
Smooth T with GMRES 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
    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 GMRES 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 GMRES by solving A T = 0 subject to nullspace and sparsity constraints.
Smooth T with GMRES by solving A T = 0 subject to nullspace and sparsity constraints .
def gmres_prolongation_smoothing(A, T, B, BtBinv, pattern, maxiter,
                                 tol, weighting='local', Cpt_params=None):
    """Smooth T with GMRES 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
    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 GMRES 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

    """
    # For non-SPD system, apply GMRES with Diagonal Preconditioning

    # Preallocate space for new search directions
    uones = np.zeros(pattern.data.shape, dtype=T.dtype)
    AV = sparse.bsr_matrix((uones, pattern.indices, pattern.indptr),
                           shape=(pattern.shape))

    # Preallocate for Givens Rotations, Hessenberg matrix and Krylov Space
    xtype = sparse.sputils.upcast(A.dtype, T.dtype, B.dtype)
    Q = []      # Givens Rotations
    V = []      # Krylov Space
    # vs = []     # vs store the pointers to each column of V for speed

    # Upper Hessenberg matrix, converted to upper tri with Givens Rots
    H = np.zeros((maxiter+1, maxiter+1), dtype=xtype)

    # 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 = 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

    # Apply diagonal preconditioner
    if weighting in ('local', '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 = np.sqrt((R.data.conjugate()*R.data).sum())
    g = np.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 = np.ravel((A*T).data)
    # print "Iteration " + str(i+1) + "   \
    # Energy = %1.3e"%np.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(pattern)
        #   with the added constraint that explicit zeros are in AP wherever
        #   AP = 0 and pattern does not
        AV.data[:] = 0.0
        amg_core.incomplete_mat_mult_bsr(A.indptr, A.indices,
                                         np.ravel(A.data),
                                         V[i].indptr, V[i].indices,
                                         np.ravel(V[i].data),
                                         AV.indptr, AV.indices,
                                         np.ravel(AV.data),
                                         int(T.shape[0]/T.blocksize[0]),
                                         int(T.shape[1]/T.blocksize[1]),
                                         A.blocksize[0], A.blocksize[1],
                                         T.blocksize[1])

        if weighting in ('local', 'diagonal'):
            AV = scale_rows(AV, Dinv)
        else:
            AV = Dinv*AV

        # Enforce AV*B = 0
        satisfy_constraints(AV, B, BtBinv)
        V.append(AV.copy())

        # Modified Gram-Schmidt
        for j in range(i+1):
            # Frobenius inner-product
            H[j, i] = (V[j].conjugate().multiply(V[i+1])).sum()
            V[i+1] = V[i+1] - H[j, i]*V[j]

        # Frobenius Norm
        H[i+1, i] = np.sqrt((V[i+1].data.conjugate()*V[i+1].data).sum())

        # Check for breakdown
        if H[i+1, i] != 0.0:
            V[i+1] = (1.0 / H[i+1, i]) * V[i+1]

        # Apply previous Givens rotations to H
        if i > 0:
            apply_givens(Q, H[:, i], i)

        # Calculate and apply next complex-valued Givens Rotation
        if H[i+1, i] != 0:
            h1 = H[i, i]
            h2 = H[i+1, i]
            h1_mag = np.abs(h1)
            h2_mag = np.abs(h2)
            if h1_mag < h2_mag:
                mu = h1/h2
                tau = np.conjugate(mu)/np.abs(mu)
            else:
                mu = h2/h1
                tau = mu/np.abs(mu)

            denom = np.sqrt(h1_mag**2 + h2_mag**2)
            c = h1_mag/denom
            s = h2_mag*tau/denom
            Qblock = np.array([[c, np.conjugate(s)], [-s, c]], dtype=xtype)
            Q.append(Qblock)

            # Apply Givens Rotation to g,
            #   the RHS for the linear system in the Krylov Subspace.
            g[i:i+2] = np.dot(Qblock, g[i:i+2])

            # Apply effect of Givens Rotation to H
            H[i, i] = np.dot(Qblock[0, :], H[i:i+2, i])
            H[i+1, i] = 0.0

        normr = np.abs(g[i+1])
        # print "Iteration " + str(i+1) + "   Normr  %1.3e"%normr
    # End while loop

    # Find best update to x in Krylov Space, V.  Solve (i x i) system.
    if i != -1:
        y = la.solve(H[0:i+1, 0:i+1], g[0:i+1])
        for j in range(i+1):
            T = T + y[j]*V[j]

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

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

    return T
def gmres_prolongation_smoothing ( A , T , B , BtBinv , pattern , maxiter , tol , weighting = 'local' , Cpt_params = None ) : # For non-SPD system, apply GMRES with Diagonal Preconditioning # Preallocate space for new search directions uones = np . zeros ( pattern . data . shape , dtype = T . dtype ) AV = sparse . bsr_matrix ( ( uones , pattern . indices , pattern . indptr ) , shape = ( pattern . shape ) ) # Preallocate for Givens Rotations, Hessenberg matrix and Krylov Space xtype = sparse . sputils . upcast ( A . dtype , T . dtype , B . dtype ) Q = [ ] # Givens Rotations V = [ ] # Krylov Space # vs = []     # vs store the pointers to each column of V for speed # Upper Hessenberg matrix, converted to upper tri with Givens Rots H = np . zeros ( ( maxiter + 1 , maxiter + 1 ) , dtype = xtype ) # 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 = 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 # Apply diagonal preconditioner if weighting in ( 'local' , '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 = np . sqrt ( ( R . data . conjugate ( ) * R . data ) . sum ( ) ) g = np . 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 = np.ravel((A*T).data) # print "Iteration " + str(i+1) + "   \ # Energy = %1.3e"%np.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(pattern) #   with the added constraint that explicit zeros are in AP wherever #   AP = 0 and pattern does not AV . data [ : ] = 0.0 amg_core . incomplete_mat_mult_bsr ( A . indptr , A . indices , np . ravel ( A . data ) , V [ i ] . indptr , V [ i ] . indices , np . ravel ( V [ i ] . data ) , AV . indptr , AV . indices , np . ravel ( AV . data ) , int ( T . shape [ 0 ] / T . blocksize [ 0 ] ) , int ( T . shape [ 1 ] / T . blocksize [ 1 ] ) , A . blocksize [ 0 ] , A . blocksize [ 1 ] , T . blocksize [ 1 ] ) if weighting in ( 'local' , 'diagonal' ) : AV = scale_rows ( AV , Dinv ) else : AV = Dinv * AV # Enforce AV*B = 0 satisfy_constraints ( AV , B , BtBinv ) V . append ( AV . copy ( ) ) # Modified Gram-Schmidt for j in range ( i + 1 ) : # Frobenius inner-product H [ j , i ] = ( V [ j ] . conjugate ( ) . multiply ( V [ i + 1 ] ) ) . sum ( ) V [ i + 1 ] = V [ i + 1 ] - H [ j , i ] * V [ j ] # Frobenius Norm H [ i + 1 , i ] = np . sqrt ( ( V [ i + 1 ] . data . conjugate ( ) * V [ i + 1 ] . data ) . sum ( ) ) # Check for breakdown if H [ i + 1 , i ] != 0.0 : V [ i + 1 ] = ( 1.0 / H [ i + 1 , i ] ) * V [ i + 1 ] # Apply previous Givens rotations to H if i > 0 : apply_givens ( Q , H [ : , i ] , i ) # Calculate and apply next complex-valued Givens Rotation if H [ i + 1 , i ] != 0 : h1 = H [ i , i ] h2 = H [ i + 1 , i ] h1_mag = np . abs ( h1 ) h2_mag = np . abs ( h2 ) if h1_mag < h2_mag : mu = h1 / h2 tau = np . conjugate ( mu ) / np . abs ( mu ) else : mu = h2 / h1 tau = mu / np . abs ( mu ) denom = np . sqrt ( h1_mag ** 2 + h2_mag ** 2 ) c = h1_mag / denom s = h2_mag * tau / denom Qblock = np . array ( [ [ c , np . conjugate ( s ) ] , [ - s , c ] ] , dtype = xtype ) Q . append ( Qblock ) # Apply Givens Rotation to g, #   the RHS for the linear system in the Krylov Subspace. g [ i : i + 2 ] = np . dot ( Qblock , g [ i : i + 2 ] ) # Apply effect of Givens Rotation to H H [ i , i ] = np . dot ( Qblock [ 0 , : ] , H [ i : i + 2 , i ] ) H [ i + 1 , i ] = 0.0 normr = np . abs ( g [ i + 1 ] ) # print "Iteration " + str(i+1) + "   Normr  %1.3e"%normr # End while loop # Find best update to x in Krylov Space, V.  Solve (i x i) system. if i != - 1 : y = la . solve ( H [ 0 : i + 1 , 0 : i + 1 ] , g [ 0 : i + 1 ] ) for j in range ( i + 1 ) : T = T + y [ j ] * V [ j ] # vect = np.ravel((A*T).data) # print "Final Iteration " + str(i) + "   \ # Energy = %1.3e"%np.sqrt( (vect.conjugate()*vect).sum() ) # Ensure identity at C-pts if Cpt_params [ 0 ] : T = Cpt_params [ 1 ] [ 'I_F' ] * T + Cpt_params [ 1 ] [ 'P_I' ] return T