diff --git a/pyomo/contrib/incidence_analysis/__init__.py b/pyomo/contrib/incidence_analysis/__init__.py
index 178bbc96db8..4dbfc66d25f 100644
--- a/pyomo/contrib/incidence_analysis/__init__.py
+++ b/pyomo/contrib/incidence_analysis/__init__.py
@@ -1,7 +1,10 @@
 from .triangularize import block_triangularize
 from .matching import maximum_matching
-from .interface import IncidenceGraphInterface
+from .interface import (
+    IncidenceGraphInterface,
+    get_bipartite_incidence_graph,
+)
 from .util import (
-        generate_strongly_connected_components,
-        solve_strongly_connected_components,
-        )
+    generate_strongly_connected_components,
+    solve_strongly_connected_components,
+)
diff --git a/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py b/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py
index 347538659d4..b861cce356a 100644
--- a/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py
+++ b/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py
@@ -27,8 +27,7 @@
 
 
 def _get_projected_digraph(bg, matching, top_nodes):
-    """
-    """
+    """ """
     digraph = DiGraph()
     digraph.add_nodes_from(top_nodes)
     for n in top_nodes:
@@ -45,8 +44,7 @@ def _get_projected_digraph(bg, matching, top_nodes):
 
 
 def _get_reachable_from(digraph, sources):
-    """
-    """
+    """ """
     _filter = set()
     reachable = []
     for node in sources:
@@ -96,14 +94,17 @@ def dulmage_mendelsohn(bg, top_nodes=None, matching=None):
     t_other = [t for t in top_nodes if t not in _filter]
     b_other = [b for b in bot_nodes if b not in _filter]
 
-    return ((
-        t_unmatched,
-        t_reachable,
-        t_matched_with_reachable,
-        t_other,
-        ), (
-        b_unmatched,
-        b_reachable,
-        b_matched_with_reachable,
-        b_other,
-        ))
+    return (
+        (
+            t_unmatched,
+            t_reachable,
+            t_matched_with_reachable,
+            t_other,
+        ),
+        (
+            b_unmatched,
+            b_reachable,
+            b_matched_with_reachable,
+            b_other,
+        ),
+    )
diff --git a/pyomo/contrib/incidence_analysis/connected.py b/pyomo/contrib/incidence_analysis/connected.py
index f268e292991..67b7e0b7392 100644
--- a/pyomo/contrib/incidence_analysis/connected.py
+++ b/pyomo/contrib/incidence_analysis/connected.py
@@ -13,8 +13,7 @@
 
 
 def get_independent_submatrices(matrix):
-    """
-    """
+    """ """
     nxc = nx.algorithms.components
     nxb = nx.algorithms.bipartite
     from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix
@@ -27,11 +26,10 @@ def get_independent_submatrices(matrix):
     # nodes have values in [N, N+M-1].
     # We could also check the "bipartite" attribute of each node...
     row_blocks = [
-        sorted([node for node in comp if node < N])
-        for comp in connected_components
+        sorted([node for node in comp if node < N]) for comp in connected_components
     ]
     col_blocks = [
-        sorted([node-N for node in comp if node >= N])
+        sorted([node - N for node in comp if node >= N])
         for comp in connected_components
     ]
     return row_blocks, col_blocks
diff --git a/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py b/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py
index 761b6ee69c2..64d77c64266 100644
--- a/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py
+++ b/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py
@@ -13,20 +13,22 @@
 from pyomo.common.dependencies import networkx as nx
 from pyomo.contrib.incidence_analysis.common.dulmage_mendelsohn import (
     dulmage_mendelsohn as dm_nx,
-    )
+)
+
 """
 This module imports the general Dulmage-Mendelsohn-on-a-graph function
 from "common" and implements an interface for coo_matrix-like objects.
 """
 
 RowPartition = namedtuple(
-        "RowPartition",
-        ["unmatched", "overconstrained", "underconstrained", "square"],
-        )
+    "RowPartition",
+    ["unmatched", "overconstrained", "underconstrained", "square"],
+)
 ColPartition = namedtuple(
-        "ColPartition",
-        ["unmatched", "underconstrained", "overconstrained", "square"],
-        )
+    "ColPartition",
+    ["unmatched", "underconstrained", "overconstrained", "square"],
+)
+
 
 def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None):
     """
@@ -45,9 +47,9 @@ def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None):
         graph = matrix_or_graph
         if top_nodes is None:
             raise ValueError(
-                    "top_nodes must be specified if a graph is provided,"
-                    "\notherwise the result is ambiguous."
-                    )
+                "top_nodes must be specified if a graph is provided,"
+                "\notherwise the result is ambiguous."
+            )
         partition = dm_nx(graph, top_nodes=top_nodes, matching=matching)
         # RowPartition and ColPartition do not make sense for a general graph.
         # However, here we assume that this graph comes from a Pyomo model,
@@ -74,17 +76,17 @@ def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None):
         # Matrix rows have bipartite=0, columns have bipartite=1
         bg = from_biadjacency_matrix(matrix)
         row_partition, col_partition = dm_nx(
-                bg,
-                top_nodes=list(range(M)),
-                matching=matching,
-                )
+            bg,
+            top_nodes=list(range(M)),
+            matching=matching,
+        )
 
         partition = (
-                row_partition,
-                tuple([n - M for n in subset] for subset in col_partition)
-                # Column nodes have values in [M, M+N-1]. Apply the offset
-                # to get values corresponding to indices in user's matrix.
-                )
+            row_partition,
+            tuple([n - M for n in subset] for subset in col_partition)
+            # Column nodes have values in [M, M+N-1]. Apply the offset
+            # to get values corresponding to indices in user's matrix.
+        )
 
     partition = (RowPartition(*partition[0]), ColPartition(*partition[1]))
     return partition
diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py
index edb4f370424..6ddb1c546d0 100644
--- a/pyomo/contrib/incidence_analysis/interface.py
+++ b/pyomo/contrib/incidence_analysis/interface.py
@@ -26,15 +26,17 @@
     get_independent_submatrices,
 )
 from pyomo.contrib.incidence_analysis.triangularize import (
+    get_scc_of_projection,
     block_triangularize,
     get_diagonal_blocks,
     get_blocks_from_maps,
-    )
+)
 from pyomo.contrib.incidence_analysis.dulmage_mendelsohn import (
     dulmage_mendelsohn,
     RowPartition,
     ColPartition,
-    )
+)
+
 if scipy_available:
     from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
     import scipy as sp
@@ -50,15 +52,27 @@ def _check_unindexed(complist):
     for comp in complist:
         if comp.is_indexed():
             raise ValueError(
-                    "Variables and constraints must be unindexed "
-                    "ComponentData objects. Got %s, which is indexed."
-                    % comp.name
-                    )
+                "Variables and constraints must be unindexed "
+                "ComponentData objects. Got %s, which is indexed." % comp.name
+            )
 
 
 def get_incidence_graph(variables, constraints, include_fixed=True):
-    """
-    This function gets the incidence graph of Pyomo variables and constraints.
+    return get_bipartite_incidence_graph(
+        variables, constraints, include_fixed=include_fixed
+    )
+
+
+def get_bipartite_incidence_graph(variables, constraints, include_fixed=True):
+    """Return the bipartite incidence graph of Pyomo variables and constraints.
+
+    Each node in the returned graph is an integer. The convention is that,
+    for a graph with N variables and M constraints, nodes 0 through M-1
+    correspond to constraints and nodes M through M+N-1 correspond to variables.
+    Nodes correspond to variables and constraints in the provided orders.
+    For consistency with NetworkX's "convention", constraint nodes are tagged
+    with `bipartite=0` while variable nodes are tagged with `bipartite=1`,
+    although these attributes are not used.
 
     Arguments:
     ----------
@@ -72,14 +86,15 @@ def get_incidence_graph(variables, constraints, include_fixed=True):
     Returns:
     --------
     NetworkX Graph
-        
+
     """
