(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, Cnodes = standard_aggregation(C, **kwargs)
    elif fn == 'naive':
        AggOp, Cnodes = naive_aggregation(C, **kwargs)
    elif fn == 'lloyd':
        AggOp, Cnodes = lloyd_aggregation(C, **kwargs)
    elif fn == 'predefined':
        AggOp = kwargs['AggOp'].tocsr()
        Cnodes = kwargs['Cnodes']
    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[:, 0:get_blocksize(A)] = T B_coarse[:, 0:get_blocksize(A)].
    T, dummy = fit_candidates(AggOp, B[:, 0:get_blocksize(A)])
    del dummy
    if A.symmetry == 'nonsymmetric':
        TH, dummyH = fit_candidates(AggOp, BH[:, 0:get_blocksize(A)])
        del dummyH

    # Create necessary root node matrices
    Cpt_params = (True, get_Cpt_params(A, Cnodes, AggOp, T))
    T = scale_T(T, Cpt_params[1]['P_I'], Cpt_params[1]['I_F'])
    if A.symmetry == 'nonsymmetric':
        TH = scale_T(TH, Cpt_params[1]['P_I'], Cpt_params[1]['I_F'])

    # Set coarse grid near nullspace modes as injected fine grid near
    # null-space modes
    B = Cpt_params[1]['P_I'].T*levels[-1].B
    if A.symmetry == 'nonsymmetric':
        BH = Cpt_params[1]['P_I'].T*levels[-1].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 == 'energy':
        P = energy_prolongation_smoother(A, T, C, B, levels[-1].B,
                                         Cpt_params=Cpt_params, **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 == 'energy':
            R = energy_prolongation_smoother(AH, TH, C, BH, levels[-1].BH,
                                             Cpt_params=Cpt_params, **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].Fpts = Cpt_params[1]['Fpts']  # Fpts
        levels[-1].P_I = Cpt_params[1]['P_I']    # Injection operator
        levels[-1].I_F = Cpt_params[1]['I_F']    # Identity on F-pts
        levels[-1].I_C = Cpt_params[1]['I_C']    # Identity on C-pts

    levels[-1].P = P                             # smoothed prolongator
    levels[-1].R = R                             # restriction operator
    levels[-1].Cpts = Cpt_params[1]['Cpts']      # Cpts (i.e., rootnodes)

    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 , Cnodes = standard_aggregation ( C , * * kwargs ) elif fn == 'naive' : AggOp , Cnodes = naive_aggregation ( C , * * kwargs ) elif fn == 'lloyd' : AggOp , Cnodes = lloyd_aggregation ( C , * * kwargs ) elif fn == 'predefined' : AggOp = kwargs [ 'AggOp' ] . tocsr ( ) Cnodes = kwargs [ 'Cnodes' ] 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[:, 0:get_blocksize(A)] = T B_coarse[:, 0:get_blocksize(A)]. T , dummy = fit_candidates ( AggOp , B [ : , 0 : get_blocksize ( A ) ] ) del dummy if A . symmetry == 'nonsymmetric' : TH , dummyH = fit_candidates ( AggOp , BH [ : , 0 : get_blocksize ( A ) ] ) del dummyH # Create necessary root node matrices Cpt_params = ( True , get_Cpt_params ( A , Cnodes , AggOp , T ) ) T = scale_T ( T , Cpt_params [ 1 ] [ 'P_I' ] , Cpt_params [ 1 ] [ 'I_F' ] ) if A . symmetry == 'nonsymmetric' : TH = scale_T ( TH , Cpt_params [ 1 ] [ 'P_I' ] , Cpt_params [ 1 ] [ 'I_F' ] ) # Set coarse grid near nullspace modes as injected fine grid near # null-space modes B = Cpt_params [ 1 ] [ 'P_I' ] . T * levels [ - 1 ] . B if A . symmetry == 'nonsymmetric' : BH = Cpt_params [ 1 ] [ 'P_I' ] . T * levels [ - 1 ] . 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 == 'energy' : P = energy_prolongation_smoother ( A , T , C , B , levels [ - 1 ] . B , Cpt_params = Cpt_params , * * 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 == 'energy' : R = energy_prolongation_smoother ( AH , TH , C , BH , levels [ - 1 ] . BH , Cpt_params = Cpt_params , * * 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 ] . Fpts = Cpt_params [ 1 ] [ 'Fpts' ] # Fpts levels [ - 1 ] . P_I = Cpt_params [ 1 ] [ 'P_I' ] # Injection operator levels [ - 1 ] . I_F = Cpt_params [ 1 ] [ 'I_F' ] # Identity on F-pts levels [ - 1 ] . I_C = Cpt_params [ 1 ] [ 'I_C' ] # Identity on C-pts levels [ - 1 ] . P = P # smoothed prolongator levels [ - 1 ] . R = R # restriction operator levels [ - 1 ] . Cpts = Cpt_params [ 1 ] [ 'Cpts' ] # Cpts (i.e., rootnodes) 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