diff --git a/.chplcheckignore b/.chplcheckignore
index 3929ed8da8e..0353f30fc43 100644
--- a/.chplcheckignore
+++ b/.chplcheckignore
@@ -15,6 +15,7 @@ CommDiagnosticsMsg.chpl
 compat
 ConcatenateMsg.chpl
 CSVMsg.chpl
+CustomCopyAggregation.chpl
 DataFrameIndexingMsg.chpl
 deprecated
 DynamicSort.chpl
diff --git a/.flake8 b/.flake8
index ce46040e01c..08d57ba88a7 100644
--- a/.flake8
+++ b/.flake8
@@ -58,7 +58,7 @@ per-file-ignores =
     arkouda/pandas/join.py: DOC105,DOC105,DOC106,DOC107,DOC203,DOC501,DOC503
     arkouda/pandas/match.py: DOC101,DOC103,DOC105,DOC106,DOC107,DOC203,DOC501,DOC502,DOC503
     arkouda/pandas/matcher.py: DOC101,DOC103,DOC501,DOC503,DOC603
-    arkouda/scipy/sparrayclass.py: DOC101,DOC103,DOC107,DOC503
+    arkouda/scipy/sparrayclass.py: DOC101,DOC103,DOC107,DOC203,DOC503
     arkouda/pandas/series.py: DOC101,DOC103,DOC105,DOC106,DOC107,DOC203,DOC201,DOC501,DOC502,DOC503,DOC601,DOC603
     arkouda/plotting.py: DOC105,DOC203,DOC503,DOC501
     arkouda/scipy/_stats_py.py: DOC105,DOC106,DOC107,DOC203,DOC501,DOC503
diff --git a/arkouda-env-dev.yml b/arkouda-env-dev.yml
index d49a2597221..c610f849a00 100644
--- a/arkouda-env-dev.yml
+++ b/arkouda-env-dev.yml
@@ -41,6 +41,7 @@ dependencies:
   - pytest-benchmark>=4.0.0
   - mathjax
   - pandas-stubs
+  - scipy-stubs
   - types-python-dateutil
   - pre-commit
   - pytest-subtests
diff --git a/arkouda/scipy/sparrayclass.py b/arkouda/scipy/sparrayclass.py
index 45df14d053c..892fdb34406 100644
--- a/arkouda/scipy/sparrayclass.py
+++ b/arkouda/scipy/sparrayclass.py
@@ -142,7 +142,7 @@ def to_pdarray(self) -> List[pdarray]:
         if dtype_name not in NumericDTypes:
             raise TypeError(f"unsupported dtype {dtype}")
         response_arrays = generic_msg(
-            cmd=f"sparse_to_pdarrays<{self.dtype},{self.layout}>", args={"matrix": self}
+            cmd=f"sparse_to_pdarrays<{self.dtype},{self.layout}>", args={"matrix": self.name}
         )
         array_list = create_pdarrays(type_cast(str, response_arrays))
         return array_list
@@ -163,6 +163,41 @@ def fill_vals(self, a: pdarray):
             args={"matrix": self, "vals": a},
         )
 
+    def to_scipy_sparse(self):
+        """
+        Convert the Arkouda sparse array to an equivalent SciPy sparse array.
+
+        Returns
+        -------
+        scipy.sparse.sparray
+            A SciPy sparse array with the same shape, values, and logical layout.
+
+        Notes
+        -----
+        SciPy sparse array classes require SciPy versions that provide
+        ``coo_array``, ``csr_array``, and ``csc_array``. If those are not
+        available, this method falls back to the corresponding sparse matrix
+        classes.
+        """
+        import scipy.sparse as sp
+
+        rows, cols, vals = (arr.to_ndarray() for arr in self.to_pdarray())
+
+        if hasattr(sp, "coo_array"):
+            coo = sp.coo_array((vals, (rows, cols)), shape=tuple(self.shape))
+            if self.layout == "CSR":
+                return sp.csr_array(coo)
+            if self.layout == "CSC":
+                return sp.csc_array(coo)
+            return coo
+
+        coo = sp.coo_matrix((vals, (rows, cols)), shape=tuple(self.shape))
+        if self.layout == "CSR":
+            return coo.tocsr()
+        if self.layout == "CSC":
+            return coo.tocsc()
+        return coo
+
 
 # creates sparray object
 #   only after:
diff --git a/benchmarks/graph_infra/arkouda.graph b/benchmarks/graph_infra/arkouda.graph
index 1415f6d6d30..c6f78591c69 100644
--- a/benchmarks/graph_infra/arkouda.graph
+++ b/benchmarks/graph_infra/arkouda.graph
@@ -159,3 +159,15 @@ graphkeys: Noop ops/s
 files: noop.dat
 graphtitle: Noop Performance
 ylabel: Performance (ops/s)
