(A, symmetry, pdef, candidate_iters, epsilon,
                        max_levels, max_coarse, aggregate, prepostsmoother,
                        smooth, strength, work, initial_candidate=None)
return x, aggregate, strength
Compute aggregation and the first near-nullspace candidate.

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

    Notes
    -----
    This follows Algorithm 3 in Brezina et al.

    References
    ----------
    .. [1] Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge
       "Adaptive Smoothed Aggregation (aSA) Multigrid"
       SIAM Review Volume 47,  Issue 2  (2005)
       http://www.cs.umn.edu/~maclach/research/aSA2.pdf
Compute aggregation and the first near-nullspace candidate.
Compute aggregation and the first near - nullspace candidate .
def initial_setup_stage(A, symmetry, pdef, candidate_iters, epsilon,
                        max_levels, max_coarse, aggregate, prepostsmoother,
                        smooth, strength, work, initial_candidate=None):
    """Compute aggregation and the first near-nullspace candidate.

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

    Notes
    -----
    This follows Algorithm 3 in Brezina et al.

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

    """
    # Define relaxation routine
    def relax(A, x):
        fn, _ = unpack_arg(prepostsmoother)
        if fn == 'gauss_seidel':
            gauss_seidel(A, x, np.zeros_like(x),
                         iterations=candidate_iters, sweep='symmetric')
        elif fn == 'gauss_seidel_nr':
            gauss_seidel_nr(A, x, np.zeros_like(x),
                            iterations=candidate_iters, sweep='symmetric')
        elif fn == 'gauss_seidel_ne':
            gauss_seidel_ne(A, x, np.zeros_like(x),
                            iterations=candidate_iters, sweep='symmetric')
        elif fn == 'jacobi':
            jacobi(A, x, np.zeros_like(x), iterations=1,
                   omega=1.0 / rho_D_inv_A(A))
        elif fn == 'richardson':
            polynomial(A, x, np.zeros_like(x), iterations=1,
                       coefficients=[1.0/approximate_spectral_radius(A)])
        elif fn == 'gmres':
            x[:] = (gmres(A, np.zeros_like(x), x0=x,
                          maxiter=candidate_iters)[0]).reshape(x.shape)
        else:
            raise TypeError('Unrecognized smoother')

    # flag for skipping steps f-i in step 4
    skip_f_to_i = True

    # step 1
    A_l = A
    if initial_candidate is None:
        x = np.random.rand(A_l.shape[0], 1).astype(A_l.dtype)
        # The following type check matches the usual 'complex' type,
        # but also numpy data types such as 'complex64', 'complex128'
        # and 'complex256'.
        if A_l.dtype.name.startswith('complex'):
            x = x + 1.0j*np.random.rand(A_l.shape[0], 1)
    else:
        x = np.array(initial_candidate, dtype=A_l.dtype)

    # step 2
    relax(A_l, x)
    work[:] += A_l.nnz * candidate_iters*2

    # step 3
    # not advised to stop the iteration here: often the first relaxation pass
    # _is_ good, but the remaining passes are poor
    # if x_A_x/x_A_x_old < epsilon:
    #    # relaxation alone is sufficient
    #    print 'relaxation alone works: %g'%(x_A_x/x_A_x_old)
    #    return x, []

    # step 4
    As = [A]
    xs = [x]
    Ps = []
    AggOps = []
    StrengthOps = []

    while A.shape[0] > max_coarse and max_levels > 1:
        # The real check to break from the while loop is below

        # Begin constructing next level
        fn, kwargs = unpack_arg(strength[len(As)-1])  # step 4b
        if fn == 'symmetric':
            C_l = symmetric_strength_of_connection(A_l, **kwargs)
            # Diagonal must be nonzero
            C_l = C_l + eye(C_l.shape[0], C_l.shape[1], format='csr')
        elif fn == 'classical':
            C_l = classical_strength_of_connection(A_l, **kwargs)
            # Diagonal must be nonzero
            C_l = C_l + eye(C_l.shape[0], C_l.shape[1], format='csr')
            if isspmatrix_bsr(A_l):
                C_l = amalgamate(C_l, A_l.blocksize[0])
        elif fn in ('ode', 'evolution'):
            C_l = evolution_strength_of_connection(A_l,
                                                   np.ones(
                                                       (A_l.shape[0], 1),
                                                       dtype=A.dtype),
                                                   **kwargs)
        elif fn == 'predefined':
            C_l = kwargs['C'].tocsr()
        elif fn is None:
            C_l = A_l.tocsr()
        else:
            raise ValueError(f'Unrecognized strength of connection method: {str(fn)}')

        # In SA, strength represents 'distance', so we take magnitude of
        # complex values
        if C_l.dtype.name.startswith('complex'):
            C_l.data = np.abs(C_l.data)

        # Create a unified strength framework so that large values represent
        # strong connections and small values represent weak connections
        if fn in ('ode', 'evolution', 'energy_based'):
            C_l.data = 1.0 / C_l.data

        # aggregation
        fn, kwargs = unpack_arg(aggregate[len(As) - 1])
        if fn == 'standard':
            AggOp = standard_aggregation(C_l, **kwargs)[0]
        elif fn == 'lloyd':
            AggOp = lloyd_aggregation(C_l, **kwargs)[0]
        elif fn == 'predefined':
            AggOp = kwargs['AggOp'].tocsr()
        else:
            raise ValueError(f'Unrecognized aggregation method {str(fn)}')

        T_l, x = fit_candidates(AggOp, x)  # step 4c

        fn, kwargs = unpack_arg(smooth[len(As)-1])  # step 4d
        if fn == 'jacobi':
            P_l = jacobi_prolongation_smoother(A_l, T_l, C_l, x, **kwargs)
        elif fn == 'richardson':
            P_l = richardson_prolongation_smoother(A_l, T_l, **kwargs)
        elif fn == 'energy':
            P_l = energy_prolongation_smoother(A_l, T_l, C_l, x, None,
                                               (False, {}), **kwargs)
        elif fn is None:
            P_l = T_l
        else:
            raise ValueError(f'Unrecognized prolongation smoother method: {str(fn)}')

        # R should reflect A's structure # step 4e
        if symmetry == 'symmetric':
            A_l = P_l.T.asformat(P_l.format) * A_l * P_l
        elif symmetry == 'hermitian':
            A_l = P_l.H.asformat(P_l.format) * A_l * P_l
        else:
            raise ValueError(f'aSA not implemented for symmetry={symmetry}.')

        StrengthOps.append(C_l)
        AggOps.append(AggOp)
        Ps.append(P_l)
        As.append(A_l)

        # skip to step 5 as in step 4e
        if (A_l.shape[0] <= max_coarse) or (len(AggOps) + 1 >= max_levels):
            break

        if not skip_f_to_i:
            x_hat = x.copy()  # step 4g
            relax(A_l, x)  # step 4h
            work[:] += A_l.nnz*candidate_iters*2
            if pdef is True:
                x_A_x = np.dot(np.conjugate(x).T, A_l*x)
                xhat_A_xhat = np.dot(np.conjugate(x_hat).T, A_l*x_hat)
                err_ratio = (x_A_x/xhat_A_xhat)**(1.0/candidate_iters)
            else:
                # use A.H A inner-product
                Ax = A_l * x
                # Axhat = A_l * x_hat
                x_A_x = np.dot(np.conjugate(Ax).T, Ax)
                xhat_A_xhat = np.dot(np.conjugate(x_hat).T, A_l*x_hat)
                err_ratio = (x_A_x/xhat_A_xhat)**(1.0/candidate_iters)

            if err_ratio < epsilon:  # step 4i
                # print 'sufficient convergence, skipping'
                skip_f_to_i = True
                if x_A_x == 0:
                    x = x_hat  # need to restore x
        else:
            # just carry out relaxation, don't check for convergence
            relax(A_l, x)  # step 4h
            work[:] += 2 * A_l.nnz * candidate_iters

        # store xs for diagnostic use and for use in step 5
        xs.append(x)

    # step 5
    # Extend coarse-level candidate to the finest level
    # --> note that we start with the x from the second coarsest level
    x = xs[-1]
    # make sure that xs[-1] has been relaxed by step 4h, i.e. relax(As[-2], x)
    for lev in range(len(Ps)-2, -1, -1):  # lev = coarsest ... finest-1
        P = Ps[lev]                     # I: lev --> lev+1
        A = As[lev]                     # A on lev+1
        x = P * x
        relax(A, x)
        work[:] += A.nnz*candidate_iters*2

    # Set predefined strength of connection and aggregation
    if len(AggOps) > 1:
        aggregate = [('predefined', {'AggOp': AggOps[i]})
                     for i in range(len(AggOps))]
        strength = [('predefined', {'C': StrengthOps[i]})
                    for i in range(len(StrengthOps))]

    return x, aggregate, strength