-    _check_unindexed(variables+constraints)
-    N, M = len(variables), len(constraints)
+    _check_unindexed(variables + constraints)
+    N = len(variables)
+    M = len(constraints)
     graph = nx.Graph()
     graph.add_nodes_from(range(M), bipartite=0)
-    graph.add_nodes_from(range(M, M+N), bipartite=1)
-    var_node_map = ComponentMap((v, M+i) for i, v in enumerate(variables))
+    graph.add_nodes_from(range(M, M + N), bipartite=1)
+    var_node_map = ComponentMap((v, M + i) for i, v in enumerate(variables))
     for i, con in enumerate(constraints):
         for var in identify_variables(con.expr, include_fixed=include_fixed):
             if var in var_node_map:
@@ -87,6 +102,62 @@ def get_incidence_graph(variables, constraints, include_fixed=True):
     return graph
 
 
+def extract_bipartite_subgraph(graph, nodes0, nodes1):
+    """Return the bipartite subgraph of a graph.
+
+    Two lists of nodes to project onto must be provided. These will correspond
+    to the "bipartite sets" in the subgraph. If the two sets provided have
+    M and N nodes, the subgraph will have nodes 0 through M+N-1, with the first
+    M corresponding to the first set provided and the last N corresponding
+    to the second set.
+
+    Parameters
+    ----------
+    graph: NetworkX Graph
+        The graph from which a subgraph is extracted
+    nodes0: list
+        A list of nodes in the original graph that will form the first
+        bipartite set of the projected graph (and have ``bipartite=0``)
+    nodes1: list
+        A list of nodes in the original graph that will form the second
+        bipartite set of the projected graph (and have ``bipartite=1``)
+
+    Returns
+    -------
+    subgraph: NetworkX Graph
+        Graph containing integer nodes corresponding to positions in the
+        provided lists, with edges where corresponding nodes are adjacent
+        in the original graph.
+
+    """
+    subgraph = nx.Graph()
+    sub_M = len(nodes0)
+    sub_N = len(nodes1)
+    subgraph.add_nodes_from(range(sub_M), bipartite=0)
+    subgraph.add_nodes_from(range(sub_M, sub_M + sub_N), bipartite=1)
+
+    old_new_map = {}
+    for i, node in enumerate(nodes0 + nodes1):
+        if node in old_new_map:
+            raise RuntimeError("Node %s provided more than once.")
+        old_new_map[node] = i
+
+    for (node1, node2) in graph.edges():
+        if node1 in old_new_map and node2 in old_new_map:
+            new_node_1 = old_new_map[node1]
+            new_node_2 = old_new_map[node2]
+            if (
+                subgraph.nodes[new_node_1]["bipartite"]
+                == subgraph.nodes[new_node_2]["bipartite"]
+            ):
+                raise RuntimeError(
+                    "Subgraph is not bipartite. Found an edge between nodes"
+                    " %s and %s (in the original graph)." % (node1, node2)
+                )
+            subgraph.add_edge(new_node_1, new_node_2)
+    return subgraph
+
+
 def _generate_variables_in_constraints(constraints, include_fixed=False):
     known_vars = ComponentSet()
     for con in constraints:
@@ -115,19 +186,21 @@ def get_structural_incidence_matrix(variables, constraints, include_fixed=True):
     Entries are 1.0.
 
     """
-    _check_unindexed(variables+constraints)
+    _check_unindexed(variables + constraints)
     N, M = len(variables), len(constraints)
     var_idx_map = ComponentMap((v, i) for i, v in enumerate(variables))
     rows = []
     cols = []
     for i, con in enumerate(constraints):
-        cols.extend(var_idx_map[v] for v in
-                identify_variables(con.expr, include_fixed=include_fixed)
-                if v in var_idx_map)
-        rows.extend([i]*(len(cols) - len(rows)))
+        cols.extend(
+            var_idx_map[v]
+            for v in identify_variables(con.expr, include_fixed=include_fixed)
+            if v in var_idx_map
+        )
+        rows.extend([i] * (len(cols) - len(rows)))
     assert len(rows) == len(cols)
-    data = [1.0]*len(rows)
-    matrix = sp.sparse.coo_matrix( (data, (rows, cols)), shape=(M, N) )
+    data = [1.0] * len(rows)
+    matrix = sp.sparse.coo_matrix((data, (rows, cols)), shape=(M, N))
     return matrix
 
 
@@ -155,21 +228,20 @@ class IncidenceGraphInterface(object):
     """
 
     def __init__(
-            self,
-            model=None,
-            active=True,
-            include_fixed=False,
-            include_inequality=True,
-            ):
-        """
-        """
+        self,
+        model=None,
+        active=True,
+        include_fixed=False,
+        include_inequality=True,
+    ):
+        """ """
         # If the user gives us a model or an NLP, we assume they want us
-        # to cache the incidence matrix for fast analysis of submatrices
-        # later on.
+        # to cache the incidence graph for fast analysis later on.
         # WARNING: This cache will become invalid if the user alters their
-        #          model.
+        # model.
         if model is None:
             self.cached = IncidenceMatrixType.NONE
+            self.incidence_graph = None
         elif isinstance(model, PyomoNLP):
             if not active:
                 raise ValueError(
@@ -186,20 +258,28 @@ def __init__(
             nlp = model
             self.cached = IncidenceMatrixType.NUMERIC
             self.variables = nlp.get_pyomo_variables()
-            self.constraints = nlp.get_pyomo_constraints()
+            self.constraints = [
+                con for con in nlp.get_pyomo_constraints()
+                if include_inequality or isinstance(con.expr, EqualityExpression)
+            ]
             self.var_index_map = ComponentMap(
-                    (var, idx) for idx, var in enumerate(self.variables))
+                (var, idx) for idx, var in enumerate(self.variables)
+            )
             self.con_index_map = ComponentMap(
-                    (con, idx) for idx, con in enumerate(self.constraints))
+                (con, idx) for idx, con in enumerate(self.constraints)
+            )
             if include_inequality:
-                self.incidence_matrix = nlp.evaluate_jacobian()
+                incidence_matrix = nlp.evaluate_jacobian()
             else:
-                self.incidence_matrix = nlp.evaluate_jacobian_eq()
+                incidence_matrix = nlp.evaluate_jacobian_eq()
+            nxb = nx.algorithms.bipartite
+            self.incidence_graph = nxb.from_biadjacency_matrix(incidence_matrix)
         elif isinstance(model, Block):
             self.cached = IncidenceMatrixType.STRUCTURAL
             self.constraints = [
-                con for con in
-                model.component_data_objects(Constraint, active=active)
+                con for con in model.component_data_objects(
+                    Constraint, active=active
+                )
                 if include_inequality or isinstance(con.expr, EqualityExpression)
             ]
             self.variables = list(
@@ -208,69 +288,124 @@ def __init__(
                 )
             )
             self.var_index_map = ComponentMap(
-                    (var, i) for i, var in enumerate(self.variables))
+                (var, i) for i, var in enumerate(self.variables)
+            )
             self.con_index_map = ComponentMap(
-                    (con, i) for i, con in enumerate(self.constraints))
-            self.incidence_matrix = get_structural_incidence_matrix(
-                    self.variables,
-                    self.constraints,
-                    )
+                (con, i) for i, con in enumerate(self.constraints)
+            )
+            self.incidence_graph = get_bipartite_incidence_graph(
+                self.variables,
+                self.constraints,
+                # Note that include_fixed is not necessary here. We have
+                # already checked this condition above.
+            )
         else:
             raise TypeError(
                 "Unsupported type for incidence matrix. Expected "
-                "%s or %s but got %s."
-                % (PyomoNLP, Block, type(model))
-                )
+                "%s or %s but got %s." % (PyomoNLP, Block, type(model))
+            )
 
