(ml, symmetry, candidate_iters, prepostsmoother,
                        smooth, eliminate_local, coarse_solver, work)
return x.reshape(-1, 1)
Compute additional candidates and improvements.

    Parameters
    ----------
    candidate_iters
        number of test relaxation iterations
    epsilon
        minimum acceptable relaxation convergence factor

    Notes
    -----
    Follows Algorithm 4 in Brezina et al.

    References
    ----------
    .. [1] Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge
       "Adaptive Smoothed Aggregation (alphaSA) Multigrid"
       SIAM Review Volume 47,  Issue 2  (2005)
       http://www.cs.umn.edu/~maclach/research/aSA2.pdf
Compute additional candidates and improvements.
Compute additional candidates and improvements .
def general_setup_stage(ml, symmetry, candidate_iters, prepostsmoother,
                        smooth, eliminate_local, coarse_solver, work):
    """Compute additional candidates and improvements.

    Parameters
    ----------
    candidate_iters
        number of test relaxation iterations
    epsilon
        minimum acceptable relaxation convergence factor

    Notes
    -----
    Follows Algorithm 4 in Brezina et al.

    References
    ----------
    .. [1] Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge
       "Adaptive Smoothed Aggregation (alphaSA) Multigrid"
       SIAM Review Volume 47,  Issue 2  (2005)
       http://www.cs.umn.edu/~maclach/research/aSA2.pdf

    """
    def make_bridge(T):
        M, N = T.shape
        K = T.blocksize[0]
        bnnz = T.indptr[-1]
        # the K+1 represents the new dof introduced by the new candidate.  the
        # bridge 'T' ignores this new dof and just maps zeros there
        data = np.zeros((bnnz, K+1, K), dtype=T.dtype)
        data[:, :-1, :] = T.data
        return bsr_matrix((data, T.indices, T.indptr),
                          shape=((K + 1) * int(M / K), N))

    def expand_candidates(B_old, nodesize):  # pylint: disable=unused-variable
        # insert a new dof that is always zero, to create NullDim+1 dofs per
        # node in B
        NullDim = B_old.shape[1]
        nnodes = int(B_old.shape[0] / nodesize)
        Bnew = np.zeros((nnodes, nodesize+1, NullDim), dtype=B_old.dtype)
        Bnew[:, :-1, :] = B_old.reshape(nnodes, nodesize, NullDim)
        return Bnew.reshape(-1, NullDim)

    levels = ml.levels

    x = np.random.rand(levels[0].A.shape[0], 1)
    if levels[0].A.dtype.name.startswith('complex'):
        x = x + 1.0j*np.random.rand(levels[0].A.shape[0], 1)
    b = np.zeros_like(x)

    x = ml.solve(b, x0=x, tol=1e-20,
                 maxiter=candidate_iters)
    work[:] += ml.operator_complexity()*ml.levels[0].A.nnz*candidate_iters*2

    T0 = levels[0].T.copy()

    # TEST FOR CONVERGENCE HERE

    for i in range(len(ml.levels) - 2):
        # alpha-SA paper does local elimination here, but after talking
        # to Marian, its not clear that this helps things
        # fn, kwargs = unpack_arg(eliminate_local)
        # if fn == True:
        #    eliminate_local_candidates(x,levels[i].AggOp,levels[i].A,
        #    levels[i].T, **kwargs)

        # add candidate to B
        B = np.hstack((levels[i].B, x.reshape(-1, 1)))

        # construct Ptent
        T, R = fit_candidates(levels[i].AggOp, B)

        levels[i].T = T
        x = R[:, -1].reshape(-1, 1)

        # smooth P
        fn, kwargs = unpack_arg(smooth[i])
        if fn == 'jacobi':
            levels[i].P = jacobi_prolongation_smoother(levels[i].A, T,
                                                       levels[i].C, R,
                                                       **kwargs)
        elif fn == 'richardson':
            levels[i].P = richardson_prolongation_smoother(levels[i].A, T,
                                                           **kwargs)
        elif fn == 'energy':
            levels[i].P = energy_prolongation_smoother(levels[i].A, T,
                                                       levels[i].C, R, None,
                                                       (False, {}), **kwargs)
            x = R[:, -1].reshape(-1, 1)
        elif fn is None:
            levels[i].P = T
        else:
            raise ValueError(f'Unrecognized prolongation smoother method: {str(fn)}')

        # construct R
        if symmetry == 'symmetric':  # R should reflect A's structure
            levels[i].R = levels[i].P.T.asformat(levels[i].P.format)
        elif symmetry == 'hermitian':
            levels[i].R = levels[i].P.H.asformat(levels[i].P.format)

        # construct coarse A
        levels[i+1].A = levels[i].R * levels[i].A * levels[i].P

        # construct bridging P
        T_bridge = make_bridge(levels[i+1].T)
        R_bridge = levels[i+2].B

        # smooth bridging P
        fn, kwargs = unpack_arg(smooth[i+1])
        if fn == 'jacobi':
            levels[i+1].P = jacobi_prolongation_smoother(levels[i+1].A,
                                                         T_bridge,
                                                         levels[i+1].C,
                                                         R_bridge, **kwargs)
        elif fn == 'richardson':
            levels[i+1].P = richardson_prolongation_smoother(levels[i+1].A,
                                                             T_bridge,
                                                             **kwargs)
        elif fn == 'energy':
            levels[i+1].P = energy_prolongation_smoother(levels[i+1].A,
                                                         T_bridge,
                                                         levels[i+1].C,
                                                         R_bridge, None,
                                                         (False, {}), **kwargs)
        elif fn is None:
            levels[i+1].P = T_bridge
        else:
            raise ValueError(f'Unrecognized prolongation smoother method {str(fn)}')

        # construct the 'bridging' R
        if symmetry == 'symmetric':  # R should reflect A's structure
            levels[i+1].R = levels[i+1].P.T.asformat(levels[i+1].P.format)
        elif symmetry == 'hermitian':
            levels[i+1].R = levels[i+1].P.H.asformat(levels[i+1].P.format)

        # run solver on candidate
        solver = MultilevelSolver(levels[i+1:], coarse_solver=coarse_solver)
        change_smoothers(solver, presmoother=prepostsmoother,
                         postsmoother=prepostsmoother)
        x = solver.solve(np.zeros_like(x), x0=x,
                         tol=1e-20,
                         maxiter=candidate_iters)
        work[:] += 2 * solver.operator_complexity() * solver.levels[0].A.nnz *\
            candidate_iters*2

        # update values on next level
        levels[i+1].B = R[:, :-1].copy()
        levels[i+1].T = T_bridge

    # note that we only use the x from the second coarsest level
    fn, kwargs = unpack_arg(prepostsmoother)
    for lvl in reversed(levels[:-2]):
        x = lvl.P * x
        work[:] += lvl.A.nnz*candidate_iters*2

        if fn == 'gauss_seidel':
            # only relax at nonzeros, so as not to mess up any locally dropped
            # candidates
            indices = np.ravel(x).nonzero()[0]
            gauss_seidel_indexed(lvl.A, x, np.zeros_like(x), indices,
                                 iterations=candidate_iters, sweep='symmetric')

        elif fn == 'gauss_seidel_ne':
            gauss_seidel_ne(lvl.A, x, np.zeros_like(x),
                            iterations=candidate_iters, sweep='symmetric')

        elif fn == 'gauss_seidel_nr':
            gauss_seidel_nr(lvl.A, x, np.zeros_like(x),
                            iterations=candidate_iters, sweep='symmetric')

        elif fn == 'jacobi':
            jacobi(lvl.A, x, np.zeros_like(x), iterations=1,
                   omega=1.0 / rho_D_inv_A(lvl.A))

        elif fn == 'richardson':
            polynomial(lvl.A, x, np.zeros_like(x), iterations=1,
                       coefficients=[1.0/approximate_spectral_radius(lvl.A)])

        elif fn == 'gmres':
            x[:] = (gmres(lvl.A, np.zeros_like(x), x0=x,
                          maxiter=candidate_iters)[0]).reshape(x.shape)
        else:
            raise TypeError('Unrecognized smoother')

    # x will be dense again, so we have to drop locally again
    elim, elim_kwargs = unpack_arg(eliminate_local)
    if elim is True:
        x = x/norm(x, 'inf')
        eliminate_local_candidates(x, levels[0].AggOp, levels[0].A, T0,
                                   **elim_kwargs)

    return x.reshape(-1, 1)
