diff --git a/R/map.R b/R/map.R
index 77e2f8ec..e254d451 100644
--- a/R/map.R
+++ b/R/map.R
@@ -218,6 +218,41 @@ map_raw <- function(.x, .f, ...) {
   .Call(map_impl, environment(), ".x", ".f", "raw")
 }
 
+#' @rdname map
+#' @export
+map_lgl_matrix <- function(.x, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map_matrix_impl, environment(), ".x", ".f", "logical", .n, .by_row)
+}
+
+#' @rdname map
+#' @export
+map_chr_matrix <- function(.x, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map_matrix_impl, environment(), ".x", ".f", "character", .n, .by_row)
+}
+
+#' @rdname map
+#' @export
+map_int_matrix <- function(.x, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map_matrix_impl, environment(), ".x", ".f", "integer", .n, .by_row)
+}
+
+#' @rdname map
+#' @export
+map_dbl_matrix <- function(.x, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map_matrix_impl, environment(), ".x", ".f", "double", .n, .by_row)
+}
+
+#' @rdname map
+#' @export
+map_raw_matrix <- function(.x, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map_matrix_impl, environment(), ".x", ".f", "raw", .n, .by_row)
+}
+
 #' @rdname map
 #' @param .id Either a string or `NULL`. If a string, the output will contain
 #'   a variable with that name, storing either the name (if `.x` is named) or
diff --git a/R/map2-pmap.R b/R/map2-pmap.R
index 76030e25..e850ac27 100644
--- a/R/map2-pmap.R
+++ b/R/map2-pmap.R
@@ -133,6 +133,38 @@ map2_raw <- function(.x, .y, .f, ...) {
   .f <- as_mapper(.f, ...)
   .Call(map2_impl, environment(), ".x", ".y", ".f", "raw")
 }