+        # Note that these attributes are no longer used in get_diagonal_blocks,
+        # and should be removed. However, they are tested, so we will leave
+        # them in until a breaking change.
         self.row_block_map = None
         self.col_block_map = None
 
     def _validate_input(self, variables, constraints):
         if variables is None:
-            if self.cached is IncidenceMatrixType.NONE:
-                raise ValueError(
-                        "Neither variables nor a model have been provided."
-                        )
+            if self.incidence_graph is None:
+                raise ValueError("Neither variables nor a model have been provided.")
             else:
                 variables = self.variables
         if constraints is None:
-            if self.cached is IncidenceMatrixType.NONE:
-                raise ValueError(
-                        "Neither constraints nor a model have been provided."
-                        )
+            if self.incidence_graph is None:
+                raise ValueError("Neither constraints nor a model have been provided.")
             else:
                 constraints = self.constraints
 
-        _check_unindexed(variables+constraints)
+        _check_unindexed(variables + constraints)
         return variables, constraints
 
-    def _extract_submatrix(self, variables, constraints):
-        # Assumes variables and constraints are valid
-        if self.cached is IncidenceMatrixType.NONE:
-            return get_structural_incidence_matrix(
-                    variables,
-                    constraints,
-                    include_fixed=False,
-                    )
+    def _extract_subgraph(self, variables, constraints):
+        if self.incidence_graph is None:
+            # Note that, as variables are explicitly specified, there
+            # is no need for an include_fixed argument.
+            return get_bipartite_incidence_graph(variables, constraints)
+        else:
+            constraint_nodes = [self.con_index_map[con] for con in constraints]
+
+            # Note that this is the number of constraints in the original graph,
+            # not the subgraph.
+            M = len(self.constraints)
+            variable_nodes = [M + self.var_index_map[var] for var in variables]
+            subgraph = extract_bipartite_subgraph(
+                self.incidence_graph, constraint_nodes, variable_nodes
+            )
+            return subgraph
+
+    @property
+    def incidence_matrix(self):
+        if self.incidence_graph is None:
+            return None
         else:
-            N, M = len(variables), len(constraints)
-            old_new_var_indices = dict((self.var_index_map[v], i)
-                    for i, v in enumerate(variables))
-            old_new_con_indices = dict((self.con_index_map[c], i)
-                    for i, c in enumerate(constraints))
-            coo = self.incidence_matrix
-            new_row = []
-            new_col = []
-            new_data = []
-            for r, c, e in zip(coo.row, coo.col, coo.data):
-                if r in old_new_con_indices and c in old_new_var_indices:
-                    new_row.append(old_new_con_indices[r])
-                    new_col.append(old_new_var_indices[c])
-                    new_data.append(e)
+            M = len(self.constraints)
+            N = len(self.variables)
+            row = []
+            col = []
+            data = []
+            # Here we assume that the incidence graph is bipartite with nodes
+            # 0 through M-1 forming one of the bipartite sets.
+            for i in range(M):
+                assert self.incidence_graph.nodes[i]["bipartite"] == 0
+                for j in self.incidence_graph[i]:
+                    assert self.incidence_graph.nodes[j]["bipartite"] == 1
+                    row.append(i)
+                    col.append(j-M)
+                    data.append(1.0)
             return sp.sparse.coo_matrix(
-                    (new_data, (new_row, new_col)),
-                    shape=(M, N),
-                    )
+                (data, (row, col)),
+                shape=(M, N),
+            )
+
+    def get_adjacent_to(self, component):
+        """Return a list of components adjacent to the provided component
+        in the cached bipartite incidence graph of variables and constraints
+
+        Parameters
+        ----------
+        component: ComponentData
+            The variable or constraint data object whose adjacent components
+            are returned
+
+        Returns
+        -------
+        list of ComponentData
+            List of constraint or variable data objects adjacent to the
+            provided component
+
+        """
+        if self.incidence_graph is None:
+            raise RuntimeError(
+                "Cannot get components adjacent to %s if an incidence graph"
+                " is not cached." % component
+            )
+        _check_unindexed([component])
+        M = len(self.constraints)
+        N = len(self.variables)
+        if component in self.var_index_map:
+            vnode = M + self.var_index_map[component]
+            adj = self.incidence_graph[vnode]
+            adj_comps = [self.constraints[i] for i in adj]
+        elif component in self.con_index_map:
+            cnode = self.con_index_map[component]
+            adj = self.incidence_graph[cnode]
+            adj_comps = [self.variables[j-M] for j in adj]
+        else:
+            raise RuntimeError(
+                "Cannot find component %s in the cached incidence graph."
+                % component
+            )
+        return adj_comps
 
     def maximum_matching(self, variables=None, constraints=None):
         """
@@ -278,13 +413,15 @@ def maximum_matching(self, variables=None, constraints=None):
         in terms of a map from constraints to variables.
         """
         variables, constraints = self._validate_input(variables, constraints)
-        matrix = self._extract_submatrix(variables, constraints)
-
-        matching = maximum_matching(matrix.tocoo())
-        # Matching maps row (constraint) indices to column (variable) indices
-
-        return ComponentMap((constraints[i], variables[j])
-                for i, j in matching.items())
+        graph = self._extract_subgraph(variables, constraints)
+        con_nodes = list(range(len(constraints)))
+        matching = maximum_matching(graph, top_nodes=con_nodes)
+        # Matching maps constraint nodes to variable nodes. Here we need to
+        # know the convention according to which the graph was constructed.
+        M = len(constraints)
+        return ComponentMap(
+            (constraints[i], variables[j-M]) for i, j in matching.items()
+        )
 
     def get_connected_components(self, variables=None, constraints=None):
         """
@@ -293,13 +430,21 @@ def get_connected_components(self, variables=None, constraints=None):
         and constraints.
         """
         variables, constraints = self._validate_input(variables, constraints)
-        matrix = self._extract_submatrix(variables, constraints)
+        graph = self._extract_subgraph(variables, constraints)
+        nxc = nx.algorithms.components
+        M = len(constraints)
+        N = len(variables)
+        connected_components = list(nxc.connected_components(graph))
+
+        con_blocks = [
+            sorted([i for i in comp if i < M]) for comp in connected_components
+        ]
+        con_blocks = [[constraints[i] for i in block] for block in con_blocks]
+        var_blocks = [
+            sorted([j for j in comp if j >= M]) for comp in connected_components
+        ]
+        var_blocks = [[variables[i-M] for i in block] for block in var_blocks]
 
-        row_blocks, col_blocks = get_independent_submatrices(matrix.tocoo())
-        con_blocks = [[constraints[i] for i in block] for block in row_blocks]
-        var_blocks = [[variables[j] for j in block] for block in col_blocks]
-        # Switch the order of the partitions here to match the method call.
-        # Hopefully this does not get too confusing...
         return var_blocks, con_blocks
 
     def block_triangularize(self, variables=None, constraints=None):
