diff --git a/.gitignore b/.gitignore
index c586de43..a6f31d68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,5 @@ renv.lock
 tmaps
 src/*.so
 src/*.o
-test.R
\ No newline at end of file
+test.R
+test.cpp
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index d15dd0d5..789ab7e6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -118,7 +118,7 @@ Suggests:
     umap,
     UCell,
     withr
-LinkingTo: Rcpp
+LinkingTo: Rcpp,RcppArmadillo
 RoxygenNote: 7.2.3
 Config/testthat/edition: 3
 URL: https://github.com/zhanghao-njmu/SCP
diff --git a/NAMESPACE b/NAMESPACE
index 563d2be0..4b0ec076 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -70,6 +70,7 @@ export(LineagePlot)
 export(ListDB)
 export(MNN_integrate)
 export(MatrixConvert)
+export(MergeRows)
 export(PAGAPlot)
 export(PrepareDB)
 export(PrepareEnv)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 258e9cc8..cd49f94c 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,7 +1,19 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-asMatrix <- function(rp, cp, z, nrows, ncols) {
-    .Call('_SCP_asMatrix', PACKAGE = 'SCP', rp, cp, z, nrows, ncols)
+MergeRowsForDense <- function(matrix, groupings) {
+    .Call('_SCP_MergeRowsForDense', PACKAGE = 'SCP', matrix, groupings)
+}
+
+MergeRowsForSparse <- function(matrix, groupings) {
+    .Call('_SCP_MergeRowsForSparse', PACKAGE = 'SCP', matrix, groupings)
+}
+
+MergeRowsC <- function(matrix, groupings) {
+    .Call('_SCP_MergeRowsC', PACKAGE = 'SCP', matrix, groupings)
+}
+
+asMatrixC <- function(rp, cp, z, nrows, ncols) {
+    .Call('_SCP_asMatrixC', PACKAGE = 'SCP', rp, cp, z, nrows, ncols)
 }
 
diff --git a/R/SCP-analysis.R b/R/SCP-analysis.R
index 1267c439..7150d39c 100644
--- a/R/SCP-analysis.R
+++ b/R/SCP-analysis.R
@@ -476,6 +476,7 @@ MatrixConvert <- function(matrix, aggregate_method = "sum", geneID_from_IDtype =
   return(matrix)
 }
 
+
 #' @export
 FindCorrelatedFeatures <- function(srt, feature, method = "ejaccard", assay = NULL, slot = "data") {
   assay <- assay %||% DefaultAssay(srt)
diff --git a/R/SCP-app.R b/R/SCP-app.R
index 94900483..e8f165e8 100644
--- a/R/SCP-app.R
+++ b/R/SCP-app.R
@@ -248,7 +248,7 @@ CreateMetaFile <- function(srt, MetaFile, name = NULL, write_tools = FALSE, writ
 #' @importFrom Seurat Reductions Assays DefaultAssay
 #' @export
 PrepareSCExplorer <- function(object,
-                              base_dir = "SCExplorer", DataFile = "Data.hdf5", MetaFile = "Meta.hdf5",
+                              base_dir = "SCExplorer",
                               assays = "RNA", slots = c("counts", "data"),
                               ignore_nlevel = 100, write_tools = FALSE, write_misc = FALSE,
                               compression_level = 6, overwrite = FALSE) {
@@ -257,8 +257,8 @@ PrepareSCExplorer <- function(object,
     message("Create SCExplorer base directory: ", base_dir)
     dir.create(base_dir, recursive = TRUE, showWarnings = FALSE)
   }
-  DataFile_full <- paste0(base_dir, "/", DataFile)
-  MetaFile_full <- paste0(base_dir, "/", MetaFile)
+  DataFile_full <- paste0(base_dir, "/Data.hdf5")
+  MetaFile_full <- paste0(base_dir, "/Meta.hdf5")
 
   if (!is.list(object)) {
     object <- list(object)
@@ -529,7 +529,6 @@ CreateSeuratObject2 <- function(counts, project = "SeuratProject", assay = "RNA"
 #' @param create_script A logical. Whether to create the SCExplorer app script. Default is TRUE.
 #' @param style_script A logical. Whether to style the SCExplorer app script. Default is TRUE.
 #' @param overwrite A logical. Whether to overwrite existing files. Default is FALSE.
-#' @param return_app A logical. Whether to return the SCExplorer app. Default is TRUE.
 #'
 #' @seealso \code{\link{CreateDataFile}} \code{\link{CreateMetaFile}} \code{\link{PrepareSCExplorer}} \code{\link{FetchH5}}
 #'
@@ -579,8 +578,6 @@ CreateSeuratObject2 <- function(counts, project = "SeuratProject", assay = "RNA"
 #' @importFrom utils packageVersion
 #' @export
 RunSCExplorer <- function(base_dir = "SCExplorer",
-                          DataFile = "Data.hdf5",
-                          MetaFile = "Meta.hdf5",
                           title = "SCExplorer",
                           initial_dataset = NULL,
                           initial_reduction = NULL,
@@ -598,13 +595,14 @@ RunSCExplorer <- function(base_dir = "SCExplorer",
                           initial_raster = NULL,
                           session_workers = 2,
                           plotting_workers = 8,
-                          create_script = TRUE,
+                          create_script = FALSE,
                           style_script = require("styler", quietly = TRUE),
-                          overwrite = FALSE,
-                          return_app = TRUE) {
+                          overwrite = FALSE) {
   check_R(c("rhdf5", "HDF5Array", "shiny@1.6.0", "ggplot2", "ragg", "htmlwidgets", "plotly", "bslib", "future", "promises", "BiocParallel"))
-  DataFile_full <- paste0(base_dir, "/", DataFile)
-  MetaFile_full <- paste0(base_dir, "/", MetaFile)
+  base_dir <- normalizePath(base_dir, mustWork = FALSE)
+  DataFile_full <- paste0(base_dir, "/Data.hdf5")
+  MetaFile_full <- paste0(base_dir, "/Meta.hdf5")
+
   if (!file.exists(DataFile_full) || !file.exists(MetaFile_full)) {
     stop("Please create the DataFile and MetaFile using PrepareSCExplorer function first!")
   }
@@ -728,7 +726,7 @@ ui <- fluidPage(
           ),
           fluidRow(
             column(
-              width = 6, align = "center",
+              width = 4, align = "center",
               radioButtons(
                 inputId = "label1",
                 label = "Label",
@@ -738,7 +736,7 @@ ui <- fluidPage(
               )
             ),
             column(
-              width = 6, align = "center",
+              width = 4, align = "center",
               radioButtons(
                 inputId = "raster1",
                 label = "Raster",
@@ -746,6 +744,16 @@ ui <- fluidPage(
                 selected = initial_raster,
                 inline = TRUE
               )
+            ),
+            column(
+              width = 4, align = "center",
+              radioButtons(
+                inputId = "dimension1",
+                label = "Dimension",
+                choices = c("2D", "3D"),
+                selected = "2D",
+                inline = TRUE
+              )
             )
           ),
           fluidRow(
@@ -757,7 +765,7 @@ ui <- fluidPage(
                 value = 1,
                 min = 0.1,
                 max = 10,
-                step = 0.5,
+                step = 0.1,
                 width = "150px"
               )
             ),
@@ -811,27 +819,13 @@ ui <- fluidPage(
         mainPanel(
           width = 9,
           fluidPage(
-            tabsetPanel(
-              tabPanel(
-                title = "2D plot",
-                column(
-                  width = 12, offset = 0, style = "padding:0px;margin:0%",
-                  div(
-                    style = "overflow-x: auto;",
-                    uiOutput("plot1")
-                  )
-                )
-              ),
-              tabPanel(
-                title = "3D plot",
-                column(
-                  width = 12, offset = 0, style = "padding:0px;margin:0%",
-                  div(
-                    style = "overflow-x: auto;",
-                    plotly::plotlyOutput("plot1_3d", height = "100%", width = "100%")
-                  )
-                )
-              )
+            conditionalPanel(
+              condition = "input.dimension1 == \'2D\'",
+              uiOutput("plot1")
+            ),
+            conditionalPanel(
+              condition = "input.dimension1 == \'3D\'",
+              plotly::plotlyOutput("plot1_3d", height = "100%", width = "100%")
             )
           )
         )
@@ -909,7 +903,7 @@ ui <- fluidPage(
           ),
           fluidRow(
             column(
-              width = 4, align = "center",
+              width = 6, align = "center",
               radioButtons(
                 inputId = "coExp2",
                 label = "Co-expression",
@@ -919,7 +913,7 @@ ui <- fluidPage(
               ),
             ),
             column(
-              width = 4, align = "center",
+              width = 6, align = "center",
               radioButtons(
                 inputId = "scale2",
                 label = "Color scale",
@@ -927,9 +921,11 @@ ui <- fluidPage(
                 selected = "feature",
                 inline = TRUE
               ),
-            ),
+            )
+          ),
+          fluidRow(
             column(
-              width = 4, align = "center",
+              width = 6, align = "center",
               radioButtons(
                 inputId = "raster2",
                 label = "Raster",
@@ -937,6 +933,16 @@ ui <- fluidPage(
                 selected = initial_raster,
                 inline = TRUE
               )
+            ),
+            column(
+              width = 6, align = "center",
+              radioButtons(
+                inputId = "dimension2",
+                label = "Dimension",
+                choices = c("2D", "3D"),
+                selected = "2D",
+                inline = TRUE
+              )
             )
           ),
           fluidRow(
@@ -948,7 +954,7 @@ ui <- fluidPage(
                 value = 1,
                 min = 0.1,
                 max = 10,
-                step = 0.5,
+                step = 0.1,
                 width = "150px"
               )
             ),
@@ -1001,26 +1007,14 @@ ui <- fluidPage(
         ),
         mainPanel(
           width = 9,
-          tabsetPanel(
-            tabPanel(
-              title = "2D plot",
-              column(
-                width = 12, offset = 0, style = "padding:0px;margin:0%",
-                div(
-                  style = "overflow-x: auto;",
-                  uiOutput("plot2")
-                )
-              )
+          fluidPage(
+            conditionalPanel(
+              condition = "input.dimension2 == \'2D\'",
+              uiOutput("plot2")
             ),
-            tabPanel(
-              title = "3D plot",
-              column(
-                width = 12, offset = 0, style = "padding:0px;margin:0%",
-                div(
-                  style = "overflow-x: auto;",
-                  plotly::plotlyOutput("plot2_3d", height = "100%", width = "100%")
-                )
-              )
+            conditionalPanel(
+              condition = "input.dimension2 == \'3D\'",
+              plotly::plotlyOutput("plot2_3d", height = "100%", width = "100%")
             )
           )
         )
@@ -1214,20 +1208,7 @@ ui <- fluidPage(
         ),
         mainPanel(
           width = 9,
-          fluidPage(
-            tabsetPanel(
-              tabPanel(
-                title = "Statistical plot",
-                column(
-                  width = 12, offset = 0, style = "padding:0px;margin:0%",
-                  div(
-                    style = "overflow-x: auto;",
-                    uiOutput("plot3")
-                  )
-                )
-              )
-            )
-          )
+          fluidPage(uiOutput("plot3"))
         )
       )
     ),
@@ -1486,20 +1467,7 @@ ui <- fluidPage(
         ),
         mainPanel(
           width = 9,
-          fluidPage(
-            tabsetPanel(
-              tabPanel(
-                title = "Statistical plot",
-                column(
-                  width = 12, offset = 0, style = "padding:0px;margin:0%",
-                  div(
-                    style = "overflow-x: auto;",
-                    uiOutput("plot4")
-                  )
-                )
-              )
-            )
-          )
+          fluidPage(uiOutput("plot4"))
         )
       )
     )
@@ -1633,8 +1601,9 @@ server <- function(input, output, session) {
     } else {
       split1 <- input$split1
     }
-    label1 <- input$label1
-    raster1 <- input$raster1
+    label1 <- as.logical(input$label1)
+    raster1 <- as.logical(input$raster1)
+    dimension1 <- input$dimension1
     palette1 <- input$palette1
     theme1 <- input$theme1
     size1 <- input$size1
@@ -1659,24 +1628,26 @@ server <- function(input, output, session) {
 
         theme1 <- get(theme1, envir = asNamespace(themes[theme1]))
 
-        # print(">>> plot:")
-        # print(system.time(
-        p1_dim <- SCP::CellDimPlot(srt_tmp,
-          group.by = group1, split.by = split1, reduction = reduction1, raster = raster1, pt.size = pt_size1,
-          label = label1, palette = palette1, theme_use = theme1,
-          ncol = ncol1, byrow = byrow1, force = TRUE
-        )
-        # ))
+        plot3d <- max(sapply(names(srt_tmp@reductions), function(r) dim(srt_tmp[[r]])[2])) >= 3 & dimension1 == "3D"
 
-        # print(">>> panel_fix:")
-        # print(system.time(
-        p1_dim <- SCP::panel_fix(SCP::slim_data(p1_dim), height = size1, units = "in", raster = panel_raster, BPPARAM = BPPARAM, verbose = FALSE)
-        # ))
-        attr(p1_dim, "dpi") <- 300
-        plot3d <- max(sapply(names(srt_tmp@reductions), function(r) dim(srt_tmp[[r]])[2])) >= 3
         if (isTRUE(plot3d)) {
+          p1_dim <- NULL
           p1_3d <- SCP::CellDimPlot3D(srt_tmp, group.by = group1, pt.size = pt_size1 * 2, reduction = reduction1, palette = palette1, force = TRUE)
         } else {
+          # print(">>> plot:")
+          # print(system.time(
+          p1_dim <- SCP::CellDimPlot(srt_tmp,
+            group.by = group1, split.by = split1, reduction = reduction1, raster = raster1, pt.size = pt_size1,
+            label = label1, palette = palette1, theme_use = theme1,
+            ncol = ncol1, byrow = byrow1, force = TRUE
+          )
+          # ))
+
+          # print(">>> panel_fix:")
+          # print(system.time(
+          p1_dim <- SCP::panel_fix(SCP::slim_data(p1_dim), height = size1, units = "in", raster = panel_raster, BPPARAM = BPPARAM, verbose = FALSE)
+          # ))
+          attr(p1_dim, "dpi") <- 300
           p1_3d <- NULL
         }
         return(list(p1_dim, p1_3d))
@@ -1686,7 +1657,7 @@ server <- function(input, output, session) {
   }) %>%
     bindCache(
       input$dataset1, input$reduction1, input$group1, input$split1,
-      input$palette1, input$theme1, input$label1, input$raster1,
+      input$palette1, input$theme1, input$label1, input$raster1, input$dimension1,
       input$pt_size1, input$size1, input$ncol1, input$arrange1
     ) %>%
     bindEvent(input$submit1, ignoreNULL = FALSE, ignoreInit = FALSE)
@@ -1697,31 +1668,36 @@ server <- function(input, output, session) {
     r1()$then(function(x) {
       promisedData[["p1_dim"]] <- x[[1]]
       promisedData[["p1_3d"]] <- x[[2]]
-      width <- get_attr(x[[1]], "width")
-      height <- get_attr(x[[1]], "height")
-      dpi <- get_attr(x[[1]], "dpi")
-
       prog$set(value = 8, message = "Render plot...", detail = "[Cell dimensional reduction plot]")
-      # print("renderPlot:")
-      # print(system.time(
-      output$plot1 <- renderUI({
-        renderPlot(
-          {
-            x[[1]]
-          },
-          width = width * 96,
-          height = height * 96,
-          res = 96
-        )
-      })
-      # ))
 
-      # print("renderPlotly:")
-      # print(system.time(
-      output$plot1_3d <- plotly::renderPlotly({
-        x[[2]]
-      })
-      # ))
+      if (!is.null(x[[1]])) {
+        width <- get_attr(x[[1]], "width")
+        height <- get_attr(x[[1]], "height")
+        dpi <- get_attr(x[[1]], "dpi")
+
+        # print("renderPlot:")
+        # print(system.time(
+        output$plot1 <- renderUI({
+          renderPlot(
+            {
+              x[[1]]
+            },
+            width = width * 96,
+            height = height * 96,
+            res = 96
+          )
+        })
+        # ))
+      }
+
+      if (!is.null(x[[2]])) {
+        # print("renderPlotly:")
+        # print(system.time(
+        output$plot1_3d <- plotly::renderPlotly({
+          x[[2]]
+        })
+        # ))
+      }
     }) %>%
       finally(~ {
         prog$set(value = 10, message = "Done.", detail = "[Cell dimensional reduction plot]")
@@ -1734,15 +1710,19 @@ server <- function(input, output, session) {
       paste0("CellDimPlot-", format(Sys.time(), "%Y%m%d%H%M%S"), ".zip")
     },
     content = function(file) {
-      width <- get_attr(promisedData[["p1_dim"]], "width")
-      height <- get_attr(promisedData[["p1_dim"]], "height")
-      dpi <- get_attr(promisedData[["p1_dim"]], "dpi")
+      if (!is.null(promisedData[["p1_dim"]])) {
+        width <- get_attr(promisedData[["p1_dim"]], "width")
+        height <- get_attr(promisedData[["p1_dim"]], "height")
+        dpi <- get_attr(promisedData[["p1_dim"]], "dpi")
 
-      temp1 <- tempfile(pattern = "CellDimPlot-", fileext = ".png")
-      ggplot2::ggsave(filename = temp1, plot = promisedData[["p1_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
+        temp1 <- tempfile(pattern = "CellDimPlot-", fileext = ".png")
+        ggplot2::ggsave(filename = temp1, plot = promisedData[["p1_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
 
-      temp2 <- tempfile(pattern = "CellDimPlot-", fileext = ".pdf")
-      ggplot2::ggsave(filename = temp2, plot = promisedData[["p1_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
+        temp2 <- tempfile(pattern = "CellDimPlot-", fileext = ".pdf")
+        ggplot2::ggsave(filename = temp2, plot = promisedData[["p1_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
+      } else {
+        temp1 <- temp2 <- NULL
+      }
 
       if (!is.null(promisedData[["p1_3d"]])) {
         temp3 <- tempfile(pattern = "CellDimPlot3D-", fileext = ".html")
@@ -1773,9 +1753,10 @@ server <- function(input, output, session) {
     slots2 <- input$slots2
     features2 <- input$features2
     feature_area2 <- input$feature_area2
-    coExp2 <- input$coExp2
+    coExp2 <- as.logical(input$coExp2)
     scale2 <- input$scale2
-    raster2 <- input$raster2
+    raster2 <- as.logical(input$raster2)
+    dimension2 <- input$dimension2
     palette2 <- input$palette2
     theme2 <- input$theme2
     size2 <- input$size2
@@ -1812,27 +1793,29 @@ server <- function(input, output, session) {
 
         theme2 <- get(theme2, envir = asNamespace(themes[theme2]))
 
-        # print(">>> plot:")
-        # print(system.time(
-        p2_dim <- SCP::FeatureDimPlot(
-          srt = srt_tmp, features = features2, split.by = split2, reduction = reduction2, slot = "data", raster = raster2, pt.size = pt_size2,
-          calculate_coexp = coExp2, keep_scale = scale2, palette = palette2, theme_use = theme2,
-          ncol = ncol2, byrow = byrow2, force = TRUE
-        )
-        # ))
+        plot3d <- max(sapply(names(srt_tmp@reductions), function(r) dim(srt_tmp[[r]])[2])) >= 3 & dimension2 == "3D"
 
-        # print(">>> panel_fix:")
-        # print(system.time(
-        p2_dim <- SCP::panel_fix(SCP::slim_data(p2_dim), height = size2, units = "in", raster = panel_raster, BPPARAM = BPPARAM, verbose = FALSE)
-        # ))
-        attr(p2_dim, "dpi") <- 300
-        plot3d <- max(sapply(names(srt_tmp@reductions), function(r) dim(srt_tmp[[r]])[2])) >= 3
         if (isTRUE(plot3d)) {
+          p2_dim <- NULL
           p2_3d <- SCP::FeatureDimPlot3D(
             srt = srt_tmp, features = features2, reduction = reduction2, pt.size = pt_size2 * 2,
             calculate_coexp = coExp2, force = TRUE
           )
         } else {
+          # print(">>> plot:")
+          # print(system.time(
+          p2_dim <- SCP::FeatureDimPlot(
+            srt = srt_tmp, features = features2, split.by = split2, reduction = reduction2, slot = "data", raster = raster2, pt.size = pt_size2,
+            calculate_coexp = coExp2, keep_scale = scale2, palette = palette2, theme_use = theme2,
+            ncol = ncol2, byrow = byrow2, force = TRUE
+          )
+          # ))
+
+          # print(">>> panel_fix:")
+          # print(system.time(
+          p2_dim <- SCP::panel_fix(SCP::slim_data(p2_dim), height = size2, units = "in", raster = panel_raster, BPPARAM = BPPARAM, verbose = FALSE)
+          # ))
+          attr(p2_dim, "dpi") <- 300
           p2_3d <- NULL
         }
         return(list(p2_dim, p2_3d))
@@ -1843,26 +1826,72 @@ server <- function(input, output, session) {
     bindCache(
       input$dataset2, input$reduction2, input$split2, input$assays2, input$slots2,
       input$features2, input$feature_area2,
-      input$palette2, input$theme2, input$coExp2, input$scale2, input$raster2,
+      input$palette2, input$theme2, input$coExp2, input$scale2, input$raster2, input$dimension2,
       input$pt_size2, input$size2, input$ncol2, input$arrange2
     ) %>%
     bindEvent(input$submit2, ignoreNULL = FALSE, ignoreInit = FALSE)
 
+  observe({
+    prog <- Progress$new(min = 1, max = 10)
+    prog$set(value = 3, message = "Fetch data...", detail = "[Feature dimensional reduction plot]")
+    r2()$then(function(x) {
+      promisedData[["p2_dim"]] <- x[[1]]
+      promisedData[["p2_3d"]] <- x[[2]]
+      prog$set(value = 8, message = "Render plot...", detail = "[Feature dimensional reduction plot]")
+
+      if (!is.null(x[[1]])) {
+        width <- get_attr(x[[1]], "width")
+        height <- get_attr(x[[1]], "height")
+        dpi <- get_attr(x[[1]], "dpi")
+
+        # print("renderPlot:")
+        # print(system.time(
+        output$plot2 <- renderUI({
+          renderPlot(
+            {
+              x[[1]]
+            },
+            width = width * 96,
+            height = height * 96,
+            res = 96
+          )
+        })
+        # ))
+      }
+
+      if (!is.null(x[[2]])) {
+        # print("renderPlotly:")
+        # print(system.time(
+        output$plot2_3d <- plotly::renderPlotly({
+          x[[2]]
+        })
+        # ))
+      }
+    }) %>%
+      finally(~ {
+        prog$set(value = 10, message = "Done.", detail = "[Feature dimensional reduction plot]")
+        prog$close()
+      })
+  }) %>% bindEvent(input$submit2, ignoreNULL = FALSE, ignoreInit = FALSE)
+
   output$download2 <- downloadHandler(
     filename = function() {
       paste0("FeatureDimPlot-", format(Sys.time(), "%Y%m%d%H%M%S"), ".zip")
     },
     content = function(file) {
-      width <- get_attr(promisedData[["p2_dim"]], "width")
-      height <- get_attr(promisedData[["p2_dim"]], "height")
-      dpi <- get_attr(promisedData[["p2_dim"]], "dpi")
+      if (!is.null(promisedData[["p2_dim"]])) {
+        width <- get_attr(promisedData[["p2_dim"]], "width")
+        height <- get_attr(promisedData[["p2_dim"]], "height")
+        dpi <- get_attr(promisedData[["p2_dim"]], "dpi")
 
-      temp1 <- tempfile(pattern = "FeatureDimPlot-", fileext = ".png")
-      ggplot2::ggsave(filename = temp1, plot = promisedData[["p2_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
-
-      temp2 <- tempfile(pattern = "FeatureDimPlot-", fileext = ".pdf")
-      ggplot2::ggsave(filename = temp2, plot = promisedData[["p2_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
+        temp1 <- tempfile(pattern = "FeatureDimPlot-", fileext = ".png")
+        ggplot2::ggsave(filename = temp1, plot = promisedData[["p2_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
 
+        temp2 <- tempfile(pattern = "FeatureDimPlot-", fileext = ".pdf")
+        ggplot2::ggsave(filename = temp2, plot = promisedData[["p2_dim"]], width = width, height = height, units = "in", dpi = dpi, limitsize = FALSE)
+      } else {
+        temp1 <- temp2 <- NULL
+      }
       if (!is.null(promisedData[["p2_3d"]])) {
         temp3 <- tempfile(pattern = "FeatureDimPlot3D-", fileext = ".html")
         htmlwidgets::saveWidget(
@@ -1879,44 +1908,6 @@ server <- function(input, output, session) {
     contentType = "application/zip"
   )
 
-  observe({
-    prog <- Progress$new(min = 1, max = 10)
-    prog$set(value = 3, message = "Fetch data...", detail = "[Feature dimensional reduction plot]")
-    r2()$then(function(x) {
-      promisedData[["p2_dim"]] <- x[[1]]
-      promisedData[["p2_3d"]] <- x[[2]]
-      width <- get_attr(x[[1]], "width")
-      height <- get_attr(x[[1]], "height")
-      dpi <- get_attr(x[[1]], "dpi")
-
-      prog$set(value = 8, message = "Render plot...", detail = "[Feature dimensional reduction plot]")
-      # print("renderPlot:")
-      # print(system.time(
-      output$plot2 <- renderUI({
-        renderPlot(
-          {
-            x[[1]]
-          },
-          width = width * 96,
-          height = height * 96,
-          res = 96
-        )
-      })
-      # ))
-
-      # print("renderPlotly:")
-      # print(system.time(
-      output$plot2_3d <- plotly::renderPlotly({
-        x[[2]]
-      })
-      # ))
-    }) %>%
-      finally(~ {
-        prog$set(value = 10, message = "Done.", detail = "[Feature dimensional reduction plot]")
-        prog$close()
-      })
-  }) %>% bindEvent(input$submit2, ignoreNULL = FALSE, ignoreInit = FALSE)
-
   # submit3  ----------------------------------------------------------------
   r3 <- reactive({
     dataset3 <- input$dataset3
@@ -1936,8 +1927,8 @@ server <- function(input, output, session) {
     } else {
       split3 <- input$split3
     }
-    label3 <- input$label3
-    flip3 <- input$flip3
+    label3 <- as.logical(input$label3)
+    flip3 <- as.logical(input$flip3)
     palette3 <- input$palette3
     theme3 <- input$theme3
     labelsize3 <- input$labelsize3
@@ -2085,15 +2076,15 @@ server <- function(input, output, session) {
     feature_area4 <- input$feature_area4
     plotby4 <- input$plotby4
     fillby4 <- input$fillby4
-    coExp4 <- input$coExp4
-    stack4 <- input$stack4
-    flip4 <- input$flip4
-    addbox4 <- input$addbox4
-    addpoint4 <- input$addpoint4
-    addtrend4 <- input$addtrend4
+    coExp4 <- as.logical(input$coExp4)
+    stack4 <- as.logical(input$stack4)
+    flip4 <- as.logical(input$flip4)
+    addbox4 <- as.logical(input$addbox4)
+    addpoint4 <- as.logical(input$addpoint4)
+    addtrend4 <- as.logical(input$addtrend4)
     palette4 <- input$palette4
     theme4 <- input$theme4
-    sameylims4 <- input$sameylims4
+    sameylims4 <- as.logical(input$sameylims4)
     size4 <- input$size4
     ncol4 <- input$ncol4
     byrow4 <- input$arrange4
@@ -2237,7 +2228,7 @@ server <- function(input, output, session) {
 
   main_code <- readLines(textConnection(main_code))
   args <- mget(names(formals()))
-  args <- args[!names(args) %in% c("base_dir", "create_script", "style_script", "overwrite", "return_app")]
+  args <- args[!names(args) %in% c("base_dir", "create_script", "style_script", "overwrite")]
   args_code <- NULL
   for (varnm in names(args)) {
     args_code <- c(args_code, paste0(varnm, "=", deparse(args[[varnm]])))
@@ -2269,11 +2260,14 @@ server <- function(input, output, session) {
       BPPARAM = MulticoreParam(workers = plotting_workers)
     }",
     "page_theme <- bs_theme(bootswatch = 'zephyr')",
+    paste0("DataFile <- \"", DataFile_full, "\""),
+    paste0("MetaFile <- \"", MetaFile_full, "\""),
     main_code,
     "shinyApp(ui = ui, server = server)"
   )
   temp <- tempfile("SCExplorer")
   writeLines(app_code, temp)
+  source(temp)
   if (isTRUE(create_script)) {
     app_file <- paste0(base_dir, "/app.R")
     if (!file.exists(app_file) || isTRUE(overwrite)) {
@@ -2290,10 +2284,7 @@ server <- function(input, output, session) {
   }
   unlink(temp)
 
-  if (isTRUE(return_app)) {
-    app <- shiny::shinyAppDir(base_dir)
-    return(app)
-  } else {
-    return(invisible(NULL))
-  }
+  # app <- shiny::shinyAppFile(app_file)
+  app <- shiny::shinyApp(ui = ui, server = server)
+  return(app)
 }
diff --git a/R/SCP-plot.R b/R/SCP-plot.R
index 3dd0f246..57eb7c88 100644
--- a/R/SCP-plot.R
+++ b/R/SCP-plot.R
@@ -1359,7 +1359,7 @@ CellDimPlot <- function(srt, group.by, reduction = NULL, dims = c(1, 2), split.b
                         velocity_density = 1, velocity_smooth = 0.5, velocity_scale = 1, velocity_min_mass = 1, velocity_cutoff_perc = 5,
                         velocity_arrow_color = "black", velocity_arrow_angle = 20,
                         streamline_L = 5, streamline_minL = 1, streamline_res = 1, streamline_n = 15,
-                        streamline_width = c(0, 0.8), streamline_alpha = 1, streamline_color = NULL, streamline_palette = "RdYlBu", streamline_palcolor = NULL,
+                        streamline_width = c(0, 0.5), streamline_alpha = 1, streamline_color = NULL, streamline_palette = "RdYlBu", streamline_palcolor = NULL,
                         streamline_bg_color = "white", streamline_bg_stroke = 0.5,
                         hex = FALSE, hex.linewidth = 0.5, hex.count = TRUE, hex.bins = 50, hex.binwidth = NULL,
                         raster = NULL, raster.dpi = c(512, 512),
@@ -3644,6 +3644,8 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
     for (g in group.by) {
       plist_g <- plist[sapply(strsplit(names(plist), ":"), function(x) x[2]) == g]
       legend <- get_legend(plist_g[[1]])
+      title <- get_title(plist_g[[1]])
+      subtitle <- get_subtitle(plist_g[[1]])
       if (isTRUE(flip)) {
         lab <- textGrob(label = ifelse(is.null(ylab), "Expression", ylab), hjust = 0.5)
         plist_g <- lapply(seq_along(plist_g), FUN = function(i) {
@@ -3711,6 +3713,8 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
         gtable <- add_grob(gtable, lab, "left", clip = "off")
         gtable <- add_grob(gtable, legend, legend.position)
       }
+      gtable <- add_grob(gtable, subtitle, "top", adjust = 0, clip = "off")
+      gtable <- add_grob(gtable, title, "top", adjust = 0, clip = "off")
       gtable <- gtable_add_padding(gtable, unit(c(1, 1, 1, 1), units = "cm"))
       plot <- wrap_plots(gtable)
       plist_stack[[g]] <- plot
@@ -11117,9 +11121,6 @@ DynamicHeatmap <- function(srt, lineages, features = NULL, use_fitted = FALSE, b
   if (isTRUE(raster_by_magick)) {
     check_R("magick")
   }
-  if (is.null(features)) {
-    stop("No feature provided.")
-  }
   if (is.list(features)) {
     if (!is.null(names(features))) {
       feature_split <- rep(names(features), sapply(features, length))
@@ -14816,20 +14817,30 @@ as_gtable <- function(plot, ...) {
   }
 }
 
-get_legend <- function(plot) {
+get_component <- function(plot, pattern) {
   plot <- as_gtable(plot)
   grob_names <- plot$layout$name
   grobs <- plot$grobs
-  grobIndex <- which(grepl("guide-box", grob_names))
+  grobIndex <- which(grepl(pattern, grob_names))
   grobIndex <- grobIndex[1]
   matched_grobs <- grobs[[grobIndex]]
   return(matched_grobs)
 }
+get_legend <- function(plot) {
+  get_component(plot, "guide-box")
+}
+get_title <- function(plot) {
+  get_component(plot, "^title")
+}
+get_subtitle <- function(plot) {
+  get_component(plot, "^subtitle")
+}
 
 #' @importFrom grid is.grob grobWidth grobHeight
 #' @importFrom gtable is.gtable gtable_add_rows gtable_add_cols gtable_add_grob
-add_grob <- function(gtable, grob, position = c("top", "bottom", "left", "right", "none"), space = NULL, clip = "on") {
+add_grob <- function(gtable, grob, position = c("top", "bottom", "left", "right", "none"), adjust = 0.5, space = NULL, clip = "on") {
   position <- match.arg(position)
+
   if (position == "none" || is.null(grob)) {
     return(gtable)
   }
@@ -14852,19 +14863,32 @@ add_grob <- function(gtable, grob, position = c("top", "bottom", "left", "right"
 
   if (position == "top") {
     gtable <- gtable_add_rows(gtable, space, 0)
-    gtable <- gtable_add_grob(gtable, grob, t = 1, l = mean(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "l"]), clip = clip)
+    gtable <- gtable_add_grob(gtable, grob, t = 1, l = adjust_pos(gtable, direction = "horizontal", adjust = adjust), clip = clip)
   }
   if (position == "bottom") {
     gtable <- gtable_add_rows(gtable, space, -1)
-    gtable <- gtable_add_grob(gtable, grob, t = dim(gtable)[1], l = mean(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "l"]), clip = clip)
+    gtable <- gtable_add_grob(gtable, grob, t = dim(gtable)[1], l = adjust_pos(gtable, direction = "horizontal", adjust = adjust), clip = clip)
   }
   if (position == "left") {
     gtable <- gtable_add_cols(gtable, space, 0)
-    gtable <- gtable_add_grob(gtable, grob, t = mean(gtable$layout[grep("panel", gtable$layout$name), "t"]), l = 1, clip = clip)
+    gtable <- gtable_add_grob(gtable, grob, t = adjust_pos(gtable, direction = "vertical", adjust = adjust), l = 1, clip = clip)
   }
   if (position == "right") {
     gtable <- gtable_add_cols(gtable, space, -1)
-    gtable <- gtable_add_grob(gtable, grob, t = mean(gtable$layout[grep("panel", gtable$layout$name), "t"]), l = dim(gtable)[2], clip = clip)
+    gtable <- gtable_add_grob(gtable, grob, t = adjust_pos(gtable, direction = "vertical", adjust = adjust), l = dim(gtable)[2], clip = clip)
   }
   return(gtable)
 }
+
+adjust_pos <- function(gtable, direction = c("horizontal", "vertical"), adjust = 0.5) {
+  direction <- match.arg(direction)
+  if (direction == "horizontal") {
+    hmax <- max(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "l"])
+    hmin <- min(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "l"])
+    return(hmin + adjust * (hmax - hmin))
+  } else {
+    vmax <- max(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "t"])
+    vmin <- min(gtable$layout[grepl(pattern = "panel", x = gtable$layout$name), "t"])
+    return(vmin + adjust * (vmax - vmin))
+  }
+}
diff --git a/R/utils.R b/R/utils.R
index 03651e32..98ddf74c 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -829,19 +829,64 @@ unnest <- function(data, cols, keep_empty = FALSE) {
 #' @importFrom Matrix as.matrix
 #' @export
 as_matrix <- function(x) {
-  if (!inherits(matrix, "dgCMatrix")) {
+  if (!inherits(x, "dgCMatrix")) {
     return(as.matrix(x))
   } else {
     row_pos <- x@i
     col_pos <- findInterval(seq_along(x@x) - 1, x@p[-1])
-    out <- asMatrix(rp = row_pos, cp = col_pos, z = x@x, nrows = x@Dim[1], ncols = x@Dim[2])
+    out <- asMatrixC(rp = row_pos, cp = col_pos, z = x@x, nrows = x@Dim[1], ncols = x@Dim[2])
     attr(out, "dimnames") <- list(x@Dimnames[[1]], x@Dimnames[[2]])
     return(out)
   }
 }
 
+#' Merge Rows of a Matrix by Groupings
+#'
+#' MergeRows is a function that sum the values of rows based on specified groupings.
+#'
+#' @param x A matrix.
+#' @param groupings A vector specifying the groupings for merging rows.
+#' @useDynLib SCP
+#' @examples
+#' library(Matrix)
+#' ncells <- 3000
+#' nfeatures <- 1000
+#' expressed <- 500
+#' n <- ncells * expressed
+#' dimnames <- list(paste0("feature", seq_len(nfeatures)), paste0("cell", seq_len(ncells)))
+#'
+#' sparse <- sparseMatrix(
+#'   i = sample(seq_len(nfeatures), size = n, replace = TRUE),
+#'   j = sample(seq_len(ncells), size = n, replace = TRUE),
+#'   x = sample(1:10, size = n, replace = TRUE),
+#'   dimnames = dimnames
+#' )
+#' dup_index <- sample(seq_len(nfeatures), size = 50)
+#' rownames(sparse)[dup_index] <- paste0("duplicated", 1:10)
+#' table(rownames(sparse)[dup_index])
+#'
+#' dense <- matrix(
+#'   data = sample(1:10, size = nfeatures * ncells, replace = TRUE),
+#'   nrow = nfeatures, ncol = ncells,
+#'   dimnames = dimnames
+#' )
+#' dup_index <- sample(seq_len(nfeatures), size = 50)
+#' rownames(dense)[dup_index] <- paste0("duplicated", 1:10)
+#' table(rownames(dense)[dup_index])
+#'
+#' system.time(MergeRows(dense, rownames(dense)))
+#' system.time(MergeRows(sparse, rownames(sparse)))
+#' system.time(aggregate(dense, by = list(rownames(dense)), FUN = sum))
+#' @export
+MergeRows <- function(x, groupings) {
+  stopifnot(inherits(x, c("dgCMatrix", "matrix")))
+  out <- MergeRowsC(x, groupings)
+  colnames(out) <- colnames(x)
+  rownames(out) <- unique(groupings)
+  return(out)
+}
+
 #' Capitalizes the characters
-#' Making the first letter uppercase
 #'
 #' @examples
 #' x <- c("dna methylation", "rRNA processing", "post-Transcriptional gene silencing")
diff --git a/man/CellDimPlot.Rd b/man/CellDimPlot.Rd
index 341b7640..eb2096f1 100644
--- a/man/CellDimPlot.Rd
+++ b/man/CellDimPlot.Rd
@@ -96,7 +96,7 @@ CellDimPlot(
   streamline_minL = 1,
   streamline_res = 1,
   streamline_n = 15,
-  streamline_width = c(0, 0.8),
+  streamline_width = c(0, 0.5),
   streamline_alpha = 1,
   streamline_color = NULL,
   streamline_palette = "RdYlBu",
diff --git a/man/MergeRows.Rd b/man/MergeRows.Rd
new file mode 100644
index 00000000..af9d83ea
--- /dev/null
+++ b/man/MergeRows.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{MergeRows}
+\alias{MergeRows}
+\title{Merge Rows of a Matrix by Groupings}
+\usage{
+MergeRows(x, groupings)
+}
+\arguments{
+\item{x}{A matrix.}
+
+\item{groupings}{A vector specifying the groupings for merging rows.}
+}
+\description{
+MergeRows is a function that sum the values of rows based on specified groupings.
+}
diff --git a/man/PrepareSCExplorer.Rd b/man/PrepareSCExplorer.Rd
index b9e9822e..bca6ff86 100644
--- a/man/PrepareSCExplorer.Rd
+++ b/man/PrepareSCExplorer.Rd
@@ -7,8 +7,6 @@
 PrepareSCExplorer(
   object,
   base_dir = "SCExplorer",
-  DataFile = "Data.hdf5",
-  MetaFile = "Meta.hdf5",
   assays = "RNA",
   slots = c("counts", "data"),
   ignore_nlevel = 100,
@@ -23,10 +21,6 @@ PrepareSCExplorer(
 
 \item{base_dir}{The base directory where the SCExplorer hdf5 files will be written. Default is "SCExplorer".}
 
-\item{DataFile}{Path to the output data file. If not provided, the file will be named "Data.hdf5" in the current directory.}
-
-\item{MetaFile}{Path to the output meta file. If not provided, the file will be named "Meta.hdf5" in the current directory.}
-
 \item{assays}{Character vector specifying the assays to include in the data file. Default is "RNA".}
 
 \item{slots}{Character vector specifying the slots to include in the data file. Default is "data".}
diff --git a/man/RunSCExplorer.Rd b/man/RunSCExplorer.Rd
index 3af7b87a..26acc27a 100644
--- a/man/RunSCExplorer.Rd
+++ b/man/RunSCExplorer.Rd
@@ -6,8 +6,6 @@
 \usage{
 RunSCExplorer(
   base_dir = "SCExplorer",
-  DataFile = "Data.hdf5",
-  MetaFile = "Meta.hdf5",
   title = "SCExplorer",
   initial_dataset = NULL,
   initial_reduction = NULL,
@@ -27,17 +25,12 @@ RunSCExplorer(
   plotting_workers = 8,
   create_script = TRUE,
   style_script = require("styler", quietly = TRUE),
-  overwrite = FALSE,
-  return_app = TRUE
+  overwrite = FALSE
 )
 }
 \arguments{
 \item{base_dir}{A string. The base directory of the SCExplorer app. Default is "SCExplorer".}
 
-\item{DataFile}{A string. The name of the HDF5 file that stores data matrices for each dataset. Default is "Data.hdf5".}
-
-\item{MetaFile}{A string. The name of the HDF5 file that stores metadata for each dataset. Default is "Meta.hdf5".}
-
 \item{title}{A string. The title of the SCExplorer app. Default is "SCExplorer".}
 
 \item{initial_dataset}{A string. The initial dataset to be loaded into the app. Default is NULL.}
@@ -78,7 +71,9 @@ RunSCExplorer(
 
 \item{overwrite}{A logical. Whether to overwrite existing files. Default is FALSE.}
 
-\item{return_app}{A logical. Whether to return the SCExplorer app. Default is TRUE.}
+\item{DataFile}{A string. The name of the HDF5 file that stores data matrices for each dataset. Default is "Data.hdf5".}
+
+\item{MetaFile}{A string. The name of the HDF5 file that stores metadata for each dataset. Default is "Meta.hdf5".}
 }
 \description{
 RunSCExplorer
diff --git a/man/capitalize.Rd b/man/capitalize.Rd
index d9b5d9fc..7e8dab23 100644
--- a/man/capitalize.Rd
+++ b/man/capitalize.Rd
@@ -2,8 +2,7 @@
 % Please edit documentation in R/utils.R
 \name{capitalize}
 \alias{capitalize}
-\title{Capitalizes the characters
-Making the first letter uppercase}
+\title{Capitalizes the characters}
 \usage{
 capitalize(x, force_tolower = FALSE)
 }
@@ -14,7 +13,6 @@ capitalize(x, force_tolower = FALSE)
 }
 \description{
 Capitalizes the characters
-Making the first letter uppercase
 }
 \examples{
 x <- c("dna methylation", "rRNA processing", "post-Transcriptional gene silencing")
diff --git a/src/MergeRowsC.cpp b/src/MergeRowsC.cpp
new file mode 100644
index 00000000..41ed6db4
--- /dev/null
+++ b/src/MergeRowsC.cpp
@@ -0,0 +1,70 @@
+#include <RcppArmadillo.h>
+#include <unordered_map>
+#include <vector>
+#include <string>
+
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+arma::mat MergeRowsForDense(const arma::mat& matrix, const std::vector<std::string>& groupings) {
+  
+  int nrow = matrix.n_rows;
+  int ncol = matrix.n_cols;
+  std::unordered_set<std::string> unique_groupings(groupings.begin(), groupings.end());
+  int n = unique_groupings.size();
+  
+  if (nrow != static_cast<int>(groupings.size())) stop("groupings must be the length of nrow(matrix)");
+  
+  std::unordered_map<std::string, int> rowIndices(n); //reserve space when declaring
+  
+  arma::mat result(n, ncol, arma::fill::zeros); // reserve space for result matrix
+  
+  int index = 0;
+  for (int i = 0; i < nrow; ++i) {
+    const std::string& row = groupings[i]; // avoid copying string
+    auto it = rowIndices.find(row);
+    if (it != rowIndices.end()) {
+      result.row(it->second) += matrix.row(i);
+    } else {
+      rowIndices.emplace(row, index);
+      result.row(index) = matrix.row(i);
+      ++index;
+    }
+  }
+  
+  return result;
+}
+
+// [[Rcpp::export]]
+arma::sp_mat MergeRowsForSparse(const arma::sp_mat& matrix, const std::vector<std::string>& groupings) {
+  // Create a mapping matrix
+  std::unordered_map<std::string, int> unique_groupings;
+  for (const auto &group : groupings) {
+    if (unique_groupings.find(group) == unique_groupings.end()) {
+      unique_groupings[group] = unique_groupings.size();
+    }
+  }
+  arma::sp_mat mapping(groupings.size(), unique_groupings.size());
+  for (size_t i = 0; i < groupings.size(); ++i) {
+    mapping(i, unique_groupings[groupings[i]]) = 1;
+  }
+  
+  // Multiply mapping by matrix
+  arma::sp_mat result = mapping.t() * matrix;
+  
+  return result;
+}
+
+// [[Rcpp::export]]
+SEXP MergeRowsC(const SEXP& matrix, const std::vector<std::string>& groupings) {
+  if(Rf_inherits(matrix, "dgCMatrix")){    
+    arma::sp_mat mat = MergeRowsForSparse(as<arma::sp_mat>(S4(matrix)), groupings);  
+    return(wrap(mat));
+  } else if (Rf_isMatrix(matrix)){
+    arma::mat mat =  MergeRowsForDense(as<arma::mat>(matrix), groupings);    
+    return(wrap(mat));  
+  } else {    
+    stop("Unsupported matrix type");    
+  }  
+}  
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index ec7932f3..8c2c824d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,6 +1,7 @@
 // Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#include <RcppArmadillo.h>
 #include <Rcpp.h>
 
 using namespace Rcpp;
@@ -10,9 +11,45 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
-// asMatrix
-NumericMatrix asMatrix(NumericVector rp, NumericVector cp, NumericVector z, int nrows, int ncols);
-RcppExport SEXP _SCP_asMatrix(SEXP rpSEXP, SEXP cpSEXP, SEXP zSEXP, SEXP nrowsSEXP, SEXP ncolsSEXP) {
+// MergeRowsForDense
+arma::mat MergeRowsForDense(const arma::mat& matrix, const std::vector<std::string>& groupings);
+RcppExport SEXP _SCP_MergeRowsForDense(SEXP matrixSEXP, SEXP groupingsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type matrix(matrixSEXP);
+    Rcpp::traits::input_parameter< const std::vector<std::string>& >::type groupings(groupingsSEXP);
+    rcpp_result_gen = Rcpp::wrap(MergeRowsForDense(matrix, groupings));
+    return rcpp_result_gen;
+END_RCPP
+}
+// MergeRowsForSparse
+arma::sp_mat MergeRowsForSparse(const arma::sp_mat& matrix, const std::vector<std::string>& groupings);
+RcppExport SEXP _SCP_MergeRowsForSparse(SEXP matrixSEXP, SEXP groupingsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::sp_mat& >::type matrix(matrixSEXP);
+    Rcpp::traits::input_parameter< const std::vector<std::string>& >::type groupings(groupingsSEXP);
+    rcpp_result_gen = Rcpp::wrap(MergeRowsForSparse(matrix, groupings));
+    return rcpp_result_gen;
+END_RCPP
+}
+// MergeRowsC
+SEXP MergeRowsC(const SEXP& matrix, const std::vector<std::string>& groupings);
+RcppExport SEXP _SCP_MergeRowsC(SEXP matrixSEXP, SEXP groupingsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const SEXP& >::type matrix(matrixSEXP);
+    Rcpp::traits::input_parameter< const std::vector<std::string>& >::type groupings(groupingsSEXP);
+    rcpp_result_gen = Rcpp::wrap(MergeRowsC(matrix, groupings));
+    return rcpp_result_gen;
+END_RCPP
+}
+// asMatrixC
+NumericMatrix asMatrixC(NumericVector rp, NumericVector cp, NumericVector z, int nrows, int ncols);
+RcppExport SEXP _SCP_asMatrixC(SEXP rpSEXP, SEXP cpSEXP, SEXP zSEXP, SEXP nrowsSEXP, SEXP ncolsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -21,13 +58,16 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP);
     Rcpp::traits::input_parameter< int >::type nrows(nrowsSEXP);
     Rcpp::traits::input_parameter< int >::type ncols(ncolsSEXP);
-    rcpp_result_gen = Rcpp::wrap(asMatrix(rp, cp, z, nrows, ncols));
+    rcpp_result_gen = Rcpp::wrap(asMatrixC(rp, cp, z, nrows, ncols));
     return rcpp_result_gen;
 END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_SCP_asMatrix", (DL_FUNC) &_SCP_asMatrix, 5},
+    {"_SCP_MergeRowsForDense", (DL_FUNC) &_SCP_MergeRowsForDense, 2},
+    {"_SCP_MergeRowsForSparse", (DL_FUNC) &_SCP_MergeRowsForSparse, 2},
+    {"_SCP_MergeRowsC", (DL_FUNC) &_SCP_MergeRowsC, 2},
+    {"_SCP_asMatrixC", (DL_FUNC) &_SCP_asMatrixC, 5},
     {NULL, NULL, 0}
 };
 
diff --git a/src/asMatrix.cpp b/src/asMatrixC.cpp
similarity index 90%
rename from src/asMatrix.cpp
rename to src/asMatrixC.cpp
index 0897a1ba..03cf002f 100644
--- a/src/asMatrix.cpp
+++ b/src/asMatrixC.cpp
@@ -2,7 +2,7 @@
 using namespace Rcpp;    
 
 // [[Rcpp::export]]    
-NumericMatrix asMatrix(NumericVector rp,
+NumericMatrix asMatrixC(NumericVector rp,
                        NumericVector cp,
                        NumericVector z,
                        int nrows,