+
+perfkeys: Average CSR time =, Average CSC time =
+graphkeys: CSR time, CSC time
+files: sparse.dat, sparse.dat
+graphtitle: Sparse Matrix Creation Time
+ylabel: Time (Seconds)
+
+perfkeys: Average CSCxCSR time =
+graphkeys: CSCxCSR time
+files: sparse.dat
+graphtitle: Sparse Matrix Multiplication Time
+ylabel: Time (Seconds)
diff --git a/benchmarks/graph_infra/sparse.perfkeys b/benchmarks/graph_infra/sparse.perfkeys
new file mode 100644
index 00000000000..b0a0819881d
--- /dev/null
+++ b/benchmarks/graph_infra/sparse.perfkeys
@@ -0,0 +1,3 @@
+Average CSR time =
+Average CSC time =
+Average CSCxCSR time =
diff --git a/benchmarks/run_benchmarks.py b/benchmarks/run_benchmarks.py
index 1b5cf4d2b9c..e7622ad300d 100755
--- a/benchmarks/run_benchmarks.py
+++ b/benchmarks/run_benchmarks.py
@@ -52,6 +52,7 @@
     "str-gather",
     "str-in1d",
     "substring_search",
+    "sparse",
     "split",
     "sort-cases",
     "multiIO",
