(levels, strength, aggregate, smooth, improve_candidates,
                      diagonal_dominance=False, keep=True)
Extend the multigrid hierarchy.

    Service routine to implement the strength of connection, aggregation,
    tentative prolongation construction, and prolongation smoothing.  Called by
    smoothed_aggregation_solver.
Extend the multigrid hierarchy.
Extend the multigrid hierarchy .
def _extend_hierarchy(levels, strength, aggregate, smooth, improve_candidates,
                      diagonal_dominance=False, keep=True):
    """Extend the multigrid hierarchy.

    Service routine to implement the strength of connection, aggregation,
    tentative prolongation construction, and prolongation smoothing.  Called by
    smoothed_aggregation_solver.

    """
    def unpack_arg(v):
        if isinstance(v, tuple):
            return v[0], v[1]
        return v, {}

    A = levels[-1].A
    B = levels[-1].B
    if A.symmetry == 'nonsymmetric':
        AH = A.H.asformat(A.format)
        BH = levels[-1].BH

    # Compute the strength-of-connection matrix C, where larger
    # C[i,j] denote stronger couplings between i and j.
    fn, kwargs = unpack_arg(strength[len(levels)-1])
    if fn == 'symmetric':
        C = symmetric_strength_of_connection(A, **kwargs)
    elif fn == 'classical':
        C = classical_strength_of_connection(A, **kwargs)
    elif fn == 'distance':
        C = distance_strength_of_connection(A, **kwargs)
    elif fn in ('ode', 'evolution'):
        if 'B' in kwargs:
            C = evolution_strength_of_connection(A, **kwargs)
        else:
            C = evolution_strength_of_connection(A, B, **kwargs)
    elif fn == 'energy_based':
        C = energy_based_strength_of_connection(A, **kwargs)
    elif fn == 'predefined':
        C = kwargs['C'].tocsr()
    elif fn == 'algebraic_distance':
        C = algebraic_distance(A, **kwargs)
    elif fn == 'affinity':
        C = affinity_distance(A, **kwargs)
    elif fn is None:
        C = A.tocsr()
    else:
        raise ValueError(f'Unrecognized strength of connection method: {str(fn)}')

    # Avoid coarsening diagonally dominant rows
    flag, kwargs = unpack_arg(diagonal_dominance)
    if flag:
        C = eliminate_diag_dom_nodes(A, C, **kwargs)

    # Compute the aggregation matrix AggOp (i.e., the nodal coarsening of A).
    # AggOp is a boolean matrix, where the sparsity pattern for the k-th column
    # denotes the fine-grid nodes agglomerated into k-th coarse-grid node.
    fn, kwargs = unpack_arg(aggregate[len(levels)-1])
    if fn == 'standard':
        AggOp = standard_aggregation(C, **kwargs)[0]
    elif fn == 'naive':
        AggOp = naive_aggregation(C, **kwargs)[0]
    elif fn == 'lloyd':
        AggOp = lloyd_aggregation(C, **kwargs)[0]
    elif fn == 'predefined':
        AggOp = kwargs['AggOp'].tocsr()
    else:
        raise ValueError(f'Unrecognized aggregation method {str(fn)}')

    # Improve near nullspace candidates by relaxing on A B = 0
    fn, kwargs = unpack_arg(improve_candidates[len(levels)-1])
    if fn is not None:
        b = np.zeros((A.shape[0], 1), dtype=A.dtype)
        B = relaxation_as_linear_operator((fn, kwargs), A, b) * B
        levels[-1].B = B
        if A.symmetry == 'nonsymmetric':
            BH = relaxation_as_linear_operator((fn, kwargs), AH, b) * BH
            levels[-1].BH = BH

    # Compute the tentative prolongator, T, which is a tentative interpolation
    # matrix from the coarse-grid to the fine-grid.  T exactly interpolates
    # B_fine = T B_coarse.
    T, B = fit_candidates(AggOp, B)
    if A.symmetry == 'nonsymmetric':
        TH, BH = fit_candidates(AggOp, BH)

    # Smooth the tentative prolongator, so that it's accuracy is greatly
    # improved for algebraically smooth error.
    fn, kwargs = unpack_arg(smooth[len(levels)-1])
    if fn == 'jacobi':
        P = jacobi_prolongation_smoother(A, T, C, B, **kwargs)
    elif fn == 'richardson':
        P = richardson_prolongation_smoother(A, T, **kwargs)
    elif fn == 'energy':
        P = energy_prolongation_smoother(A, T, C, B, None, (False, {}),
                                         **kwargs)
    elif fn is None:
        P = T
    else:
        raise ValueError(f'Unrecognized prolongation smoother method {str(fn)}')

    # Compute the restriction matrix, R, which interpolates from the fine-grid
    # to the coarse-grid.  If A is nonsymmetric, then R must be constructed
    # based on A.H.  Otherwise R = P.H or P.T.
    symmetry = A.symmetry
    if symmetry == 'hermitian':
        R = P.H
    elif symmetry == 'symmetric':
        R = P.T
    elif symmetry == 'nonsymmetric':
        fn, kwargs = unpack_arg(smooth[len(levels)-1])
        if fn == 'jacobi':
            R = jacobi_prolongation_smoother(AH, TH, C, BH, **kwargs).H
        elif fn == 'richardson':
            R = richardson_prolongation_smoother(AH, TH, **kwargs).H
        elif fn == 'energy':
            R = energy_prolongation_smoother(AH, TH, C, BH, None, (False, {}),
                                             **kwargs)
            R = R.H
        elif fn is None:
            R = T.H
        else:
            raise ValueError(f'Unrecognized prolongation smoother method {str(fn)}')

    if keep:
        levels[-1].C = C  # strength of connection matrix
        levels[-1].AggOp = AggOp  # aggregation operator
        levels[-1].T = T  # tentative prolongator

    levels[-1].P = P  # smoothed prolongator
    levels[-1].R = R  # restriction operator

    levels.append(MultilevelSolver.Level())
    A = R * A * P              # Galerkin operator
    A.symmetry = symmetry
    levels[-1].A = A
    levels[-1].B = B           # right near nullspace candidates

    if A.symmetry == 'nonsymmetric':
        levels[-1].BH = BH