def general_setup_stage ( ml , symmetry , candidate_iters , prepostsmoother , smooth , eliminate_local , coarse_solver , work ) : def make_bridge ( T ) : M , N = T . shape K = T . blocksize [ 0 ] bnnz = T . indptr [ - 1 ] # the K+1 represents the new dof introduced by the new candidate.  the # bridge 'T' ignores this new dof and just maps zeros there data = np . zeros ( ( bnnz , K + 1 , K ) , dtype = T . dtype ) data [ : , : - 1 , : ] = T . data return bsr_matrix ( ( data , T . indices , T . indptr ) , shape = ( ( K + 1 ) * int ( M / K ) , N ) ) def expand_candidates ( B_old , nodesize ) : # pylint: disable=unused-variable # insert a new dof that is always zero, to create NullDim+1 dofs per # node in B NullDim = B_old . shape [ 1 ] nnodes = int ( B_old . shape [ 0 ] / nodesize ) Bnew = np . zeros ( ( nnodes , nodesize + 1 , NullDim ) , dtype = B_old . dtype ) Bnew [ : , : - 1 , : ] = B_old . reshape ( nnodes , nodesize , NullDim ) return Bnew . reshape ( - 1 , NullDim ) levels = ml . levels x = np . random . rand ( levels [ 0 ] . A . shape [ 0 ] , 1 ) if levels [ 0 ] . A . dtype . name . startswith ( 'complex' ) : x = x + 1.0j * np . random . rand ( levels [ 0 ] . A . shape [ 0 ] , 1 ) b = np . zeros_like ( x ) x = ml . solve ( b , x0 = x , tol = 1e-20 , maxiter = candidate_iters ) work [ : ] += ml . operator_complexity ( ) * ml . levels [ 0 ] . A . nnz * candidate_iters * 2 T0 = levels [ 0 ] . T . copy ( ) # TEST FOR CONVERGENCE HERE for i in range ( len ( ml . levels ) - 2 ) : # alpha-SA paper does local elimination here, but after talking # to Marian, its not clear that this helps things # fn, kwargs = unpack_arg(eliminate_local) # if fn == True: #    eliminate_local_candidates(x,levels[i].AggOp,levels[i].A, #    levels[i].T, **kwargs) # add candidate to B B = np . hstack ( ( levels [ i ] . B , x . reshape ( - 1 , 1 ) ) ) # construct Ptent T , R = fit_candidates ( levels [ i ] . AggOp , B ) levels [ i ] . T = T x = R [ : , - 1 ] . reshape ( - 1 , 1 ) # smooth P fn , kwargs = unpack_arg ( smooth [ i ] ) if fn == 'jacobi' : levels [ i ] . P = jacobi_prolongation_smoother ( levels [ i ] . A , T , levels [ i ] . C , R , * * kwargs ) elif fn == 'richardson' : levels [ i ] . P = richardson_prolongation_smoother ( levels [ i ] . A , T , * * kwargs ) elif fn == 'energy' : levels [ i ] . P = energy_prolongation_smoother ( levels [ i ] . A , T , levels [ i ] . C , R , None , ( False , { } ) , * * kwargs ) x = R [ : , - 1 ] . reshape ( - 1 , 1 ) elif fn is None : levels [ i ] . P = T else : raise ValueError ( f'Unrecognized prolongation smoother method: {str(fn)}' ) # construct R if symmetry == 'symmetric' : # R should reflect A's structure levels [ i ] . R = levels [ i ] . P . T . asformat ( levels [ i ] . P . format ) elif symmetry == 'hermitian' : levels [ i ] . R = levels [ i ] . P . H . asformat ( levels [ i ] . P . format ) # construct coarse A levels [ i + 1 ] . A = levels [ i ] . R * levels [ i ] . A * levels [ i ] . P # construct bridging P T_bridge = make_bridge ( levels [ i + 1 ] . T ) R_bridge = levels [ i + 2 ] . B # smooth bridging P fn , kwargs = unpack_arg ( smooth [ i + 1 ] ) if fn == 'jacobi' : levels [ i + 1 ] . P = jacobi_prolongation_smoother ( levels [ i + 1 ] . A , T_bridge , levels [ i + 1 ] . C , R_bridge , * * kwargs ) elif fn == 'richardson' : levels [ i + 1 ] . P = richardson_prolongation_smoother ( levels [ i + 1 ] . A , T_bridge , * * kwargs ) elif fn == 'energy' : levels [ i + 1 ] . P = energy_prolongation_smoother ( levels [ i + 1 ] . A , T_bridge , levels [ i + 1 ] . C , R_bridge , None , ( False , { } ) , * * kwargs ) elif fn is None : levels [ i + 1 ] . P = T_bridge else : raise ValueError ( f'Unrecognized prolongation smoother method {str(fn)}' ) # construct the 'bridging' R if symmetry == 'symmetric' : # R should reflect A's structure levels [ i + 1 ] . R = levels [ i + 1 ] . P . T . asformat ( levels [ i + 1 ] . P . format ) elif symmetry == 'hermitian' : levels [ i + 1 ] . R = levels [ i + 1 ] . P . H . asformat ( levels [ i + 1 ] . P . format ) # run solver on candidate solver = MultilevelSolver ( levels [ i + 1 : ] , coarse_solver = coarse_solver ) change_smoothers ( solver , presmoother = prepostsmoother , postsmoother = prepostsmoother ) x = solver . solve ( np . zeros_like ( x ) , x0 = x , tol = 1e-20 , maxiter = candidate_iters ) work [ : ] += 2 * solver . operator_complexity ( ) * solver . levels [ 0 ] . A . nnz * candidate_iters * 2 # update values on next level levels [ i + 1 ] . B = R [ : , : - 1 ] . copy ( ) levels [ i + 1 ] . T = T_bridge # note that we only use the x from the second coarsest level fn , kwargs = unpack_arg ( prepostsmoother ) for lvl in reversed ( levels [ : - 2 ] ) : x = lvl . P * x work [ : ] += lvl . A . nnz * candidate_iters * 2 if fn == 'gauss_seidel' : # only relax at nonzeros, so as not to mess up any locally dropped # candidates indices = np . ravel ( x ) . nonzero ( ) [ 0 ] gauss_seidel_indexed ( lvl . A , x , np . zeros_like ( x ) , indices , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'gauss_seidel_ne' : gauss_seidel_ne ( lvl . A , x , np . zeros_like ( x ) , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'gauss_seidel_nr' : gauss_seidel_nr ( lvl . A , x , np . zeros_like ( x ) , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'jacobi' : jacobi ( lvl . A , x , np . zeros_like ( x ) , iterations = 1 , omega = 1.0 / rho_D_inv_A ( lvl . A ) ) elif fn == 'richardson' : polynomial ( lvl . A , x , np . zeros_like ( x ) , iterations = 1 , coefficients = [ 1.0 / approximate_spectral_radius ( lvl . A ) ] ) elif fn == 'gmres' : x [ : ] = ( gmres ( lvl . A , np . zeros_like ( x ) , x0 = x , maxiter = candidate_iters ) [ 0 ] ) . reshape ( x . shape ) else : raise TypeError ( 'Unrecognized smoother' ) # x will be dense again, so we have to drop locally again elim , elim_kwargs = unpack_arg ( eliminate_local ) if elim is True : x = x / norm ( x , 'inf' ) eliminate_local_candidates ( x , levels [ 0 ] . AggOp , levels [ 0 ] . A , T0 , * * elim_kwargs ) return x . reshape ( - 1 , 1 )