@@ -310,21 +455,30 @@ def block_triangularize(self, variables=None, constraints=None):
         of the incidence matrix.
         """
         variables, constraints = self._validate_input(variables, constraints)
-        matrix = self._extract_submatrix(variables, constraints)
+        graph = self._extract_subgraph(variables, constraints)
 
-        row_block_map, col_block_map = block_triangularize(matrix.tocoo())
+        M = len(constraints)
+        con_nodes = list(range(M))
+        sccs = get_scc_of_projection(graph, con_nodes)
+        row_idx_map = {r: idx for idx, scc in enumerate(sccs) for r, _ in scc}
+        col_idx_map = {c-M: idx for idx, scc in enumerate(sccs) for _, c in scc}
         # Cache maps in case we want to get diagonal blocks quickly in the
         # future.
+        row_block_map = row_idx_map
+        col_block_map = {j-M: idx for j, idx in col_idx_map.items()}
         self.row_block_map = row_block_map
         self.col_block_map = col_block_map
-        con_block_map = ComponentMap((constraints[i], idx)
-                for i, idx in row_block_map.items())
-        var_block_map = ComponentMap((variables[j], idx)
-                for j, idx in col_block_map.items())
+        con_block_map = ComponentMap(
+            (constraints[i], idx) for i, idx in row_block_map.items()
+        )
+        var_block_map = ComponentMap(
+            (variables[j], idx) for j, idx in col_block_map.items()
+        )
         # Switch the order of the maps here to match the method call.
         # Hopefully this does not get too confusing...
         return var_block_map, con_block_map
 
+    # TODO: Update this method with a new name and deprecate old name.
     def get_diagonal_blocks(self, variables=None, constraints=None):
         """
         Returns the diagonal blocks in a block triangularization of the
@@ -339,16 +493,12 @@ def get_diagonal_blocks(self, variables=None, constraints=None):
 
         """
         variables, constraints = self._validate_input(variables, constraints)
-        matrix = self._extract_submatrix(variables, constraints)
-
-        if self.row_block_map is None or self.col_block_map is None:
-            block_rows, block_cols = get_diagonal_blocks(matrix)
-        else:
-            block_rows, block_cols = get_blocks_from_maps(
-                self.row_block_map, self.col_block_map
-            )
-        block_cons = [[constraints[i] for i in block] for block in block_rows]
-        block_vars = [[variables[i] for i in block] for block in block_cols]
+        graph = self._extract_subgraph(variables, constraints)
+        M = len(constraints)
+        con_nodes = list(range(M))
+        sccs = get_scc_of_projection(graph, con_nodes)
+        block_cons = [[constraints[i] for i, _ in scc] for scc in sccs]
+        block_vars = [[variables[j-M] for _, j in scc] for scc in sccs]
         return block_vars, block_cons
 
     def dulmage_mendelsohn(self, variables=None, constraints=None):
@@ -366,15 +516,18 @@ def dulmage_mendelsohn(self, variables=None, constraints=None):
 
         """
         variables, constraints = self._validate_input(variables, constraints)
-        matrix = self._extract_submatrix(variables, constraints)
-
-        row_partition, col_partition = dulmage_mendelsohn(matrix.tocoo())
+        graph = self._extract_subgraph(variables, constraints)
+        M = len(constraints)
+        top_nodes = list(range(M))
+        row_partition, col_partition = dulmage_mendelsohn(
+            graph, top_nodes=top_nodes
+        )
         con_partition = RowPartition(
-                *[[constraints[i] for i in subset] for subset in row_partition]
-                )
+            *[[constraints[i] for i in subset] for subset in row_partition]
+        )
         var_partition = ColPartition(
-                *[[variables[i] for i in subset] for subset in col_partition]
-                )
+            *[[variables[i-M] for i in subset] for subset in col_partition]
+        )
         # Switch the order of the maps here to match the method call.
         # Hopefully this does not get too confusing...
         return var_partition, con_partition
@@ -402,7 +555,7 @@ def remove_nodes(self, nodes, constraints=None):
         """
         if constraints is None:
             constraints = []
-        if self.cached is IncidenceMatrixType.NONE:
+        if self.incidence_graph is None:
             raise RuntimeError(
                 "Attempting to remove variables and constraints from cached "
                 "incidence matrix,\nbut no incidence matrix has been cached."
@@ -411,13 +564,11 @@ def remove_nodes(self, nodes, constraints=None):
         to_exclude.update(constraints)
         vars_to_include = [v for v in self.variables if v not in to_exclude]
         cons_to_include = [c for c in self.constraints if c not in to_exclude]
-        incidence_matrix = self._extract_submatrix(
-            vars_to_include, cons_to_include
-        )
+        incidence_graph = self._extract_subgraph(vars_to_include, cons_to_include)
         # update attributes
         self.variables = vars_to_include
         self.constraints = cons_to_include
-        self.incidence_matrix = incidence_matrix
+        self.incidence_graph = incidence_graph
         self.var_index_map = ComponentMap(
             (var, i) for i, var in enumerate(self.variables)
         )
diff --git a/pyomo/contrib/incidence_analysis/matching.py b/pyomo/contrib/incidence_analysis/matching.py
index 25906958103..e9ba4e3fb00 100644
--- a/pyomo/contrib/incidence_analysis/matching.py
+++ b/pyomo/contrib/incidence_analysis/matching.py
@@ -11,34 +11,78 @@
 
 from pyomo.common.dependencies import networkx as nx
 
-def maximum_matching(matrix):
-    """
-    Returns a maximum matching of the rows and columns of the 
-    matrix as a dict from row indices to column indices.
+
+def maximum_matching(matrix_or_graph, top_nodes=None):
+    """Return a maximum cardinality matching of the provided matrix or
+    bipartite graph
+
+    If a matrix is provided, the matching is returned as a map from row
+    indices to column indices. If a bipartite graph is provided, a list of
+    "top nodes" must be provided as well. These correspond to one of the
+    "bipartite sets". The matching is then returned as a map from "top nodes"
+    to the other set of nodes.
+
+    Parameters
+    ----------
+    matrix_or_graph: SciPy sparse matrix or NetworkX Graph
+        The matrix or graph whose maximum matching will be computed
+    top_nodes: list
+        Integer nodes representing a bipartite set in a graph. Must be provided
+        if and only if a NetworkX Graph is provided.
+
+    Returns
+    -------
+    max_matching: dict
+        Dict mapping from integer nodes in the first bipartite set (row
+        indices) to nodes in the second (column indices).
+
     """
     nxb = nx.algorithms.bipartite
     nxc = nx.algorithms.components
     from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix
 
