diff --git a/src/MathProgIncidence.jl b/src/MathProgIncidence.jl
index 350fdf5..cdcdfac 100644
--- a/src/MathProgIncidence.jl
+++ b/src/MathProgIncidence.jl
@@ -34,6 +34,8 @@ include("dulmage_mendelsohn.jl")
 include("block_triangularize.jl")
 include("interface.jl")
 
+include("visualize.jl")
+
 # Methods to get incidence matrices (as SparseMatrixCSC) from
 # JuMP models or incidence graphs
 include("incidence_matrix.jl")
diff --git a/src/interface.jl b/src/interface.jl
index 2b4662f..764d50e 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -277,6 +277,8 @@ function maximum_matching(
     return maximum_matching(igraph)
 end
 
+# TODO: These types should be parameterized to handle integers as well
+# as variables/constraints.
 const DMConPartition = NamedTuple{
     (:underconstrained, :square, :overconstrained, :unmatched),
     Tuple{
@@ -297,6 +299,15 @@ const DMVarPartition = NamedTuple{
     },
 }
 
+const DulmageMendelsohnDecomposition = NamedTuple{
+    (:con, :var),
+    # NOTE: These fields are not fully typed as we need to support ints as
+    # wells as variables/constraints.
+    # TODO: Parameterize this type? Or just always use ints, then keep the vars/cons
+    # somewhere else?
+    Tuple{NamedTuple,NamedTuple},
+}
+
 """
     dulmage_mendelsohn(igraph::IncidenceGraphInterface)
 
@@ -394,7 +405,7 @@ julia> # As there are no unmatched constraints, the overconstrained subsystem is
 """
 function dulmage_mendelsohn(
     igraph::IncidenceGraphInterface
-)
+)::DulmageMendelsohnDecomposition
     ncon = length(igraph._con_node_map)
     con_node_set = Set(1:ncon)
     con_dmp, var_dmp = dulmage_mendelsohn(igraph._graph, con_node_set)
@@ -413,13 +424,13 @@ function dulmage_mendelsohn(
         square = [nodes[n] for n in var_dmp[4]],
         overconstrained = [nodes[n] for n in var_dmp[3]],
     )
-    return con_dmp, var_dmp
+    return DulmageMendelsohnDecomposition((con_dmp, var_dmp))
 end
 
 function dulmage_mendelsohn(
     constraints::Vector{<:JuMP.ConstraintRef},
     variables::Vector{JuMP.VariableRef},
-)::Tuple{DMConPartition, DMVarPartition}
+)::DulmageMendelsohnDecomposition
     igraph = IncidenceGraphInterface(constraints, variables)
     return dulmage_mendelsohn(igraph)
 end