diff --git a/benchmarks/sparse.py b/benchmarks/sparse.py
new file mode 100644
index 00000000000..c927018d1f4
--- /dev/null
+++ b/benchmarks/sparse.py
@@ -0,0 +1,194 @@
+#!/usr/bin/env python3
+
+import arkouda as ak
+import argparse
+import time
+import numpy as np
+from scipy.sparse import coo_array, csr_matrix, find
+from arkouda.scipy.sparsematrix import create_sparse_matrix, sparse_matrix_matrix_mult
+
+# TODO: Add support for 'float64'
+TYPES = ("int64",)
+
+
+def compare_scipy(left, right, rtol=1e-9, atol=0.0, equal_nan=False):
+    if left.shape != right.shape:
+        return False
+
+    lr, lc, lv = find(left)
+    rr, rc, rv = find(right)
+
+    if len(lr) != len(rr) or not np.all(lr == rr):
+        return False
+    if len(lc) != len(rc) or not np.all(lc == rc):
+        return False
+    if not np.allclose(lv, rv, rtol=rtol, atol=atol, equal_nan=equal_nan):
+        return False
+
+    return True
+
+
+def time_ak_sparse(N, trials, dtype, seed):
+    print(">>> arkouda {} sparse".format(dtype))
+    cfg = ak.get_config()
+    print("numLocales = {}, numNodes {}, N = {:,}".format(cfg["numLocales"], cfg["numNodes"], N))
+
+    nnz = N * 10
+
+    rows = ak.randint(1, N - 1, nnz, seed=seed)
+    cols = ak.randint(1, N - 1, nnz, seed=seed * 2)
+    rows, cols = ak.unique([rows, cols])
+
+    if dtype == "int64":
+        vals = ak.randint(1, len(rows), len(rows), seed=seed * 3)
+    elif dtype == "float64":
+        vals = ak.uniform(len(rows), seed=seed * 3) + 0.5
+
+    csr_times = []
+    csc_times = []
+    multiplication_times = []
+
+    for i in range(trials):
+        start = time.time()
+        csr = create_sparse_matrix(N, rows, cols, vals, layout="CSR")
+        csr_times.append(time.time() - start)
+
+        start = time.time()
+        csc = create_sparse_matrix(N, cols, rows, vals, layout="CSC")
+        csc_times.append(time.time() - start)
+
+        start = time.time()
+        result = sparse_matrix_matrix_mult(csc, csr)
+        multiplication_times.append(time.time() - start)
+
+    print("Average CSR time = {:.4f} seconds".format(np.mean(csr_times)))
+    print("Average CSC time = {:.4f} seconds".format(np.mean(csc_times)))
+    print("Average CSCxCSR time = {:.4f} seconds".format(np.mean(multiplication_times)))
+
+
+def time_np_sparse(N, trials, dtype, seed):
+    print(">>> numpy {} sparse".format(dtype))
+    print("N = {:,}".format(N))
+
+    nnz = N * 10
+
+    rows = ak.randint(1, N - 1, nnz, seed=seed)
+    cols = ak.randint(1, N - 1, nnz, seed=seed * 2)
+    rows, cols = ak.unique([rows, cols])
+
+    if dtype == "int64":
+        vals = ak.randint(1, len(rows), len(rows), seed=seed * 3)
+    elif dtype == "float64":
+        vals = ak.uniform(len(rows), seed=seed * 3) + 0.5
+
+    rows = rows.to_ndarray()
+    cols = cols.to_ndarray()
+    vals = vals.to_ndarray()
+
+    csr_times = []
+    csc_times = []
+    multiplication_times = []
+
+    for i in range(trials):
+        start = time.time()
+        csr = coo_array((vals, (rows, cols)), shape=(N, N)).tocsr()
+        csr_times.append(time.time() - start)
+
+        start = time.time()
+        csc = coo_array((vals, (cols, rows)), shape=(N, N)).tocsc()
+        csc_times.append(time.time() - start)
+
+        start = time.time()
+        result = csc.dot(csr).tocsr()
+        multiplication_times.append(time.time() - start)
+
+    print("Average CSR time = {:.4f} seconds".format(np.mean(csr_times)))
+    print("Average CSC time = {:.4f} seconds".format(np.mean(csc_times)))
+    print("Average CSCxCSR time = {:.4f} seconds".format(np.mean(multiplication_times)))
+
+
+def check_correctness(dtype):
+    seed = 1234
+    N = 10_000
+    nnz = N * 10
+
+    rows = ak.randint(1, N - 1, nnz, seed=seed)
+    cols = ak.randint(1, N - 1, nnz, seed=seed * 2)
+    rows, cols = ak.unique([rows, cols])
+
+    if dtype == "int64":
+        vals = ak.randint(1, len(rows), len(rows), seed=seed * 3)
+    elif dtype == "float64":
+        vals = ak.uniform(len(rows), seed=seed * 3) + 0.5
+
+    csr = create_sparse_matrix(N, rows, cols, vals, layout="CSR")
+    scipy_csr = coo_array(
+        (vals.to_ndarray(), (rows.to_ndarray(), cols.to_ndarray())), shape=(N, N)
+    ).tocsr()
+    ak_csr = csr.to_scipy_sparse()
+    assert compare_scipy(scipy_csr, ak_csr), "CSR matrices do not match"
+
+    csc = create_sparse_matrix(N, cols, rows, vals, layout="CSC")
+    scipy_csc = coo_array(
+        (vals.to_ndarray(), (cols.to_ndarray(), rows.to_ndarray())), shape=(N, N)
+    ).tocsc()
+    ak_csc = csc.to_scipy_sparse()
+    assert compare_scipy(scipy_csc, ak_csc), "CSC matrices do not match"
+
+    result = sparse_matrix_matrix_mult(csc, csr)
+    scipy_result = scipy_csc.dot(scipy_csr).tocsr()
+    ak_result = result.to_scipy_sparse()
+    assert compare_scipy(scipy_result, ak_result), "Multiplication results do not match"
+
+
+def create_parser():
+    parser = argparse.ArgumentParser(description="Benchmark sparse matrix creation and multiplication.")
+    parser.add_argument("hostname", type=str, help="Name of the Arkouda server")
+    parser.add_argument("port", type=int, help="Port of the Arkouda server")
+    parser.add_argument("-n", "--size", type=int, default=(10**6), help="Size of the sparse matrices")
+    parser.add_argument("-t", "--trials", type=int, default=3, help="Number of trials for benchmarking")
+    parser.add_argument(
+        "-d", "--dtype", default="int64", help="Dtype of array ({})".format(", ".join(TYPES))
+    )
+    parser.add_argument(
+        "--numpy",
+        default=False,
+        action="store_true",
+        help="Run the same operation in NumPy to compare performance.",
+    )
+    parser.add_argument(
+        "--correctness-only",
+        default=False,
+        action="store_true",
+        help="Only check correctness, not performance.",
+    )
+    parser.add_argument(
+        "-s", "--seed", default=1234, type=int, help="Value to initialize random number generator"
+    )
+    return parser
+
+
+if __name__ == "__main__":
+    import sys
+
+    args = create_parser().parse_args()
+
+    if args.dtype not in TYPES:
+        raise ValueError("Dtype must be {}, not {}".format("/".join(TYPES), args.dtype))
+    ak.verbose = False
+    ak.connect(args.hostname, args.port)
+
+    if args.correctness_only:
+        for dtype in TYPES:
+            check_correctness(dtype)
+        sys.exit(0)
+
+    print("N = {:,}".format(args.size))
+    print("number of trials = ", args.trials)
+
+    time_ak_sparse(args.size, args.trials, args.dtype, args.seed)
+
+    if args.numpy:
+        time_np_sparse(args.size, args.trials, args.dtype, args.seed)
+
+    sys.exit(0)
diff --git a/docker/ci/almalinux-with-arkouda-deps/requirements.txt b/docker/ci/almalinux-with-arkouda-deps/requirements.txt
index 71e3893ee64..30b9b03f514 100644
--- a/docker/ci/almalinux-with-arkouda-deps/requirements.txt
+++ b/docker/ci/almalinux-with-arkouda-deps/requirements.txt
@@ -31,6 +31,7 @@ pytest-json-report
 pytest-benchmark
 mathjax
 pandas-stubs
+scipy-stubs
 types-python-dateutil
 blosc2>=2.3.0
 numexpr>=2.6.2