-    M, N = matrix.shape
-    bg = from_biadjacency_matrix(matrix)
-
-    # Check assumptions regarding from_biadjacency_matrix function:
-    for i in range(M):
-        # First M nodes in graph correspond to rows
-        assert bg.nodes[i]['bipartite'] == 0
-
-    for j in range(M, M+N):
-        # Last N nodes in graph correspond to columns
-        assert bg.nodes[j]['bipartite'] == 1
-
-    # If the matrix is block diagonal, the graph will be disconnected.
-    # This is fine, but we need to separate into connected components
-    # for NetworkX to not complain.
-    conn_comp = [bg.subgraph(c) for c in nxc.connected_components(bg)]
-    matchings = [nxb.maximum_matching(c) for c in conn_comp]
-    # If n0 < M, then n1 >= M. n0 is the row index, n1-M is the column index
-    max_matching = {
-            n0: n1-M for m in matchings for n0, n1 in m.items() if n0 < M
-            }
+    if isinstance(matrix_or_graph, nx.Graph):
+        graph_provided = True
+        if top_nodes is None:
+            raise RuntimeError(
+                "top_nodes argument must be set if a graph is provided."
+            )
+        M = len(top_nodes)
+        N = len(matrix_or_graph.nodes) - M
+        bg = matrix_or_graph
+        if not nxb.is_bipartite(bg):
+            raise RuntimeError("Provided graph is not bipartite.")
+    else:
+        graph_provided = False
+        # Assume something SciPy-sparse compatible was provided.
+        if top_nodes is not None:
+            raise RuntimeError(
+                "top_nodes argument cannot be used if a matrix is provided"
+            )
+        M, N = matrix_or_graph.shape
+        top_nodes = list(range(M))
+        bg = from_biadjacency_matrix(matrix_or_graph)
+
+        # Check assumptions regarding from_biadjacency_matrix function:
+        for i in range(M):
+            # First M nodes in graph correspond to rows
+            assert bg.nodes[i]["bipartite"] == 0
+
+        for j in range(M, M + N):
+            # Last N nodes in graph correspond to columns
+            assert bg.nodes[j]["bipartite"] == 1
+
+    matching = nxb.maximum_matching(bg, top_nodes=top_nodes)
+    if graph_provided:
+        top_node_set = set(top_nodes)
+        # If a graph was provided, we return a mapping from "top nodes"
+        # to their matched "bottom nodes"
+        max_matching = {
+            n0: n1 for n0, n1 in matching.items() if n0 in top_node_set
+        }
+    else:
+        # If a matrix was provided, we return a mapping from row indices
+        # to column indices.
+        max_matching = {n0: n1-M for n0, n1 in matching.items() if n0 < M}
+
     return max_matching
diff --git a/pyomo/contrib/incidence_analysis/tests/test_interface.py b/pyomo/contrib/incidence_analysis/tests/test_interface.py
index c4f4a1e0e8a..4668e6344a7 100644
--- a/pyomo/contrib/incidence_analysis/tests/test_interface.py
+++ b/pyomo/contrib/incidence_analysis/tests/test_interface.py
@@ -18,6 +18,8 @@
     get_structural_incidence_matrix,
     get_numeric_incidence_matrix,
     get_incidence_graph,
+    get_bipartite_incidence_graph,
+    extract_bipartite_subgraph,
 )
 from pyomo.contrib.incidence_analysis.matching import maximum_matching
 from pyomo.contrib.incidence_analysis.triangularize import block_triangularize
@@ -32,6 +34,7 @@
 if scipy_available:
     from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
 if networkx_available:
+    import networkx as nx
     from networkx.algorithms.bipartite.matrix import from_biadjacency_matrix
 from pyomo.contrib.pynumero.asl import AmplInterface
 
@@ -1238,5 +1241,239 @@ def test_include_inequality_nlp(self):
         self.assertEqual(igraph.incidence_matrix.shape, (12, 8))
 
 