+
+#' @export
+#' @rdname map2
+map2_lgl_matrix <- function(.x, .y, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map2_matrix_impl, environment(), ".x", ".y", ".f", "logical", .n, .by_row)
+}
+#' @export
+#' @rdname map2
+map2_int_matrix <- function(.x, .y, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map2_matrix_impl, environment(), ".x", ".y", ".f", "integer", .n, .by_row)
+}
+#' @export
+#' @rdname map2
+map2_dbl_matrix <- function(.x, .y, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map2_matrix_impl, environment(), ".x", ".y", ".f", "double", .n, .by_row)
+}
+#' @export
+#' @rdname map2
+map2_chr_matrix <- function(.x, .y, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map2_matrix_impl, environment(), ".x", ".y", ".f", "character", .n, .by_row)
+}
+#' @export
+#' @rdname map2
+map2_raw_matrix <- function(.x, .y, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  .Call(map2_matrix_impl, environment(), ".x", ".y", ".f", "raw", .n, .by_row)
+}
+
 #' @rdname map2
 #' @export
 map2_dfr <- function(.x, .y, .f, ..., .id = NULL) {
@@ -228,6 +260,61 @@ pmap_raw <- function(.l, .f, ...) {
   .Call(pmap_impl, environment(), ".l", ".f", "raw")
 }
 
+#' @export
+#' @rdname map2
+pmap_lgl_matrix <- function(.l, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  if (is.data.frame(.l)) {
+    .l <- as.list(.l)
+  }
+
+  .Call(pmap_matrix_impl, environment(), ".l", ".f", "logical", .n, .by_row)
+}
+
+#' @export
+#' @rdname map2
+pmap_int_matrix <- function(.l, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  if (is.data.frame(.l)) {
+    .l <- as.list(.l)
+  }
+
+  .Call(pmap_matrix_impl, environment(), ".l", ".f", "integer", .n, .by_row)
+}
+
+#' @export
+#' @rdname map2
+pmap_dbl_matrix <- function(.l, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  if (is.data.frame(.l)) {
+    .l <- as.list(.l)
+  }
+
+  .Call(pmap_matrix_impl, environment(), ".l", ".f", "double", .n, .by_row)
+}
+
+#' @export
+#' @rdname map2
+pmap_chr_matrix <- function(.l, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  if (is.data.frame(.l)) {
+    .l <- as.list(.l)
+  }
+
+  .Call(pmap_matrix_impl, environment(), ".l", ".f", "character", .n, .by_row)
+}
+
+#' @export
+#' @rdname map2
+pmap_raw_matrix <- function(.l, .f, .n, ..., .by_row = FALSE) {
+  .f <- as_mapper(.f, ...)
+  if (is.data.frame(.l)) {
+    .l <- as.list(.l)
+  }
+
+  .Call(pmap_matrix_impl, environment(), ".l", ".f", "raw", .n, .by_row)
+}
+
 #' @rdname map2
 #' @export
 pmap_dfr <- function(.l, .f, ..., .id = NULL) {
diff --git a/src/coerce.c b/src/coerce.c
index 86ea0055..66f6a28b 100644
--- a/src/coerce.c
+++ b/src/coerce.c
@@ -103,7 +103,6 @@ void set_vector_value(SEXP to, int i, SEXP from, int j) {
   }
 }
 
-
 SEXP coerce_impl(SEXP x, SEXP type_) {
   int n = Rf_length(x);
 
diff --git a/src/init.c b/src/init.c
index 7c08b805..871b777d 100644
--- a/src/init.c
+++ b/src/init.c
@@ -15,8 +15,11 @@ extern SEXP coerce_impl(SEXP, SEXP);
 extern SEXP pluck_impl(SEXP, SEXP, SEXP, SEXP);
 extern SEXP flatten_impl(SEXP);
 extern SEXP map_impl(SEXP, SEXP, SEXP, SEXP);
+extern SEXP map_matrix_impl(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP map2_impl(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP map2_matrix_impl(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP pmap_impl(SEXP, SEXP, SEXP, SEXP);
+extern SEXP pmap_matrix_impl(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP transpose_impl(SEXP, SEXP);
 extern SEXP vflatten_impl(SEXP, SEXP);
 
@@ -26,8 +29,11 @@ static const R_CallMethodDef CallEntries[] = {
   {"pluck_impl",            (DL_FUNC) &pluck_impl,   4},
   {"flatten_impl",          (DL_FUNC) &flatten_impl,   1},
   {"map_impl",              (DL_FUNC) &map_impl,       4},
+  {"map_matrix_impl",       (DL_FUNC) &map_matrix_impl, 6},
   {"map2_impl",             (DL_FUNC) &map2_impl,      5},
+  {"map2_matrix_impl",      (DL_FUNC) &map2_matrix_impl, 7},
   {"pmap_impl",             (DL_FUNC) &pmap_impl,      4},
+  {"pmap_matrix_impl",      (DL_FUNC) &pmap_matrix_impl, 6},
   {"transpose_impl",        (DL_FUNC) &transpose_impl, 2},
   {"vflatten_impl",         (DL_FUNC) &vflatten_impl,  2},
   {"purrr_eval",            (DL_FUNC) &Rf_eval,        2},
diff --git a/src/map.c b/src/map.c
index e14b2af0..65cfaea6 100644
--- a/src/map.c
+++ b/src/map.c
@@ -57,6 +57,49 @@ SEXP call_loop(SEXP env, SEXP call, int n, SEXPTYPE type, int force_args) {
   return out;
 }
 
+// call must involve i
+SEXP call_loop_matrix(SEXP env, SEXP call, int n, int m, SEXPTYPE type, int force_args, int by_row) {
+  // Create variable "i" and map to scalar integer
+  SEXP i_val = PROTECT(Rf_ScalarInteger(1));
+  SEXP i = Rf_install("i");
+  Rf_defineVar(i, i_val, env);
+
+  SEXP out;
+  if (by_row) {
+    out = PROTECT(Rf_allocMatrix(type, n, m));
+  } else {
+    out = PROTECT(Rf_allocMatrix(type, m, n));
+  }
+  for (int i = 0; i < n; ++i) {
+    if (i % 1024 == 0)
+      R_CheckUserInterrupt();
+
+    INTEGER(i_val)[0] = i + 1;
+
+#if defined(R_VERSION) && R_VERSION >= R_Version(3, 2, 3)
+    SEXP res = PROTECT(R_forceAndCall(call, force_args, env));
+#else
+    SEXP res = PROTECT(Rf_eval(call, env));
+#endif
+    if (type != VECSXP && Rf_length(res) != m) {
+      SEXP ptype = PROTECT(Rf_allocVector(type, 0));
+      stop_bad_element_vector(res, i + 1, ptype, m, "Result", NULL, false);
+    }
+
+    for (int j = 0; j < m; j++) {
+      if (by_row) {
+        set_vector_value(out, n * j + i, res, j);
+      } else {
+        set_vector_value(out, m * i + j, res, j);
+      }
+    }
+    UNPROTECT(1);
+  }
+
+  UNPROTECT(2);
+  return out;
+}
+
 SEXP map_impl(SEXP env, SEXP x_name_, SEXP f_name_, SEXP type_) {
   const char* x_name = CHAR(Rf_asChar(x_name_));
   const char* f_name = CHAR(Rf_asChar(f_name_));
@@ -91,6 +134,42 @@ SEXP map_impl(SEXP env, SEXP x_name_, SEXP f_name_, SEXP type_) {
   return out;
 }
 
+SEXP map_matrix_impl(SEXP env, SEXP x_name_, SEXP f_name_, SEXP type_, SEXP n_elem_, SEXP by_row_) {
+  const char* x_name = CHAR(Rf_asChar(x_name_));
+  const char* f_name = CHAR(Rf_asChar(f_name_));
+  const int n_elem = Rf_asInteger(n_elem_);
+  const int by_row = Rf_asLogical(by_row_);
+
+  SEXP x = Rf_install(x_name);
+  SEXP f = Rf_install(f_name);
+  SEXP i = Rf_install("i");
+  SEXPTYPE type = Rf_str2type(CHAR(Rf_asChar(type_)));
+
+  SEXP x_val = PROTECT(Rf_eval(x, env));
+  check_vector(x_val, ".x");
+
+  int n = Rf_length(x_val);
+  if (n == 0) {
+    SEXP out = PROTECT(Rf_allocVector(type, 0));
+    copy_names(x_val, out);
+    UNPROTECT(2);
+    return out;
+  }
+
+  // Constructs a call like f(x[[i]], ...) - don't want to substitute
+  // actual values for f or x, because they may be long, which creates
+  // bad tracebacks()
+  SEXP Xi = PROTECT(Rf_lang3(R_Bracket2Symbol, x, i));
+  SEXP f_call = PROTECT(Rf_lang3(f, Xi, R_DotsSymbol));
+
+  SEXP out = PROTECT(call_loop_matrix(env, f_call, n, n_elem, type, 1, by_row));
+  copy_names(x_val, out);
+
+  UNPROTECT(4);
+
+  return out;
+}
+
 SEXP map2_impl(SEXP env, SEXP x_name_, SEXP y_name_, SEXP f_name_, SEXP type_) {
   const char* x_name = CHAR(Rf_asChar(x_name_));
   const char* y_name = CHAR(Rf_asChar(y_name_));
@@ -137,6 +216,54 @@ SEXP map2_impl(SEXP env, SEXP x_name_, SEXP y_name_, SEXP f_name_, SEXP type_) {
   return out;
 }
 
+SEXP map2_matrix_impl(SEXP env, SEXP x_name_, SEXP y_name_, SEXP f_name_, SEXP type_, SEXP n_elem_, SEXP by_row_) {
+  const char* x_name = CHAR(Rf_asChar(x_name_));
+  const char* y_name = CHAR(Rf_asChar(y_name_));
+  const char* f_name = CHAR(Rf_asChar(f_name_));
+  const int n_elem = Rf_asInteger(n_elem_);
+  const int by_row = Rf_asLogical(by_row_);
+
+  SEXP x = Rf_install(x_name);
+  SEXP y = Rf_install(y_name);
+  SEXP f = Rf_install(f_name);
+  SEXP i = Rf_install("i");
+  SEXPTYPE type = Rf_str2type(CHAR(Rf_asChar(type_)));
+
+  SEXP x_val = PROTECT(Rf_eval(x, env));
+  check_vector(x_val, ".x");
+  SEXP y_val = PROTECT(Rf_eval(y, env));
+  check_vector(y_val, ".y");
+
+  int nx = Rf_length(x_val), ny = Rf_length(y_val);
+  if (nx == 0 || ny == 0) {
+    SEXP out = PROTECT(Rf_allocVector(type, 0));
+    copy_names(x_val, out);
+    UNPROTECT(3);
+    return out;
+  }
+  if (nx != ny && !(nx == 1 || ny == 1)) {
+    Rf_errorcall(R_NilValue,
+                 "Mapped vectors must have consistent lengths:\n"
+                 "* `.x` has length %d\n"
+                 "* `.y` has length %d",
+                 nx,
+                 ny);
+  }
+  int n = (nx > ny) ? nx : ny;
+
+  // Constructs a call like f(x[[i]], y[[i]], ...)
+  SEXP one = PROTECT(Rf_ScalarInteger(1));
+  SEXP Xi = PROTECT(Rf_lang3(R_Bracket2Symbol, x, nx == 1 ? one : i));
+  SEXP Yi = PROTECT(Rf_lang3(R_Bracket2Symbol, y, ny == 1 ? one : i));
+  SEXP f_call = PROTECT(Rf_lang4(f, Xi, Yi, R_DotsSymbol));
+
+  SEXP out = PROTECT(call_loop_matrix(env, f_call, n, n_elem, type, 2, by_row));
+  copy_names(x_val, out);
+
+  UNPROTECT(7);
+  return out;
+}
+
 SEXP pmap_impl(SEXP env, SEXP l_name_, SEXP f_name_, SEXP type_) {
   const char* l_name = CHAR(Rf_asChar(l_name_));
   SEXP l = Rf_install(l_name);
@@ -229,3 +356,98 @@ SEXP pmap_impl(SEXP env, SEXP l_name_, SEXP f_name_, SEXP type_) {
   UNPROTECT(5);
   return out;
 }
+
+SEXP pmap_matrix_impl(SEXP env, SEXP l_name_, SEXP f_name_, SEXP type_, SEXP n_elem_, SEXP by_row_) {
+  const char* l_name = CHAR(Rf_asChar(l_name_));
+  SEXP l = Rf_install(l_name);
+  SEXP l_val = PROTECT(Rf_eval(l, env));
+  SEXPTYPE type = Rf_str2type(CHAR(Rf_asChar(type_)));
+  const int n_elem = Rf_asInteger(n_elem_);
+  const int by_row = Rf_asLogical(by_row_);
+
+  if (!Rf_isVectorList(l_val)) {
+    stop_bad_type(l_val, "a list", NULL, l_name);
+  }
+
+  // Check all elements are lists and find maximum length
+  int m = Rf_length(l_val);
+  int n = 0;
+  for (int j = 0; j < m; ++j) {
+    SEXP j_val = VECTOR_ELT(l_val, j);
+
+    if (!Rf_isVector(j_val) && !Rf_isNull(j_val)) {
+      stop_bad_element_type(j_val, j + 1, "a vector", NULL, l_name);
+    }
+
+    int nj = Rf_length(j_val);
+
+    if (nj == 0) {
+      SEXP out = PROTECT(Rf_allocVector(type, 0));
+      copy_names(j_val, out);
+      UNPROTECT(2);
+      return out;
+    }
+
+    if (nj > n) {
+      n = nj;
+    }
+
+  }
+
+  // Check length of all elements
+  for (int j = 0; j < m; ++j) {
+    SEXP j_val = VECTOR_ELT(l_val, j);
+    int nj = Rf_length(j_val);
+
+    if (nj != 1 && nj != n) {
+      stop_bad_element_length(j_val, j + 1, n, NULL, ".l", true);
+    }
+  }
+
+  SEXP l_names = PROTECT(Rf_getAttrib(l_val, R_NamesSymbol));
+  int has_names = !Rf_isNull(l_names);
+
+  const char* f_name = CHAR(Rf_asChar(f_name_));
+  SEXP f = Rf_install(f_name);
+  SEXP i = Rf_install("i");
+  SEXP one = PROTECT(Rf_ScalarInteger(1));
+
+  // Construct call like f(.l[[1]][[i]], .l[[2]][[i]], ...)
+  //
+  // Currently accessing S3 vectors in a list like .l[[c(1, i)]] will not
+  // preserve the class (cf. #358).
+  //
+  // We construct the call backwards because can only add to the front of a
+  // linked list. That makes PROTECTion tricky because we need to update it
+  // each time to point to the start of the linked list.
+
+  SEXP f_call = Rf_lang1(R_DotsSymbol);
+  PROTECT_INDEX fi;
+  PROTECT_WITH_INDEX(f_call, &fi);
+
+  for (int j = m - 1; j >= 0; --j) {
+    int nj = Rf_length(VECTOR_ELT(l_val, j));
+
+    // Construct call like .l[[j]][[i]]
+    SEXP j_ = PROTECT(Rf_ScalarInteger(j + 1));
+    SEXP l_j = PROTECT(Rf_lang3(R_Bracket2Symbol, l, j_));
+    SEXP l_ji = PROTECT(Rf_lang3(R_Bracket2Symbol, l_j, nj == 1 ? one : i));
+
+    REPROTECT(f_call = Rf_lcons(l_ji, f_call), fi);
+    if (has_names && CHAR(STRING_ELT(l_names, j))[0] != '\0')
+      SET_TAG(f_call, Rf_install(CHAR(STRING_ELT(l_names, j))));
+
+    UNPROTECT(3);
+  }
+
+  REPROTECT(f_call = Rf_lcons(f, f_call), fi);
+
+  SEXP out = PROTECT(call_loop_matrix(env, f_call, n, n_elem, type, m, by_row));
+
+  if (Rf_length(l_val)) {
+    copy_names(VECTOR_ELT(l_val, 0), out);
+  }
+
+  UNPROTECT(5);
+  return out;
+}