diff --git a/pydoc/requirements.txt b/pydoc/requirements.txt
index 7f69d04c073..8188c3745a3 100644
--- a/pydoc/requirements.txt
+++ b/pydoc/requirements.txt
@@ -41,4 +41,5 @@ pytest-json-report
 pytest-benchmark>=4.0.0
 mathjax
 pandas-stubs
+scipy-stubs
 types-python-dateutil
diff --git a/pydoc/setup/REQUIREMENTS.md b/pydoc/setup/REQUIREMENTS.md
index 5a5894f1aff..b4f7f962eb1 100644
--- a/pydoc/setup/REQUIREMENTS.md
+++ b/pydoc/setup/REQUIREMENTS.md
@@ -62,6 +62,7 @@ The dependencies listed here are only required if you will be doing development
 - `pytest-benchmark>=4.0.0`
 - `mathjax`
 - `pandas-stubs`
+- `scipy-stubs`
 - `types-python-dateutil`
 
 ### Installing/Updating Python Dependencies
diff --git a/pyproject.toml b/pyproject.toml
index fb943e54f47..1d5f52fe997 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -57,6 +57,7 @@ dev = [
   "sphinx-design",
   "sphinx-autodoc-typehints",
   "pandas-stubs",
+  "scipy-stubs",
   "types-python-dateutil",
   "ipython",
   "pre-commit",
diff --git a/src/CustomCopyAggregation.chpl b/src/CustomCopyAggregation.chpl
new file mode 100644
index 00000000000..81985c7d07f
--- /dev/null
+++ b/src/CustomCopyAggregation.chpl
@@ -0,0 +1,229 @@
+module CustomCopyAggregation {
+  use AggregationPrimitives;
+  use ChplConfig;
+  use CTypes;
+  use CopyAggregation;
+
+  private config param verboseAggregation = false;
+
+  private param defaultBuffSize = if CHPL_TARGET_PLATFORM == "hpe-cray-ex" then 1024
+                                                                  else if CHPL_COMM == "ugni" then 4096
+                                                                  else 8192;
+
+  private const yieldFrequency = getEnvInt("CHPL_AGGREGATION_YIELD_FREQUENCY", 1024);
+  private const dstBuffSize = getEnvInt("CHPL_AGGREGATION_DST_BUFF_SIZE", defaultBuffSize);
+  private const srcBuffSize = getEnvInt("CHPL_AGGREGATION_SRC_BUFF_SIZE", defaultBuffSize);
+  private config param aggregate = CHPL_COMM != "none";
+
+  record remoteHandler {
+    type t;
+    var original: _to_nilable(t);
+    var localHandler: _to_nilable(_to_unmanaged(original!.sourceCopy().type));
+    var loc: int;
+
+    proc init(type t) {
+      this.t = t;
+      this.localHandler = nil;
+    }
+
+    proc init(type t, ref original, loc) {
+      this.t = t;
+      this.original = original;
+      this.localHandler = nil;
+      this.loc = loc;
+    }
+
+    proc deinit() {
+      delete localHandler;
+    }
+
+    inline proc ref get() ref {
+      if localHandler == nil {
+        on Locales[loc] {
+          localHandler = original!.sourceCopy();
+        }
+      }
+      return localHandler;
+    }
+  }
+
+  /*
+    Aggregates ``copy(ref dst, src)``. Optimized for when src is local.
+    Not parallel safe and is expected to be created on a per-task basis.
+    High memory usage since there are per-destination buffers.
+  */
+  record CustomDstAggregator {
+    type elemType;
+    param custom = false;
+
+    @chpldoc.nodoc
+    var agg: if aggregate then CustomDstAggregatorImpl(elemType, ?) else nothing;
+
+    proc init(type elemType) {
+      this.elemType = elemType;
+      if aggregate then this.agg = new CustomDstAggregatorImpl(elemType, 0);
+    }
+
+    proc init(ref handler) {
+      this.elemType = handler.elemType;
+      this.custom = true;
+      if aggregate then this.agg = new CustomDstAggregatorImpl(elemType, handler);
+    }
+
+    /*
+      Sets ``dst = srcVal`` in a way that aggregates such updates
+      to improve communication efficiency assuming that ``dst`` is remote
+      and ``srcVal`` is local.
+    */
+    inline proc ref copy(ref dst: elemType, const in srcVal: elemType) {
+      if aggregate then agg.copy(dst, srcVal);
+                   else dst = srcVal;
+    }
+
+    /*
+      Copy method for custom aggregator.
+    */
+    inline proc ref copy(const in srcVal: elemType) {
+      compilerAssert(aggregate);
+      agg.copy(srcVal);
+    }
+
+    /*
+      Flushes the aggregator & completes the updates queued up from the
+      :proc:`DstAggregator.copy` calls.
+
+      :arg freeBuffers: if ``true``, deallocates buffers used by this
+                        aggregator. If ``false``, the buffers will remain
+                        allocated after this ``flush`` (to support further
+                        :proc:`DstAggregator.copy` calls) and deallocated
+                        when the aggregator variable is deinitialized.
+    */
+    inline proc ref flush(freeBuffers=false) {
+      if aggregate then agg.flush(freeBuffers=freeBuffers);
+    }
+  }
+
+  @chpldoc.nodoc
+  record CustomDstAggregatorImpl {
+    type elemType;
+    var handler;
+    type aggType = (c_ptr(elemType), elemType);
+    const bufferSize = dstBuffSize;
+    const myLocaleSpace = 0..<numLocales;
+    var lastLocale: int;
+    var opsUntilYield = yieldFrequency;
+    var lBuffers: c_ptr(c_ptr(aggType));
+    var rBuffers: [myLocaleSpace] remoteBuffer(aggType);
+    var rHandlers: [myLocaleSpace] remoteHandler(handler.type);
+    var bufferIdxs: c_ptr(int);
+
+    proc ref postinit() {
+      lBuffers = allocate(c_ptr(aggType), numLocales);
+      bufferIdxs = bufferIdxAlloc();
+      for loc in myLocaleSpace {
+        lBuffers[loc] = allocate(aggType, bufferSize);
+        bufferIdxs[loc] = 0;
+        rBuffers[loc] = new remoteBuffer(aggType, bufferSize, loc);
+        rHandlers[loc] = new remoteHandler(handler.type, handler, loc);
+      }
+    }
+
+    proc ref deinit() {
+      flush(freeBuffers=true);
+      for loc in myLocaleSpace {
+        deallocate(lBuffers[loc]);
+      }
+      deallocate(lBuffers);
+      deallocate(bufferIdxs);
+    }
+
+    proc ref flush(freeBuffers: bool) {
+      for offsetLoc in myLocaleSpace + lastLocale {
+        const loc = offsetLoc % numLocales;
+        _flushBuffer(loc, bufferIdxs[loc], freeData=freeBuffers);
+      }
+    }
+
+    inline proc ref copy(ref dst: elemType, const in srcVal: elemType) {
+      if verboseAggregation {
+        writeln("DstAggregator.copy is called");
+      }
+      // Get the locale of dst and the local address on that locale
+      const loc = dst.locale.id;
+      lastLocale = loc;
+      const dstAddr = getAddr(dst);
+
+      // Get our current index into the buffer for dst's locale
+      ref bufferIdx = bufferIdxs[loc];
+
+      // Buffer the address and desired value
+      lBuffers[loc][bufferIdx] = (dstAddr, srcVal);
+      bufferIdx += 1;
+
+      // Flush our buffer if it's full. If it's been a while since we've let
+      // other tasks run, yield so that we're not blocking remote tasks from
+      // flushing their buffers.
+      if bufferIdx == bufferSize {
+        _flushBuffer(loc, bufferIdx, freeData=false);
+        opsUntilYield = yieldFrequency;
+      } else if opsUntilYield == 0 {
+        currentTask.yieldExecution();
+        opsUntilYield = yieldFrequency;
+      } else {
+        opsUntilYield -= 1;
+      }
+    }
+
+    inline proc ref copy(const in srcVal: elemType) {
+      if verboseAggregation {
+        writeln("DstAggregator.copy is called");
+      }
+      // Get the locale of srcVal and set dstAddr to nil
+      const loc = handler.getDestinationLocale(srcVal).id;
+      lastLocale = loc;
+      const dstAddr = nil;
+
+      // Get our current index into the buffer for dst's locale
+      ref bufferIdx = bufferIdxs[loc];
+
+      // Buffer the address and desired value
+      lBuffers[loc][bufferIdx] = (dstAddr,srcVal);
+      bufferIdx += 1;
+
+      // Flush our buffer if it's full. If it's been a while since we've let
+      // other tasks run, yield so that we're not blocking remote tasks from
+      // flushing their buffers.
+      if bufferIdx == bufferSize {
+        _flushBuffer(loc, bufferIdx, freeData=false);
+        opsUntilYield = yieldFrequency;
+      } else if opsUntilYield == 0 {
+        currentTask.yieldExecution();
+        opsUntilYield = yieldFrequency;
+      } else {
+        opsUntilYield -= 1;
+      }
+    }
+
+    proc ref _flushBuffer(loc: int, ref bufferIdx, freeData) {
+      const myBufferIdx = bufferIdx;
+      if myBufferIdx == 0 then return;
+
+      // Allocate a remote buffer
+      ref rBuffer = rBuffers[loc];
+      const remBufferPtr = rBuffer.cachedAlloc();
+      ref rHandler = rHandlers[loc].get();
+
+      // Copy local buffer to remote buffer
+      rBuffer.PUT(lBuffers[loc], myBufferIdx);
+
+      // Process remote buffer
+      on Locales[loc] {
+        if !isNothing(handler) then rHandler!.flush(rBuffer, remBufferPtr, myBufferIdx);
+      }
+      if freeData {
+        rBuffer.markFreed();
+      }
+      bufferIdx = 0;
+    }
+  }
+}
diff --git a/src/MultiTypeSymEntry.chpl b/src/MultiTypeSymEntry.chpl
index 5710da91d17..10499272988 100644
--- a/src/MultiTypeSymEntry.chpl
+++ b/src/MultiTypeSymEntry.chpl
@@ -97,8 +97,8 @@ module MultiTypeSymEntry
         }
 
         /**
-         * Formats and returns data in this entry up to the specified threshold. 
-         * Arrays of size less than threshold will be printed in their entirety. 
+         * Formats and returns data in this entry up to the specified threshold.
+         * Arrays of size less than threshold will be printed in their entirety.
          * Arrays of size greater than or equal to threshold will print the first 3 and last 3 elements
          *
          * :arg thresh: threshold for data to return
@@ -123,7 +123,7 @@ module MultiTypeSymEntry
     }
 
     /* Casts a GenSymEntry to the specified type and returns it.
-       
+
        :arg gse: generic sym entry
        :type gse: borrowed GenSymEntry
 
@@ -134,11 +134,11 @@ module MultiTypeSymEntry
         return gse.toSymEntry(etype, dimensions);
     }
 
-    /* 
+    /*
         This is a dummy class to avoid having to talk about specific
-        instantiations of SymEntry. 
+        instantiations of SymEntry.
         GenSymEntries can contain multiple SymEntries, but they represent a singular object.
-        For example, SegArray contains the offsets and values array, but only the values are 
+        For example, SegArray contains the offsets and values array, but only the values are
         considered data.
     */
     class GenSymEntry:AbstractSymEntry