+@unittest.skipUnless(networkx_available, "networkx is not available.")
+@unittest.skipUnless(scipy_available, "scipy is not available.")
+class TestGetIncidenceGraph(unittest.TestCase):
+
+    def make_test_model(self):
+        m = pyo.ConcreteModel()
+        m.I = pyo.Set(initialize=[1, 2, 3, 4])
+        m.v = pyo.Var(m.I, bounds=(0, None))
+        m.eq1 = pyo.Constraint(expr=m.v[1]**2 + m.v[2]**2 == 1.0)
+        m.eq2 = pyo.Constraint(expr=m.v[1] + 2.0 == m.v[3])
+        m.ineq1 = pyo.Constraint(expr=m.v[2] - m.v[3]**0.5 + m.v[4]**2 <= 1.0)
+        m.ineq2 = pyo.Constraint(expr=m.v[2]*m.v[4] >= 1.0)
+        m.ineq3 = pyo.Constraint(expr=m.v[1] >= m.v[4]**4)
+        m.obj = pyo.Objective(expr=-m.v[1] - m.v[2] + m.v[3]**2 + m.v[4]**2)
+        return m
+
+    def test_bipartite_incidence_graph(self):
+        m = self.make_test_model()
+        constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3]
+        variables = list(m.v.values())
+        graph = get_bipartite_incidence_graph(variables, constraints)
+
+        # Nodes:
+        #   0: m.eq1
+        #   1: m.eq2
+        #   2: m.ineq1
+        #   3: m.ineq2
+        #   4: m.ineq3
+        #   5: m.v[1]
+        #   6: m.v[2]
+        #   7: m.v[3]
+        #   8: m.v[4]
+
+        # Assert some basic structure
+        self.assertEqual(len(graph.nodes), 9)
+        self.assertEqual(len(graph.edges), 11)
+        self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph))
+
+        # Assert that the "adjacency list" is what we expect
+        self.assertEqual(set(graph[0]), {5, 6})
+        self.assertEqual(set(graph[1]), {5, 7})
+        self.assertEqual(set(graph[2]), {6, 7, 8})
+        self.assertEqual(set(graph[3]), {6, 8})
+        self.assertEqual(set(graph[4]), {5, 8})
+        self.assertEqual(set(graph[5]), {0, 1, 4})
+        self.assertEqual(set(graph[6]), {0, 2, 3})
+        self.assertEqual(set(graph[7]), {1, 2})
+        self.assertEqual(set(graph[8]), {2, 3, 4})
+
+    def test_unused_var(self):
+        m = self.make_test_model()
+        constraints = [m.eq1, m.eq2]
+        variables = list(m.v.values())
+        graph = get_bipartite_incidence_graph(variables, constraints)
+
+        # Nodes:
+        #   0: m.eq1
+        #   1: m.eq2
+        #   2: m.v[1]
+        #   3: m.v[2]
+        #   4: m.v[3]
+        #   5: m.v[4]
+
+        self.assertEqual(len(graph.nodes), 6)
+        self.assertEqual(len(graph.edges), 4)
+        self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph))
+
+        # Assert that the "adjacency list" is what we expect
+        self.assertEqual(set(graph[0]), {2, 3})
+        self.assertEqual(set(graph[1]), {2, 4})
+        self.assertEqual(set(graph[2]), {0, 1})
+        self.assertEqual(set(graph[3]), {0})
+        self.assertEqual(set(graph[4]), {1})
+        self.assertEqual(set(graph[5]), set())
+
+    def test_fixed_vars(self):
+        m = self.make_test_model()
+        constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3]
+        variables = list(m.v.values())
+        m.v[1].fix()
+        m.v[4].fix()
+
+        # Slightly odd situation where we provide fixed variables, but
+        # then tell the graph to not include them. Nodes will be created
+        # for these vars, but they will not have any edges.
+        graph = get_bipartite_incidence_graph(
+            variables, constraints, include_fixed=False
+        )
+
+        # Nodes:
+        #   0: m.eq1
+        #   1: m.eq2
+        #   2: m.ineq1
+        #   3: m.ineq2
+        #   4: m.ineq3
+        #   5: m.v[1]
+        #   6: m.v[2]
+        #   7: m.v[3]
+        #   8: m.v[4]
+
+        self.assertEqual(len(graph.nodes), 9)
+        self.assertEqual(len(graph.edges), 5)
+        self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph))
+
+        # Assert that the "adjacency list" is what we expect
+        self.assertEqual(set(graph[0]), {6})
+        self.assertEqual(set(graph[1]), {7})
+        self.assertEqual(set(graph[2]), {6, 7})
+        self.assertEqual(set(graph[3]), {6})
+        self.assertEqual(set(graph[4]), set())
+        self.assertEqual(set(graph[5]), set())
+        self.assertEqual(set(graph[6]), {0, 2, 3})
+        self.assertEqual(set(graph[7]), {1, 2})
+        self.assertEqual(set(graph[8]), set())
+
+    def test_extract_subgraph(self):
+        m = self.make_test_model()
+        constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3]
+        variables = list(m.v.values())
+        graph = get_bipartite_incidence_graph(variables, constraints)
+
+        sg_cons = [0, 2]
+        sg_vars = [i + len(constraints) for i in [2, 0, 3]]
+
+        subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars)
+
+        # Subgraph nodes:
+        #   0: m.eq1
+        #   1: m.ineq1
+        #   2: m.v[3]
+        #   3: m.v[1]
+        #   4: m.v[4]
+
+        self.assertEqual(len(subgraph.nodes), 5)
+        self.assertEqual(len(subgraph.edges), 3)
+        self.assertTrue(nx.algorithms.bipartite.is_bipartite(subgraph))
+
+        self.assertEqual(set(subgraph[0]), {3})
+        self.assertEqual(set(subgraph[1]), {2, 4})
+        self.assertEqual(set(subgraph[2]), {1})
+        self.assertEqual(set(subgraph[3]), {0})
+        self.assertEqual(set(subgraph[4]), {1})
+
+    def test_extract_exceptions(self):
+        m = self.make_test_model()
+        constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3]
+        variables = list(m.v.values())
+        graph = get_bipartite_incidence_graph(variables, constraints)
+
+        sg_cons = [0, 2, 5]
+        sg_vars = [i + len(constraints) for i in [2, 3]]
+        msg = "Subgraph is not bipartite"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars)
+
+        sg_cons = [0, 2, 5]
+        sg_vars = [i + len(constraints) for i in [2, 0, 3]]
+        msg = "provided more than once"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars)
+
+
+@unittest.skipUnless(networkx_available, "networkx is not available.")
+@unittest.skipUnless(scipy_available, "scipy is not available.")
+class TestGetAdjacent(unittest.TestCase):
+
+    def test_get_adjacent_to_var(self):
+        m = make_degenerate_solid_phase_model()
+        igraph = IncidenceGraphInterface(m)
+        adj_cons = igraph.get_adjacent_to(m.rho)
+        self.assertEqual(
+            ComponentSet(adj_cons),
+            ComponentSet([
+                m.holdup_eqn[1],
+                m.holdup_eqn[2],
+                m.holdup_eqn[3],
+                m.density_eqn,
+            ]),
+        )
+
+    def test_get_adjacent_to_con(self):
+        m = make_degenerate_solid_phase_model()
+        igraph = IncidenceGraphInterface(m)
+        adj_vars = igraph.get_adjacent_to(m.density_eqn)
+        self.assertEqual(
+            ComponentSet(adj_vars),
+            ComponentSet([
+                m.x[1],
+                m.x[2],
+                m.x[3],
+                m.rho,
+            ]),
+        )
+
+    def test_get_adjacent_exceptions(self):
+        m = make_degenerate_solid_phase_model()
+        igraph = IncidenceGraphInterface()
+        msg = "Cannot get components adjacent to"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            adj_vars = igraph.get_adjacent_to(m.density_eqn)
+
+        m.x[1].fix()
+        igraph = IncidenceGraphInterface(m, include_fixed=False)
+        msg = "Cannot find component"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            adj_cons = igraph.get_adjacent_to(m.x[1])
+
+
+@unittest.skipUnless(networkx_available, "networkx is not available.")
+@unittest.skipUnless(scipy_available, "scipy is not available.")
+class TestInterface(unittest.TestCase):
+
+    def test_subgraph_with_fewer_var_or_con(self):
+        m = pyo.ConcreteModel()
+        m.I = pyo.Set(initialize=[1, 2])
+        m.v = pyo.Var(m.I)
+        m.eq1 = pyo.Constraint(expr=m.v[1] + m.v[2] == 1)
+        m.ineq1 = pyo.Constraint(expr=m.v[1] - m.v[2] <= 2)
+
+        # Defensively set include_inequality=True, which is the current
+        # default, in case this default changes.
+        igraph = IncidenceGraphInterface(m, include_inequality=True)
+
+        variables = list(m.v.values())
+        constraints = [m.ineq1]
+        matching = igraph.maximum_matching(variables, constraints)
+        self.assertEqual(len(matching), 1)
+
+        variables = [m.v[2]]
+        constraints = [m.eq1, m.ineq1]
+        matching = igraph.maximum_matching(variables, constraints)
+        self.assertEqual(len(matching), 1)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/pyomo/contrib/incidence_analysis/tests/test_triangularize.py b/pyomo/contrib/incidence_analysis/tests/test_triangularize.py
index f7542fcb6dc..9165f26f462 100644
--- a/pyomo/contrib/incidence_analysis/tests/test_triangularize.py
+++ b/pyomo/contrib/incidence_analysis/tests/test_triangularize.py
@@ -12,20 +12,118 @@
 import random
 from pyomo.contrib.incidence_analysis.matching import maximum_matching
 from pyomo.contrib.incidence_analysis.triangularize import (
-        block_triangularize,
-        get_diagonal_blocks,
-        )
+    get_scc_of_projection,
+    block_triangularize,
+    get_diagonal_blocks,
+)
 from pyomo.common.dependencies import (
-        scipy,
-        scipy_available,
-        networkx_available,
-        )
+    scipy,
+    scipy_available,
+    networkx as nx,
+    networkx_available,
+)
 if scipy_available:
     sps = scipy.sparse
+if networkx_available:
+    nxb = nx.algorithms.bipartite
 
 import pyomo.common.unittest as unittest
 
 