def _extend_hierarchy ( levels , strength , aggregate , smooth , improve_candidates , diagonal_dominance = False , keep = True ) : def unpack_arg ( v ) : if isinstance ( v , tuple ) : return v [ 0 ] , v [ 1 ] return v , { } A = levels [ - 1 ] . A B = levels [ - 1 ] . B if A . symmetry == 'nonsymmetric' : AH = A . H . asformat ( A . format ) BH = levels [ - 1 ] . BH # Compute the strength-of-connection matrix C, where larger # C[i,j] denote stronger couplings between i and j. fn , kwargs = unpack_arg ( strength [ len ( levels ) - 1 ] ) if fn == 'symmetric' : C = symmetric_strength_of_connection ( A , * * kwargs ) elif fn == 'classical' : C = classical_strength_of_connection ( A , * * kwargs ) elif fn == 'distance' : C = distance_strength_of_connection ( A , * * kwargs ) elif fn in ( 'ode' , 'evolution' ) : if 'B' in kwargs : C = evolution_strength_of_connection ( A , * * kwargs ) else : C = evolution_strength_of_connection ( A , B , * * kwargs ) elif fn == 'energy_based' : C = energy_based_strength_of_connection ( A , * * kwargs ) elif fn == 'predefined' : C = kwargs [ 'C' ] . tocsr ( ) elif fn == 'algebraic_distance' : C = algebraic_distance ( A , * * kwargs ) elif fn == 'affinity' : C = affinity_distance ( A , * * kwargs ) elif fn is None : C = A . tocsr ( ) else : raise ValueError ( f'Unrecognized strength of connection method: {str(fn)}' ) # Avoid coarsening diagonally dominant rows flag , kwargs = unpack_arg ( diagonal_dominance ) if flag : C = eliminate_diag_dom_nodes ( A , C , * * kwargs ) # Compute the aggregation matrix AggOp (i.e., the nodal coarsening of A). # AggOp is a boolean matrix, where the sparsity pattern for the k-th column # denotes the fine-grid nodes agglomerated into k-th coarse-grid node. fn , kwargs = unpack_arg ( aggregate [ len ( levels ) - 1 ] ) if fn == 'standard' : AggOp = standard_aggregation ( C , * * kwargs ) [ 0 ] elif fn == 'naive' : AggOp = naive_aggregation ( C , * * kwargs ) [ 0 ] elif fn == 'lloyd' : AggOp = lloyd_aggregation ( C , * * kwargs ) [ 0 ] elif fn == 'predefined' : AggOp = kwargs [ 'AggOp' ] . tocsr ( ) else : raise ValueError ( f'Unrecognized aggregation method {str(fn)}' ) # Improve near nullspace candidates by relaxing on A B = 0 fn , kwargs = unpack_arg ( improve_candidates [ len ( levels ) - 1 ] ) if fn is not None : b = np . zeros ( ( A . shape [ 0 ] , 1 ) , dtype = A . dtype ) B = relaxation_as_linear_operator ( ( fn , kwargs ) , A , b ) * B levels [ - 1 ] . B = B if A . symmetry == 'nonsymmetric' : BH = relaxation_as_linear_operator ( ( fn , kwargs ) , AH , b ) * BH levels [ - 1 ] . BH = BH # Compute the tentative prolongator, T, which is a tentative interpolation # matrix from the coarse-grid to the fine-grid.  T exactly interpolates # B_fine = T B_coarse. T , B = fit_candidates ( AggOp , B ) if A . symmetry == 'nonsymmetric' : TH , BH = fit_candidates ( AggOp , BH ) # Smooth the tentative prolongator, so that it's accuracy is greatly # improved for algebraically smooth error. fn , kwargs = unpack_arg ( smooth [ len ( levels ) - 1 ] ) if fn == 'jacobi' : P = jacobi_prolongation_smoother ( A , T , C , B , * * kwargs ) elif fn == 'richardson' : P = richardson_prolongation_smoother ( A , T , * * kwargs ) elif fn == 'energy' : P = energy_prolongation_smoother ( A , T , C , B , None , ( False , { } ) , * * kwargs ) elif fn is None : P = T else : raise ValueError ( f'Unrecognized prolongation smoother method {str(fn)}' ) # Compute the restriction matrix, R, which interpolates from the fine-grid # to the coarse-grid.  If A is nonsymmetric, then R must be constructed # based on A.H.  Otherwise R = P.H or P.T. symmetry = A . symmetry if symmetry == 'hermitian' : R = P . H elif symmetry == 'symmetric' : R = P . T elif symmetry == 'nonsymmetric' : fn , kwargs = unpack_arg ( smooth [ len ( levels ) - 1 ] ) if fn == 'jacobi' : R = jacobi_prolongation_smoother ( AH , TH , C , BH , * * kwargs ) . H elif fn == 'richardson' : R = richardson_prolongation_smoother ( AH , TH , * * kwargs ) . H elif fn == 'energy' : R = energy_prolongation_smoother ( AH , TH , C , BH , None , ( False , { } ) , * * kwargs ) R = R . H elif fn is None : R = T . H else : raise ValueError ( f'Unrecognized prolongation smoother method {str(fn)}' ) if keep : levels [ - 1 ] . C = C # strength of connection matrix levels [ - 1 ] . AggOp = AggOp # aggregation operator levels [ - 1 ] . T = T # tentative prolongator levels [ - 1 ] . P = P # smoothed prolongator levels [ - 1 ] . R = R # restriction operator levels . append ( MultilevelSolver . Level ( ) ) A = R * A * P # Galerkin operator A . symmetry = symmetry levels [ - 1 ] . A = A levels [ - 1 ] . B = B # right near nullspace candidates if A . symmetry == 'nonsymmetric' : levels [ - 1 ] . BH = BH