diff --git a/DESCRIPTION b/DESCRIPTION
index 6dc3309b..273806ce 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,5 +1,5 @@
 Package: optmatch
-Version: 0.10.8.9001
+Version: 0.10.8.9002
 Title: Functions for Optimal Matching
 Description: Distance based bipartite matching using minimum cost flow, oriented
     to matching of treatment and control groups in observational studies ('Hansen'
diff --git a/R/InfinitySparseMatrix.R b/R/InfinitySparseMatrix.R
index 91d831fe..c30825cd 100644
--- a/R/InfinitySparseMatrix.R
+++ b/R/InfinitySparseMatrix.R
@@ -1103,3 +1103,34 @@ as.list.DenseMatrix <- function(x, ...) {
            })
   return(x)
 }
+
+##' This matches the syntax and semantics of
+##' subset for matrices.
+##'
+##' @title Subsetting for BlockedInfinitySparseMatrices
+##' @param x BlockedInfinitySparseMatrix to be subset or bound.
+##' @param subset Logical expression indicating rows to keep.
+##' @param select Logical expression indicating columns to keep.
+##' @param ... Other arguments are ignored.
+##' @return If groups has names, a BlockedInfinitySparseMatrix with only
+##'         the selected elements, otherwise an InfinitySparsematrix with
+##'         only the selected elements
+##' @rdname bism.subset
+##' @export
+subset.BlockedInfinitySparseMatrix <- function(x, subset, select, ...) {
+    subIsm <- callGeneric(as(x, "InfinitySparseMatrix"), subset, select)
+    oldNames <- names(x@groups)
+    if (!is.null(oldNames)) {   # we can use the groups names to subset groups
+        subNames <- oldNames[which((oldNames %in% subIsm@rownames) |
+                                   (oldNames %in% subIsm@colnames))]
+        subGroups <- x@groups[subNames]
+        subObj <- new("BlockedInfinitySparseMatrix",
+                      subIsm, groups = subGroups)
+    } else {
+        # since groups doesn't have names, we can't meaningfully subset it
+        # groups is meaningless for the subsetted matrix
+        # demote object to ISM
+        subObj <- subIsm
+    }
+    return(subObj)
+}
diff --git a/R/summary.ism.R b/R/summary.ism.R
index 03836c58..47903ddc 100644
--- a/R/summary.ism.R
+++ b/R/summary.ism.R
@@ -102,7 +102,7 @@ summary.BlockedInfinitySparseMatrix <- function(object, ...,
   out <- lapply(levels(object@groups),
                 function(x) {
                   thisgroup <- names(object@groups[object@groups == x])
-                  ism <- subset(object,
+                  ism <- subset.InfinitySparseMatrix(object,
                                 subset=object@rownames %in% thisgroup,
                                 select=object@colnames %in% thisgroup)
                   s <- summary(ism, ..., distanceSummary=distanceSummary)
diff --git a/tests/testthat/test.InfinitySparseMatrix.R b/tests/testthat/test.InfinitySparseMatrix.R
index 52095bee..6ec90c8b 100644
--- a/tests/testthat/test.InfinitySparseMatrix.R
+++ b/tests/testthat/test.InfinitySparseMatrix.R
@@ -264,7 +264,7 @@ test_that("#190: agreement in dimension names", {
 
 })
 
-test_that("Subsetting", {
+test_that("ISM Subsetting", {
   m <- matrix(c(1,Inf, 2, 3), nrow = 2, ncol = 2)
   rownames(m) <- c("A", "B")
   colnames(m) <- c("C", "D")
@@ -400,6 +400,23 @@ test_that("BlockedISM addition", {
     "BlockedInfinitySparseMatrix")
 })
 
+test_that("BlockedISM Subsetting", {
+  Z <- rep(c(0,1), 8)
+  B <- rep(1:4, each = 4)
+
+  res.b <- exactMatch(Z ~ B)
+  sub.b <- subset(res.b,
+                   c(TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE),
+                   c(TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE))
+
+  expect_is(sub.b, "BlockedInfinitySparseMatrix")
+  expect_false(is.null(sub.b@groups))
+  expect_equal(names(sub.b@groups),
+               c("1", "2", "5", "6", "9", "10", "13", "14"))
+  expect_equal(rownames(sub.b), c("2", "6", "10", "14"))
+  expect_equal(colnames(sub.b), c("1", "5", "9", "13"))
+})
+
 test_that("Get subproblem size of each block", {
   Z <- rep(c(0,1), 8)
   B1 <- c(rep('a',3),rep('b', 3), rep('c', 6), rep('d', 4))