+@unittest.skipUnless(networkx_available, "networkx is not available")
+@unittest.skipUnless(scipy_available, "scipy is not available")
+class TestGetSCCOfProjection(unittest.TestCase):
+
+    def test_graph_decomposable_tridiagonal_shuffled(self):
+        """
+        This is the same graph as in test_decomposable_tridiagonal_shuffled
+        below, but now we convert the matrix into a bipartite graph and
+        use get_scc_of_projection.
+
+        The matrix decomposes into 2x2 blocks:
+        |x x      |
+        |x x      |
+        |  x x x  |
+        |    x x  |
+        |      x x|
+        """
+        N = 11
+        row = []
+        col = []
+        data = []
+
+        # Diagonal
+        row.extend(range(N))
+        col.extend(range(N))
+        data.extend(1 for _ in range(N))
+
+        # Below diagonal
+        row.extend(range(1, N))
+        col.extend(range(N-1))
+        data.extend(1 for _ in range(N-1))
+
+        # Above diagonal
+        row.extend(i for i in range(N-1) if not i%2)
+        col.extend(i+1 for i in range(N-1) if not i%2)
+        data.extend(1 for i in range(N-1) if not i%2)
+
+        # Same results hold after applying a random permutation.
+        row_perm = list(range(N))
+        col_perm = list(range(N))
+        random.shuffle(row_perm)
+        random.shuffle(col_perm)
+
+        row = [row_perm[i] for i in row]
+        col = [col_perm[j] for j in col]
+
+        matrix = sps.coo_matrix((data, (row, col)), shape=(N, N))
+        graph = nxb.matrix.from_biadjacency_matrix(matrix)
+        row_nodes = list(range(N))
+        sccs = get_scc_of_projection(graph, row_nodes)
+
+        self.assertEqual(len(sccs), (N+1)//2)
+
+        for i in range((N+1)//2):
+            # Note that these rows and cols are in the permuted space
+            rows = set(r for r, _ in sccs[i])
+            cols = set(c-N for _, c in sccs[i])
+
+            pred_rows = {row_perm[2*i]}
+            pred_cols = {col_perm[2*i]}
+
+            if 2*i+1 < N:
+                pred_rows.add(row_perm[2*i+1])
+                pred_cols.add(col_perm[2*i+1])
+
+            self.assertEqual(pred_rows, rows)
+            self.assertEqual(pred_cols, cols)
+
+    def test_scc_exceptions(self):
+        graph = nx.Graph()
+        graph.add_nodes_from(range(3))
+        graph.add_edges_from([(0, 1), (0, 2), (1, 2)])
+        top_nodes = [0]
+        msg = "graph is not bipartite"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            sccs = get_scc_of_projection(graph, top_nodes=top_nodes)
+
+        graph = nx.Graph()
+        graph.add_nodes_from(range(3))
+        graph.add_edges_from([(0, 1), (0, 2)])
+        top_nodes[0]
+        msg = "bipartite sets of different cardinalities"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            sccs = get_scc_of_projection(graph, top_nodes=top_nodes)
+
+        graph = nx.Graph()
+        graph.add_nodes_from(range(4))
+        graph.add_edges_from([(0, 1), (0, 2)])
+        top_nodes = [0, 3]
+        msg = "without a perfect matching"
+        with self.assertRaisesRegex(RuntimeError, msg):
+            sccs = get_scc_of_projection(graph, top_nodes=top_nodes)
+
+
 @unittest.skipUnless(networkx_available, "networkx is not available")
 @unittest.skipUnless(scipy_available, "scipy is not available")
 class TestTriangularize(unittest.TestCase):
@@ -37,7 +135,7 @@ def test_low_rank_exception(self):
 
         matrix = sps.coo_matrix((data, (row, col)), shape=(N, N))
 
-        with self.assertRaises(ValueError) as exc:
+        with self.assertRaises(RuntimeError) as exc:
             row_block_map, col_block_map = block_triangularize(matrix)
         self.assertIn('perfect matching', str(exc.exception))
 
diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py
index 7c5d8865c52..8007cd8cb9b 100644
--- a/pyomo/contrib/incidence_analysis/triangularize.py
+++ b/pyomo/contrib/incidence_analysis/triangularize.py
@@ -10,94 +10,137 @@
 #  ___________________________________________________________________________
 
 from pyomo.contrib.incidence_analysis.matching import maximum_matching
+from pyomo.contrib.incidence_analysis.common.dulmage_mendelsohn import (
+    # TODO: The fact that we import this function here suggests it should be
+    # promoted.
+    _get_projected_digraph,
+)
 from pyomo.common.dependencies import networkx as nx
 
 
-def block_triangularize(matrix, matching=None):
-    """
-    Computes the necessary information to permute a matrix to block-lower
-    triangular form, i.e. a partition of rows and columns into an ordered
-    set of diagonal blocks in such a permutation.
+def get_scc_of_projection(graph, top_nodes, matching=None):
+    """Return the topologically ordered strongly connected components of a
+    bipartite graph, projected with respect to a perfect matching
 
-    Arguments
-    ---------
-    matrix: A SciPy sparse matrix
-    matching: A perfect matching of rows and columns, in the form of a dict
-              mapping row indices to column indices
+    The provided undirected bipartite graph is projected into a directed graph
+    on the set of "top nodes" by treating "matched edges" as out-edges and
+    "unmatched edges" as in-edges. Then the strongly connected components of
+    the directed graph are computed. These strongly connected components are
+    unique, regardless of the choice of perfect matching. The strongly connected
+    components form a directed acyclic graph, and are returned in a topological
+    order. The order is unique, as ambiguities are resolved "lexicographically".
+
+    The "direction" of the projection (where matched edges are out-edges)
+    leads to a block *lower* triangular permutation when the top nodes
+    correspond to *rows* in the bipartite graph of a matrix.
+
+    Parameters
+    ----------
+    graph: NetworkX Graph
+        A bipartite graph
+    top_nodes: list
+        One of the bipartite sets in the graph
+    matching: dict
+        Maps each node in ``top_nodes`` to its matched node
 
     Returns
     -------
-    Two dicts. The first maps each row index to the index of its block in a
-    block-lower triangular permutation of the matrix. The second maps each
-    column index to the index of its block in a block-lower triangular
-    permutation of the matrix.
+    list of lists
+        The outer list is a list of strongly connected components. Each
+        strongly connected component is a list of tuples of matched nodes.
+        The first node is a "top node", and the second is an "other node".
+
     """
     nxb = nx.algorithms.bipartite
     nxc = nx.algorithms.components
     nxd = nx.algorithms.dag
-    from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix
-
-    M, N = matrix.shape
+    if not nxb.is_bipartite(graph):
+        raise RuntimeError("Provided graph is not bipartite.")
+    M = len(top_nodes)
+    N = len(graph.nodes) - M
     if M != N:
-        raise ValueError("block_triangularize does not currently "
-           "support non-square matrices. Got matrix with shape %s."
-           % (matrix.shape,)
-           )
-    bg = from_biadjacency_matrix(matrix)
-
+        raise RuntimeError(
+            "get_scc_of_projection does not support bipartite graphs with"
+            " bipartite sets of different cardinalities. Got sizes %s and"
+            " %s." % (M, N)
+        )
     if matching is None:
-        matching = maximum_matching(matrix)
-
-    len_matching = len(matching)
-    if len_matching != M:
-        raise ValueError("block_triangularize only supports matrices "
-                "that have a perfect matching of rows and columns. "
-                "Cardinality of maximal matching is %s" % len_matching
-                )
-
-    # Construct directed graph of rows
-    dg = nx.DiGraph()
-    dg.add_nodes_from(range(M))
-    for n in dg.nodes:
-        col_idx = matching[n]
-        col_node = col_idx + M
-        # For all rows that share this column
-        for neighbor in bg[col_node]:
-            if neighbor != n:
-                # Add an edge towards this column's matched row
-                dg.add_edge(neighbor, n)
-
-    # Partition the rows into strongly connected components (diagonal blocks)
+        # This matching maps top nodes to "other nodes" *and* other nodes
+        # back to top nodes.
+        matching = nxb.maximum_matching(graph, top_nodes=top_nodes)
+    if len(matching) != 2*M:
+        raise RuntimeError(
+            "get_scc_of_projection does not support bipartite graphs without"
+            " a perfect matching. Got a graph with %s nodes per bipartite set"
+            " and a matching of cardinality %s." % (M, (len(matching)/2))
+        )
+
+    # _get_projected_digraph treats matched edges as "in-edges", so we
+    # reverse the direction of edges here.
+    dg = _get_projected_digraph(graph, matching, top_nodes).reverse()
+
     scc_list = list(nxc.strongly_connected_components(dg))
+    n_scc = len(scc_list)
     node_scc_map = {n: idx for idx, scc in enumerate(scc_list) for n in scc}
 
     # Now we need to put the SCCs in the right order. We do this by performing
     # a topological sort on the DAG of SCCs.
     dag = nx.DiGraph()
-    for i, c in enumerate(scc_list):
-        dag.add_node(i)
+    dag.add_nodes_from(range(n_scc))
     for n in dg.nodes:
         source_scc = node_scc_map[n]
         for neighbor in dg[n]:
             target_scc = node_scc_map[neighbor]
             if target_scc != source_scc:
-                dag.add_edge(target_scc, source_scc)
-                # Reverse direction of edge. This corresponds to creating
-                # a block lower triangular matrix.
+                dag.add_edge(source_scc, target_scc)
 
     scc_order = list(nxd.lexicographical_topological_sort(dag))
 
-    scc_block_map = {c: i for i, c in enumerate(scc_order)}
-    row_block_map = {n: scc_block_map[c] for n, c in node_scc_map.items()}
-    # ^ This maps row indices to the blocks they belong to.
+    # The "natural" return type, here, is a list of lists. Each inner list
+    # is an SCC, and contains tuples of nodes. The "top node", and its matched
+    # "bottom node".
+    ordered_node_subsets = [
+        sorted([(i, matching[i]) for i in scc_list[scc_idx]])
+        for scc_idx in scc_order
+    ]
+    return ordered_node_subsets
 
-    # Invert the matching to map row indices to column indices
-    col_row_map = {c: r for r, c in matching.items()}
-    assert len(col_row_map) == M
 
-    col_block_map = {c: row_block_map[col_row_map[c]] for c in range(N)}
+def block_triangularize(matrix, matching=None):
+    """
+    Computes the necessary information to permute a matrix to block-lower
+    triangular form, i.e. a partition of rows and columns into an ordered
+    set of diagonal blocks in such a permutation.
+
+    Arguments
+    ---------
+    matrix: A SciPy sparse matrix
+    matching: A perfect matching of rows and columns, in the form of a dict
+              mapping row indices to column indices
 
-    return row_block_map, col_block_map
+    Returns
+    -------
+    Two dicts. The first maps each row index to the index of its block in a
+    block-lower triangular permutation of the matrix. The second maps each
+    column index to the index of its block in a block-lower triangular
+    permutation of the matrix.
+    """
+    nxb = nx.algorithms.bipartite
+    nxc = nx.algorithms.components
+    nxd = nx.algorithms.dag
+    from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix
+    M, N = matrix.shape
+    if M != N:
+        raise ValueError(
+            "block_triangularize does not currently support non-square"
+            " matrices. Got matrix with shape %s." % ((M, N),)
+        )
+    graph = from_biadjacency_matrix(matrix)
+    row_nodes = list(range(M))
+    sccs = get_scc_of_projection(graph, row_nodes, matching=matching)
+    row_idx_map = {r: idx for idx, scc in enumerate(sccs) for r, _ in scc}
+    col_idx_map = {c-M: idx for idx, scc in enumerate(sccs) for _, c in scc}
+    return row_idx_map, col_idx_map
 
 
 def get_blocks_from_maps(row_block_map, col_block_map):
@@ -109,7 +152,7 @@ def get_blocks_from_maps(row_block_map, col_block_map):
     Arguments
     ---------
     row_block_map: dict
-        Dict mapping each row coordinate to the coordinate of the 
+        Dict mapping each row coordinate to the coordinate of the
         block it belongs to
 
     col_block_map: dict
@@ -160,8 +203,6 @@ def get_diagonal_blocks(matrix, matching=None):
         diagonal blocks.
 
     """
-    row_block_map, col_block_map = block_triangularize(
-        matrix, matching=matching
-    )
+    row_block_map, col_block_map = block_triangularize(matrix, matching=matching)
     block_rows, block_cols = get_blocks_from_maps(row_block_map, col_block_map)
     return block_rows, block_cols
diff --git a/pyomo/contrib/incidence_analysis/util.py b/pyomo/contrib/incidence_analysis/util.py
index 9e592dad922..0e5a7d05cf9 100644
--- a/pyomo/contrib/incidence_analysis/util.py
+++ b/pyomo/contrib/incidence_analysis/util.py
@@ -15,19 +15,19 @@
 from pyomo.core.expr.visitor import identify_variables
 from pyomo.util.calc_var_value import calculate_variable_from_constraint
 from pyomo.util.subsystems import (
-        create_subsystem_block,
-        TemporarySubsystemManager,
-        generate_subsystem_blocks,
-        )
+    create_subsystem_block,
+    TemporarySubsystemManager,
+    generate_subsystem_blocks,
+)
 from pyomo.contrib.incidence_analysis.interface import IncidenceGraphInterface
 
 
 def generate_strongly_connected_components(
-        constraints,
-        variables=None,
-        include_fixed=False,
-        ):
-    """ Performs a block triangularization of the incidence matrix
+    constraints,
+    variables=None,
+    include_fixed=False,
+):
+    """Performs a block triangularization of the incidence matrix
     of the provided constraints and variables, and yields a block that
     contains the constraints and variables of each diagonal block
     (strongly connected component).
@@ -55,9 +55,9 @@ def generate_strongly_connected_components(
         variables = []
         for con in constraints:
             for var in identify_variables(
-                    con.expr,
-                    include_fixed=include_fixed,
-                    ):
+                con.expr,
+                include_fixed=include_fixed,
+            ):
                 if var not in var_set:
                     variables.append(var)
                     var_set.add(var)
@@ -65,9 +65,9 @@ def generate_strongly_connected_components(
     assert len(variables) == len(constraints)
     igraph = IncidenceGraphInterface()
     var_block_map, con_block_map = igraph.block_triangularize(
-            variables=variables,
-            constraints=constraints,
-            )
+        variables=variables,
+        constraints=constraints,
+    )
     blocks = set(var_block_map.values())
     n_blocks = len(blocks)
     var_blocks = [[] for b in range(n_blocks)]
@@ -78,21 +78,21 @@ def generate_strongly_connected_components(
         con_blocks[b].append(con)
     subsets = list(zip(con_blocks, var_blocks))
     for block, inputs in generate_subsystem_blocks(
-            subsets,
-            include_fixed=include_fixed,
-            ):
+        subsets,
+        include_fixed=include_fixed,
+    ):
         # TODO: How does len scale for reference-to-list?
         assert len(block.vars) == len(block.cons)
         yield (block, inputs)
 
 
 def solve_strongly_connected_components(
-        block,
-        solver=None,
-        solve_kwds=None,
-        calc_var_kwds=None,
-        ):
-    """ This function solves a square block of variables and equality
+    block,
+    solver=None,
+    solve_kwds=None,
+    calc_var_kwds=None,
+):
+    """This function solves a square block of variables and equality
     constraints by solving strongly connected components individually.
     Strongly connected components (of the directed graph of constraints
     obtained from a perfect matching of variables and constraints) are
@@ -139,9 +139,9 @@ def solve_strongly_connected_components(
 
     res_list = []
     for scc, inputs in generate_strongly_connected_components(
-            constraints,
-            variables,
-            ):
+        constraints,
+        variables,
+    ):
         with TemporarySubsystemManager(to_fix=inputs):
             if len(scc.vars) == 1:
                 results = calculate_variable_from_constraint(
@@ -158,9 +158,8 @@ def solve_strongly_connected_components(
                     raise RuntimeError(
                         "An external solver is required if block has strongly\n"
                         "connected components of size greater than one (is not "
-                        "a DAG).\nGot an SCC with components: \n%s\n%s"
-                        % (vars, cons)
-                        )
+                        "a DAG).\nGot an SCC with components: \n%s\n%s" % (vars, cons)
+                    )
                 results = solver.solve(scc, **solve_kwds)
                 res_list.append(results)
     return res_list