@@ -189,9 +189,9 @@ module MultiTypeSymEntry
             }
         }
 
-        /* 
-          Formats and returns data in this entry up to the specified threshold. 
-          Arrays of size less than threshold will be printed in their entirety. 
+        /*
+          Formats and returns data in this entry up to the specified threshold.
+          Arrays of size less than threshold will be printed in their entirety.
           Arrays of size greater than or equal to threshold will print the first 3 and last 3 elements
 
             :arg thresh: threshold for data to return
@@ -319,8 +319,8 @@ module MultiTypeSymEntry
         }
 
         /*
-        Formats and returns data in this entry up to the specified threshold. 
-        Arrays of size less than threshold will be printed in their entirety. 
+        Formats and returns data in this entry up to the specified threshold.
+        Arrays of size less than threshold will be printed in their entirety.
         Arrays of size greater than or equal to threshold will print the first 3 and last 3 elements
 
             :arg thresh: threshold for data to return
@@ -456,7 +456,7 @@ module MultiTypeSymEntry
      * Factory method for creating a typed SymEntry and checking mem limits
      * :arg len: the number of elements to allocate
      * :type len: int
-     * 
+     *
      * :arg t: the element type
      * :type t: type
     */
@@ -498,8 +498,8 @@ module MultiTypeSymEntry
         }
 
         /**
-         * Formats and returns data in this entry up to the specified threshold. 
-         * Arrays of size less than threshold will be printed in their entirety. 
+         * Formats and returns data in this entry up to the specified threshold.
+         * Arrays of size less than threshold will be printed in their entirety.
          * Arrays of size greater than or equal to threshold will print the first 3 and last 3 elements
          *
          * :arg thresh: threshold for data to return
@@ -738,7 +738,6 @@ module MultiTypeSymEntry
         }
     }
 
-
     class GeneratorSymEntry:AbstractSymEntry {
         type etype;
         var generator: randomStream(etype);
diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl
index a1fa325908b..77367085b57 100644
--- a/src/SparseMatrix.chpl
+++ b/src/SparseMatrix.chpl
@@ -3,7 +3,13 @@ module SparseMatrix {
   public use SpsMatUtil;
   use ArkoudaSparseMatrixCompat;
   use BlockDist;
+  use CompressedSparseLayout;
   use CommAggregation;
+  use CustomCopyAggregation;
+  use PrivateDist;
+  use ChapelLocks;
+
+  config param aggregatedSparseMatrixCreation = true;
 
   // Quick and dirty, not permanent
   proc fillSparseMatrix(ref spsMat, const A: [?D] ?eltType, param l: Layout) throws {
@@ -411,10 +417,12 @@ module SparseMatrix {
             sparseMatMatMult(aBlk, bBlk, spsData);
           }
         }
-
         // Get my locale's local indices and create a sparse matrix
         // using them and the spsData computed above.
         //
+        // TODO: the 'make*SparseMat' functions create an array of size 'NNZ'.
+        // If we need to improve memory usage, could we use a SparseIndexBuffer
+        // instead?
         const locInds = A.domain.parentDom.localSubdomain();
         var cBlk = makeSparseMat(locInds, spsData);
 
@@ -424,6 +432,7 @@ module SparseMatrix {
         C.setLocalSubarray(cBlk);
       }
     }
+
     return C;
   }
 
@@ -470,7 +479,6 @@ module SparseMatrix {
   proc sparseMatFromArrays(rows, cols, vals, shape: 2*int, param layout,
                            type eltType) throws {
     import SymArrayDmap.makeSparseDomain;
-    var (SD, dense) = makeSparseDomain(shape, layout);
 
     const minRow = min reduce rows;
     const maxRow = max reduce rows;
@@ -497,9 +505,9 @@ module SparseMatrix {
         errorClass="InvalidArgumentError"
         );
 
+    var (SD, dense) = makeSparseDomain(shape, layout);
     var A: [SD] eltType;
     addElementsToSparseArray(A, SD, rows, cols, vals);
-
     return A;
   }
 
@@ -520,36 +528,124 @@ module SparseMatrix {
     }
   }
 
+  //
+  // We need a lock per locale to ensure that the flush() method of
+  // DestinationHandler is thread-safe. Wrapping a PrivateDist array in a
+  // class is the most expedient way to achieve this.
+  //
+  // This may result in an extra GET or two when accessing the 'data' field.
+  //
+  class LockHelper {
+    var data : [PrivateSpace] chpl_LocalSpinlock;
+  }
+
+  use ChplConfig;
+  config param bufSize = 1024;
+
+  class DestinationHandler {
+    var domVal;
+    var arrVal;
+    var lockObj;
+
+    proc init(domVal, arrVal, lockObj) {
+      this.domVal = domVal;
+      this.arrVal = arrVal;
+      this.lockObj = lockObj;
+    }
+
+    inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) {
+      ref lock = lockObj.data[here.id];
+      lock.lock();
+      defer lock.unlock();
+
+      const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here);
+      var locDomVal = this.domVal.locDoms[locid]!.mySparseBlock._value;
+      var locIdxBuf = locDomVal.dsiCreateIndexBuffer(bufSize,false,false);
+      for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) {
+        assert(dstAddr == nil);
+        var (i,j,_) = srcVal;
+        locIdxBuf.add((i, j));
+      }
+      locIdxBuf.commit();
+
+      // We use a lock around the entire 'flush' call because another parallel
+      // 'flush' might be inserting indices, which will make the returned
+      // results from 'find' invalid.
+      //
+      // TODO:
+      // - try 'forall' loop here
+      // - can we create a version of bulkAdd that handles the values as well?
+      for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) {
+        assert(dstAddr == nil);
+        var (i,j,v) = srcVal;
+        var (_,loc) = locDomVal.find((i,j));
+        this.arrVal.locArr[locid]!.myElems._value.data[loc] = v;
+      }
+    }
+  }
+
+  class SourceHandler {
+    var domVal;
+    var arrVal;
+    var lockObj;
+    type elemType = (int,int,int);
+
+    proc init(D, A, locks) {
+      this.domVal = D._value;
+      this.arrVal = A._value;
+      this.lockObj = locks;
+    }
+
+    proc sourceCopy() {
+      return new unmanaged DestinationHandler(domVal,arrVal, lockObj);
+    }
+
+    proc getDestinationLocale(val: elemType) {
+      // Since elemType is a tuple of (i,j,v) then we only need (i,j)
+      var (i,j,_) = val;
+      return domVal.dist.dsiIndexToLocale((i,j));
+    }
+  }
 
   proc addElementsToSparseArray(ref A, ref SD, const ref rows, const ref cols,
                                 const ref vals) throws where
                                 !A.chpl_isNonDistributedArray() {
-    coforall (loc, locDom) in zip(getGrid(A),
-                                  SD._value.locDoms) {
-      on loc {
-        for _srcLocId in loc.id..#numLocales {
-          const srcLocId = _srcLocId % numLocales;
-          var rowChunk = rows[rows.localSubdomain(Locales[srcLocId])];
-          var colChunk = cols[rows.localSubdomain(Locales[srcLocId])];
-          var valChunk = vals[rows.localSubdomain(Locales[srcLocId])];
-          for (r,c,v) in zip(rowChunk, colChunk, valChunk) {
-            if locDom!.parentDom.contains(r,c) {
-              if locDom!.mySparseBlock.contains(r,c) then
-                throw getErrorWithContext(
-                  msg="Duplicate index (%i, %i) in sparse matrix".format(r, c),
-                  lineNumber=getLineNumber(),
-                  routineName=getRoutineName(),
-                  moduleName=getModuleName(),
-                  errorClass="InvalidArgumentError"
-                  );
-
-
-              locDom!.mySparseBlock += (r,c);
-              A[r,c] = v;
+    
+    if !aggregatedSparseMatrixCreation {
+      coforall (loc, locDom) in zip(getGrid(A),
+                                    SD._value.locDoms) {
+        on loc {
+          for _srcLocId in loc.id..#numLocales {
+            const srcLocId = _srcLocId % numLocales;
+            var rowChunk = rows[rows.localSubdomain(Locales[srcLocId])];
+            var colChunk = cols[rows.localSubdomain(Locales[srcLocId])];
+            var valChunk = vals[rows.localSubdomain(Locales[srcLocId])];
+            for (r,c,v) in zip(rowChunk, colChunk, valChunk) {
+              if locDom!.parentDom.contains(r,c) {
+                if locDom!.mySparseBlock.contains(r,c) then
+                  throw getErrorWithContext(
+                    msg="Duplicate index (%i, %i) in sparse matrix".format(r, c),
+                    lineNumber=getLineNumber(),
+                    routineName=getRoutineName(),
+                    moduleName=getModuleName(),
+                    errorClass="InvalidArgumentError"
+                    );
+
+
+                locDom!.mySparseBlock += (r,c);
+                A[r,c] = v;
+              }
             }
           }
         }
       }
+    } else {
+      var locks = new unmanaged LockHelper();
+      defer delete locks;
+
+      forall (i,j,v) in zip(rows, cols, vals)
+        with (var agg = new CustomDstAggregator(new shared SourceHandler(SD, A, locks))) do
+          agg.copy((i,j,v));
     }
 
   }
@@ -624,17 +720,16 @@ module SparseMatrix {
 
       sort(inds);
 
-      for ij in inds do
-        CDom += ij;
+      CDom.bulkAdd(inds, true, true);
 
       var C: [CDom] int;
+      // TODO: can this be parallel?
       for ij in inds do
-        try! C[ij] += spsData[ij];  // TODO: Should this really throw?
+        try! C[ij] = spsData[ij];  // TODO: Should this really throw?
 
       return C;
     }
 
-
     // create a new sparse matrix from a collection of nonzero indices
     // (nnzs) and values (vals)
     //
diff --git a/src/SparseMatrixMsg.chpl b/src/SparseMatrixMsg.chpl
index 587e05b666a..0e50627b75f 100644
--- a/src/SparseMatrixMsg.chpl
+++ b/src/SparseMatrixMsg.chpl
@@ -46,7 +46,8 @@ module SparseMatrixMsg {
     proc sparseMatrixtoPdarray(cmd: string, msgArgs: borrowed MessageArgs, st: borrowed SymTab,
                                type SparseSymEntry_etype, param SparseSymEntry_matLayout: Layout
     ): MsgTuple throws {
-        const e = st[msgArgs["matrix"]]: borrowed SparseSymEntry(SparseSymEntry_etype, 2, SparseSymEntry_matLayout);
+        type castTo = borrowed SparseSymEntry(SparseSymEntry_etype, 2, SparseSymEntry_matLayout);
+        const e = st[msgArgs["matrix"]]: castTo;
 
         const size = e.nnz;
         var rows = makeDistArray(size, int),
@@ -86,6 +87,7 @@ module SparseMatrixMsg {
               shape = msgArgs["shape"].toScalarTuple(int, 2); // Hardcode 2D for now
 
         const aV = sparseMatFromArrays(rows.a, cols.a, vals.a, shape, SparseSymEntry_matLayout, SparseSymEntry_etype);
+
         return st.insert(new shared SparseSymEntry(aV, SparseSymEntry_matLayout));
     }
 
diff --git a/tests/scipy/sparse_test.py b/tests/scipy/sparse_test.py
index 799edea8ed1..0e2ae209d2a 100644
--- a/tests/scipy/sparse_test.py
+++ b/tests/scipy/sparse_test.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import scipy.sparse as sp
 
 import arkouda as ak
 
@@ -206,3 +207,44 @@ def test_creation_csr(self):
         assert np.all(vals == vals_)
         # Check the layout is correct
         assert mat.layout == layout
+
+    def test_to_scipy_sparse(self):
+        rows = ak.array([1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 8])
+        cols = ak.array([4, 5, 6, 1, 3, 4, 2, 3, 1, 2, 8, 1, 2, 6, 1, 6, 7, 8])
+        vals = ak.array(
+            [
+                3,
+                20,
+                30,
+                10,
+                40,
+                50,
+                60,
+                70,
+                80,
+                90,
+                100,
+                110,
+                120,
+                130,
+                140,
+                150,
+                160,
+                170,
+            ]
+        )
+
+        mat = create_sparse_matrix(10, rows, cols, vals, "CSR")
+        conv = mat.to_scipy_sparse()
+        info = (vals.to_ndarray(), (rows.to_ndarray(), cols.to_ndarray()))
+        expected = sp.csr_array(info, shape=(10, 10))
+
+        assert sp.issparse(conv)
+        assert conv.shape == expected.shape
+        assert conv.format == expected.format
+
+        conv_coo = conv.tocoo()
+        expected_coo = expected.tocoo()
+        assert np.array_equal(conv_coo.row, expected_coo.row)
+        assert np.array_equal(conv_coo.col, expected_coo.col)
+        assert np.array_equal(conv_coo.data, expected_coo.data)