def initial_setup_stage ( A , symmetry , pdef , candidate_iters , epsilon , max_levels , max_coarse , aggregate , prepostsmoother , smooth , strength , work , initial_candidate = None ) : # Define relaxation routine def relax ( A , x ) : fn , _ = unpack_arg ( prepostsmoother ) if fn == 'gauss_seidel' : gauss_seidel ( A , x , np . zeros_like ( x ) , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'gauss_seidel_nr' : gauss_seidel_nr ( A , x , np . zeros_like ( x ) , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'gauss_seidel_ne' : gauss_seidel_ne ( A , x , np . zeros_like ( x ) , iterations = candidate_iters , sweep = 'symmetric' ) elif fn == 'jacobi' : jacobi ( A , x , np . zeros_like ( x ) , iterations = 1 , omega = 1.0 / rho_D_inv_A ( A ) ) elif fn == 'richardson' : polynomial ( A , x , np . zeros_like ( x ) , iterations = 1 , coefficients = [ 1.0 / approximate_spectral_radius ( A ) ] ) elif fn == 'gmres' : x [ : ] = ( gmres ( A , np . zeros_like ( x ) , x0 = x , maxiter = candidate_iters ) [ 0 ] ) . reshape ( x . shape ) else : raise TypeError ( 'Unrecognized smoother' ) # flag for skipping steps f-i in step 4 skip_f_to_i = True # step 1 A_l = A if initial_candidate is None : x = np . random . rand ( A_l . shape [ 0 ] , 1 ) . astype ( A_l . dtype ) # The following type check matches the usual 'complex' type, # but also numpy data types such as 'complex64', 'complex128' # and 'complex256'. if A_l . dtype . name . startswith ( 'complex' ) : x = x + 1.0j * np . random . rand ( A_l . shape [ 0 ] , 1 ) else : x = np . array ( initial_candidate , dtype = A_l . dtype ) # step 2 relax ( A_l , x ) work [ : ] += A_l . nnz * candidate_iters * 2 # step 3 # not advised to stop the iteration here: often the first relaxation pass # _is_ good, but the remaining passes are poor # if x_A_x/x_A_x_old < epsilon: #    # relaxation alone is sufficient #    print 'relaxation alone works: %g'%(x_A_x/x_A_x_old) #    return x, [] # step 4 As = [ A ] xs = [ x ] Ps = [ ] AggOps = [ ] StrengthOps = [ ] while A . shape [ 0 ] > max_coarse and max_levels > 1 : # The real check to break from the while loop is below # Begin constructing next level fn , kwargs = unpack_arg ( strength [ len ( As ) - 1 ] ) # step 4b if fn == 'symmetric' : C_l = symmetric_strength_of_connection ( A_l , * * kwargs ) # Diagonal must be nonzero C_l = C_l + eye ( C_l . shape [ 0 ] , C_l . shape [ 1 ] , format = 'csr' ) elif fn == 'classical' : C_l = classical_strength_of_connection ( A_l , * * kwargs ) # Diagonal must be nonzero C_l = C_l + eye ( C_l . shape [ 0 ] , C_l . shape [ 1 ] , format = 'csr' ) if isspmatrix_bsr ( A_l ) : C_l = amalgamate ( C_l , A_l . blocksize [ 0 ] ) elif fn in ( 'ode' , 'evolution' ) : C_l = evolution_strength_of_connection ( A_l , np . ones ( ( A_l . shape [ 0 ] , 1 ) , dtype = A . dtype ) , * * kwargs ) elif fn == 'predefined' : C_l = kwargs [ 'C' ] . tocsr ( ) elif fn is None : C_l = A_l . tocsr ( ) else : raise ValueError ( f'Unrecognized strength of connection method: {str(fn)}' ) # In SA, strength represents 'distance', so we take magnitude of # complex values if C_l . dtype . name . startswith ( 'complex' ) : C_l . data = np . abs ( C_l . data ) # Create a unified strength framework so that large values represent # strong connections and small values represent weak connections if fn in ( 'ode' , 'evolution' , 'energy_based' ) : C_l . data = 1.0 / C_l . data # aggregation fn , kwargs = unpack_arg ( aggregate [ len ( As ) - 1 ] ) if fn == 'standard' : AggOp = standard_aggregation ( C_l , * * kwargs ) [ 0 ] elif fn == 'lloyd' : AggOp = lloyd_aggregation ( C_l , * * kwargs ) [ 0 ] elif fn == 'predefined' : AggOp = kwargs [ 'AggOp' ] . tocsr ( ) else : raise ValueError ( f'Unrecognized aggregation method {str(fn)}' ) T_l , x = fit_candidates ( AggOp , x ) # step 4c fn , kwargs = unpack_arg ( smooth [ len ( As ) - 1 ] ) # step 4d if fn == 'jacobi' : P_l = jacobi_prolongation_smoother ( A_l , T_l , C_l , x , * * kwargs ) elif fn == 'richardson' : P_l = richardson_prolongation_smoother ( A_l , T_l , * * kwargs ) elif fn == 'energy' : P_l = energy_prolongation_smoother ( A_l , T_l , C_l , x , None , ( False , { } ) , * * kwargs ) elif fn is None : P_l = T_l else : raise ValueError ( f'Unrecognized prolongation smoother method: {str(fn)}' ) # R should reflect A's structure # step 4e if symmetry == 'symmetric' : A_l = P_l . T . asformat ( P_l . format ) * A_l * P_l elif symmetry == 'hermitian' : A_l = P_l . H . asformat ( P_l . format ) * A_l * P_l else : raise ValueError ( f'aSA not implemented for symmetry={symmetry}.' ) StrengthOps . append ( C_l ) AggOps . append ( AggOp ) Ps . append ( P_l ) As . append ( A_l ) # skip to step 5 as in step 4e if ( A_l . shape [ 0 ] <= max_coarse ) or ( len ( AggOps ) + 1 >= max_levels ) : break if not skip_f_to_i : x_hat = x . copy ( ) # step 4g relax ( A_l , x ) # step 4h work [ : ] += A_l . nnz * candidate_iters * 2 if pdef is True : x_A_x = np . dot ( np . conjugate ( x ) . T , A_l * x ) xhat_A_xhat = np . dot ( np . conjugate ( x_hat ) . T , A_l * x_hat ) err_ratio = ( x_A_x / xhat_A_xhat ) ** ( 1.0 / candidate_iters ) else : # use A.H A inner-product Ax = A_l * x # Axhat = A_l * x_hat x_A_x = np . dot ( np . conjugate ( Ax ) . T , Ax ) xhat_A_xhat = np . dot ( np . conjugate ( x_hat ) . T , A_l * x_hat ) err_ratio = ( x_A_x / xhat_A_xhat ) ** ( 1.0 / candidate_iters ) if err_ratio < epsilon : # step 4i # print 'sufficient convergence, skipping' skip_f_to_i = True if x_A_x == 0 : x = x_hat # need to restore x else : # just carry out relaxation, don't check for convergence relax ( A_l , x ) # step 4h work [ : ] += 2 * A_l . nnz * candidate_iters # store xs for diagnostic use and for use in step 5 xs . append ( x ) # step 5 # Extend coarse-level candidate to the finest level # --> note that we start with the x from the second coarsest level x = xs [ - 1 ] # make sure that xs[-1] has been relaxed by step 4h, i.e. relax(As[-2], x) for lev in range ( len ( Ps ) - 2 , - 1 , - 1 ) : # lev = coarsest ... finest-1 P = Ps [ lev ] # I: lev --> lev+1 A = As [ lev ] # A on lev+1 x = P * x relax ( A , x ) work [ : ] += A . nnz * candidate_iters * 2 # Set predefined strength of connection and aggregation if len ( AggOps ) > 1 : aggregate = [ ( 'predefined' , { 'AggOp' : AggOps [ i ] } ) for i in range ( len ( AggOps ) ) ] strength = [ ( 'predefined' , { 'C' : StrengthOps [ i ] } ) for i in range ( len ( StrengthOps ) ) ] return x , aggregate , strength