@@ -542,6 +553,15 @@ connected_components(model::JuMP.Model) = connected_components(IncidenceGraphInt
 connected_components(matrix::SparseMatrixCSC) = connected_components(IncidenceGraphInterface(matrix))
 connected_components(matrix::Matrix) = connected_components(IncidenceGraphInterface(matrix))
 
+const Subsystem = NamedTuple{
+    (:con, :var),
+    # NOTE: These fields are not fully typed as we need to support ints as
+    # wells as variables/constraints.
+    # TODO: Parameterize this type? Or just always use ints, then keep the vars/cons
+    # somewhere else?
+    Tuple{Vector,Vector},
+}
+
 """
     block_triangularize(igraph::IncidenceGraphInterface)::Vector{Tuple{Vector, Vector}}
 
@@ -631,7 +651,7 @@ function block_triangularize(
     matching = maximum_matching(igraph._graph, connodeset)
     nmatch = length(matching)
     if nvar != ncon || nvar != nmatch
-        @error (
+        error(
             "block_triangularize only supports square systems with perfect matchings."
             * "Got nvar=$nvar, ncon=$ncon, n. matched = $nmatch"
         )
@@ -639,7 +659,7 @@ function block_triangularize(
     connode_blocks = _block_triangularize(igraph._graph, matching)
     # node_blocks is a vector of tuples of vectors: [([cons], [vars]),...]
     blocks = [
-        ([igraph._nodes[n] for n in b], [igraph._nodes[matching[n]] for n in b])
+        Subsystem(([igraph._nodes[n] for n in b], [igraph._nodes[matching[n]] for n in b]))
         for b in connode_blocks
     ]
     return blocks
@@ -705,7 +725,7 @@ x[3], (x[1] * (x[3] ^ 0.5)) - 3.0 = 0
 function block_triangularize(
     constraints::Vector{<:JuMP.ConstraintRef},
     variables::Vector{JuMP.VariableRef},
-)::Vector{Tuple{Vector{JuMP.ConstraintRef}, Vector{JuMP.VariableRef}}}
+)::Vector{Subsystem}
     igraph = IncidenceGraphInterface(constraints, variables)
     return block_triangularize(igraph)
 end
diff --git a/src/visualize.jl b/src/visualize.jl
new file mode 100644
index 0000000..d0f2643
--- /dev/null
+++ b/src/visualize.jl
@@ -0,0 +1,89 @@
+#  ___________________________________________________________________________
+#
+#  MathProgIncidence.jl: Math Programming Incidence Graph Analysis
+#  Copyright (c) 2023. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001
+#  for Los Alamos National Laboratory (LANL), which is operated by Triad
+#  National Security, LLC for the U.S. Department of Energy/National Nuclear
+#  Security Administration. All rights in the program are reserved by Triad
+#  National Security, LLC, and the U.S. Department of Energy/National Nuclear
+#  Security Administration. The Government is granted for itself and others
+#  acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+#  in this material to reproduce, prepare derivative works, distribute copies
+#  to the public, perform publicly and display publicly, and to permit others
+#  to do so.
+#
+#  This software is distributed under the 3-clause BSD license.
+#  ___________________________________________________________________________
+
+"""
+Visualization of the decompositions implemented. These depend on the types
+defined in interface.jl
+"""
+
+function Base.show(io::IO, dm::DulmageMendelsohnDecomposition)
+    nvar = sum(map(length, dm.var))
+    ncon = sum(map(length, dm.con))
+    #uc_var = vcat(dm.var.underconstrained, dm.var.unmatched)
+    #uc_con = dm.con.underconstrained
+    #oc_var = dm.var.overconstrained
+    #oc_con = vcat(dm.con.overconstrained, dm.con.unmatched)
+    uc = Subsystem((dm.con.underconstrained, vcat(dm.var.unmatched, dm.var.underconstrained)))
+    oc = Subsystem((vcat(dm.con.overconstrained, dm.con.unmatched), dm.var.overconstrained))
+    wc = Subsystem((dm.con.square, dm.var.square))
+    msg = "Dulmage-Mendelsohn decomposition with $ncon constraints and $nvar variables"
+    nchar = length(msg)
+    println(io, msg)
+    println(io, repeat("-", nchar))
+    println(io, "Under-constrained subsystem: $(length(uc.con)) constraints x $(length(uc.var)) variables")
+    println(io, uc)
+    println(io, "Over-constrained subsystem:  $(length(oc.con)) constraints x $(length(oc.var)) variables")
+    println(io, oc)
+    println(io, "Well-constrained subsystem:  $(length(wc.con)) constraints x $(length(wc.var)) variables")
+    println(io, wc)
+    println(io, repeat("-", nchar))
+    return
+end
+
+function Base.show(io::IO, subsystem::Subsystem)
+    nvar = length(subsystem.var)
+    ncon = length(subsystem.con)
+    msg = "Subsystem with $nvar variables and $ncon constraints"
+    println(io, msg)
+    println(io, repeat("-", length(msg)))
+    println(io, "Variables:")
+    for var in subsystem.var
+        println(io, "  $var")
+    end
+    println(io, "Constraints:")
+    for con in subsystem.con
+        println(io, "  $con")
+    end
+    return
+end
+
+function Base.show(io::IO, subsystems::Vector{Subsystem})
+    n_subs = length(subsystems)
+    msg = "Set of $n_subs subsystems"
+    println(io, msg)
+    println(io, repeat("=", length(msg)))
+    for (i, sub) in enumerate(subsystems)
+        println()
+        nvar = length(sub.var)
+        ncon = length(sub.con)
+        msg = "Subsystem $i: $nvar variables, $ncon constraints"
+        println(io, msg)
+        println(io, repeat("-", length(msg)))
+        println(io, "Variables:")
+        for var in sub.var
+            println(io, "  $var")
+        end
+        println(io, "Constraints:")
+        for con in sub.con
+            println(io, "  $con")
+        end
+    end
+    println(io, repeat("=", length(msg)))
+    return
+end