From e42c4d48d1ed750cbf9a5ddbaafd608f4c574813 Mon Sep 17 00:00:00 2001 From: Oliver Alvarado Rodriguez Date: Fri, 12 Sep 2025 16:16:45 -0500 Subject: [PATCH 01/17] initial agg sparse matrix creation Signed-off-by: Oliver Alvarado Rodriguez --- src/CustomCopyAggregation.chpl | 237 +++++++++++++++++++++++++++++++++ src/MultiTypeSymEntry.chpl | 175 ++++++++++++++++++++++++ src/SparseMatrix.chpl | 147 ++++++++++++++++---- src/SparseMatrixMsg.chpl | 26 +++- 4 files changed, 551 insertions(+), 34 deletions(-) create mode 100644 src/CustomCopyAggregation.chpl diff --git a/src/CustomCopyAggregation.chpl b/src/CustomCopyAggregation.chpl new file mode 100644 index 00000000000..34e2c147d51 --- /dev/null +++ b/src/CustomCopyAggregation.chpl @@ -0,0 +1,237 @@ +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..= thresh { + var count = 0; + // Normal iteration like this is more efficient than + // dense iteration, so we prefer that for the first elements + for (_, (i, j)) in zip(1..3, sparseDom) { + const idxStr = " (%?, %?)".format(i, j); // Padding to match SciPy + s += "%<16s%?\n".format(idxStr, this.a[i,j]); + } + + s += " : :\n"; // Dot dot seperator, but vertical + + // For the last elements, we iterate in dense order + // Since sparseArrays cant be strided by -1 + // We also have to do some i,j swaps for CSC vs CSR differences + count = 0; + var backString = ""; + for (i, j) in denseDom by -1 { + var row = i, col = j; + if this.matLayout==Layout.CSC { + row = j; // Iterate in Col Major Order for CSC + col = i; // To match SciPy behavior + } + if !sparseDom.contains(row, col) then continue; + const idxStr = " (%?, %?)".format(row, col); // Padding to match SciPy + backString = "%<16s%?\n".format(idxStr, this.a[row,col]) + backString; + count += 1; + if count == 3 then break; + } + s+=backString; + } else { + for (i,j) in sparseDom { + const idxStr = " (%?, %?)".format(i, j); // Padding to match SciPy + s += "%<16s%?\n".format(idxStr, this.a[i,j]); + } + } + + if this.etype == bool { + s = s.replace("true","True"); + s = s.replace("false","False"); + } + + return s; + } + + /* + Verbose flag utility method + */ + proc deinit() { + if logLevel == LogLevel.DEBUG {writeln("deinit SparseSymEntry");try! stdout.flush();} + } + } + class GeneratorSymEntry:AbstractSymEntry { type etype; diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index a1fa325908b..6d0c3f3f44f 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -3,7 +3,11 @@ module SparseMatrix { public use SpsMatUtil; use ArkoudaSparseMatrixCompat; use BlockDist; + use CompressedSparseLayout; use CommAggregation; + use CustomCopyAggregation; + + config param aggregatedSparseMatrixCreation = true; // Quick and dirty, not permanent proc fillSparseMatrix(ref spsMat, const A: [?D] ?eltType, param l: Layout) throws { @@ -416,7 +420,8 @@ module SparseMatrix { // using them and the spsData computed above. // const locInds = A.domain.parentDom.localSubdomain(); - var cBlk = makeSparseMat(locInds, spsData); + var cBlk = if !aggregatedSparseMatrixCreation then makeSparseMat(locInds, spsData) + else makeParSafeSparseMat(locInds, spsData); // Stitch the local portions back together into the global-view // @@ -470,7 +475,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,10 +501,17 @@ module SparseMatrix { errorClass="InvalidArgumentError" ); - var A: [SD] eltType; - addElementsToSparseArray(A, SD, rows, cols, vals); - - return A; + if !aggregatedSparseMatrixCreation { + var (SD, dense) = makeSparseDomain(shape, layout); + var A: [SD] eltType; + addElementsToSparseArray(A, SD, rows, cols, vals); + return A; + } else { + var (SD, dense) = makeParSafeSparseDomain(shape, layout); + var A: [SD] eltType; + addElementsToSparseArray(A, SD, rows, cols, vals); + return A; + } } proc addElementsToSparseArray(ref A, ref SD, const ref rows,const ref cols, @@ -520,36 +531,93 @@ module SparseMatrix { } } + use ChplConfig; + config param bufSize = 1024; + + class DestinationHandler { + var domVal; + var arrVal; + + proc init(domVal, arrVal) { + this.domVal = domVal; + this.arrVal = arrVal; + } + + inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) { + const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here); + var locIdxBuf = this.domVal.locDoms[locid]!.mySparseBlock._value.createIndexBuffer(bufSize,true,true); + for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) { + assert(dstAddr == nil); + var (i,j,_) = srcVal; + locIdxBuf.add((i, j)); + } + locIdxBuf.commit(); + for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) { + assert(dstAddr == nil); + var (i,j,v) = srcVal; + var (_,loc) = this.domVal.locDoms[locid]!.mySparseBlock._value.find((i,j)); + this.arrVal.locArr[locid]!.myElems._value.data[loc] = v; + } + } + } + + class SourceHandler { + var domVal; + var arrVal; + type elemType = (int,int,int); + + proc init(D, A) { + this.domVal = D._value; + this.arrVal = A._value; + } + + proc sourceCopy() { + return new unmanaged DestinationHandler(domVal,arrVal); + } + + 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 { + forall (i,j,v) in zip(rows, cols, vals) + with (var agg = new CustomDstAggregator(new shared SourceHandler(SD, A))) do + agg.copy((i,j,v)); } } @@ -634,6 +702,29 @@ module SparseMatrix { return C; } + // create a new sparse matrix from a map from sparse indices to values + // + proc makeParSafeSparseMat(parentDom, spsData) { + use ArkoudaSparseMatrixCompat; + use Sort; + + var CDom: sparse subdomain(parentDom) dmapped getParSafeSparseDom(Layout.CSR); + var inds: [0.. Date: Thu, 2 Oct 2025 09:02:46 -0500 Subject: [PATCH 02/17] `isSorted` and `isUnique` turned off for `createIndexBuffer` Signed-off-by: Oliver Alvarado Rodriguez --- src/SparseMatrix.chpl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index 6d0c3f3f44f..c2f8788dc32 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -545,7 +545,7 @@ module SparseMatrix { inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) { const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here); - var locIdxBuf = this.domVal.locDoms[locid]!.mySparseBlock._value.createIndexBuffer(bufSize,true,true); + var locIdxBuf = this.domVal.locDoms[locid]!.mySparseBlock._value.createIndexBuffer(bufSize,false,false); for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) { assert(dstAddr == nil); var (i,j,_) = srcVal; From 3fd9039529075ef8a7e95f03b08774b7197dbc20 Mon Sep 17 00:00:00 2001 From: Oliver Alvarado Rodriguez Date: Tue, 21 Oct 2025 14:55:24 -0500 Subject: [PATCH 03/17] added rmat creator Signed-off-by: Oliver Alvarado Rodriguez --- arkouda/scipy/sparsematrix.py | 54 ++++++++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/arkouda/scipy/sparsematrix.py b/arkouda/scipy/sparsematrix.py index fc6f6b38196..7c9bdb06610 100644 --- a/arkouda/scipy/sparsematrix.py +++ b/arkouda/scipy/sparsematrix.py @@ -11,13 +11,65 @@ from arkouda.numpy.dtypes import NumericDTypes, int64 from arkouda.numpy.dtypes import dtype as akdtype from arkouda.numpy.pdarrayclass import pdarray +from arkouda.numpy.pdarraysetops import concatenate from arkouda.scipy.sparrayclass import create_sparray, sparray +from arkouda.numpy.pdarraycreation import zeros +from arkouda.random import randint +from arkouda.sorting import argsort +from arkouda.pandas.groupbyclass import unique -__all__ = ["random_sparse_matrix", "sparse_matrix_matrix_mult", "create_sparse_matrix"] +__all__ = [ + "rmat", + "random_sparse_matrix", + "sparse_matrix_matrix_mult", + "create_sparse_matrix" + ] logger = get_arkouda_logger(name="sparsematrix") +@typechecked +def rmat( + scale: int, + a = 0.57, b = 0.19, c = 0.19, d = 0.05, + edge_factor = 16 +) -> tuple[pdarray,pdarray]: + p = (a,b,c,d) + n = 2 ** scale # number vertices + m = n * edge_factor # number edges + + if isinstance(p, float) and 0 <= p <= 1: + a = p + b = c = d = (1.0 - p) / 3.0 + elif isinstance(p, tuple) and all(0 <= x <= 1 for x in p) and sum(p) == 1: + a, b, c, d = p + else: + raise ValueError(f"p = {p} doesn't represent valid RMAT probability.") + + ab, cNorm, aNorm = a + b, c / (c + d), a / (a + b) + + U, V = zeros(m, dtype="int64"), zeros(m, dtype="int64") + + for s in range(0, scale): + uMask = randint(0, 1, m, dtype="float64") > ab + vMask = (randint(0, 1, m, dtype="float64") > + (cNorm * uMask + aNorm * (~uMask)) + ) + U += uMask * (2 ** s) + V += vMask * (2 ** s) + + pi = argsort(randint(0, 1, n, dtype="float64")) + U, V = pi[U], pi[V] + + Us = concatenate([U,V]) + Vs = concatenate([V,U]) + + U,V = unique([Us, Vs]) + U += 1 + V += 1 + + return (U,V) + @typechecked def random_sparse_matrix( From 6217db1ce99b81c478d142aedebb8597849fc3df Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Wed, 1 Jul 2026 14:27:54 -0700 Subject: [PATCH 04/17] Update to avoid deprecated imports --- arkouda/scipy/sparsematrix.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/arkouda/scipy/sparsematrix.py b/arkouda/scipy/sparsematrix.py index 7c9bdb06610..3cc5e7bc7d8 100644 --- a/arkouda/scipy/sparsematrix.py +++ b/arkouda/scipy/sparsematrix.py @@ -14,8 +14,8 @@ from arkouda.numpy.pdarraysetops import concatenate from arkouda.scipy.sparrayclass import create_sparray, sparray from arkouda.numpy.pdarraycreation import zeros -from arkouda.random import randint -from arkouda.sorting import argsort +from arkouda.numpy.random import randint +from arkouda.numpy.sorting import argsort from arkouda.pandas.groupbyclass import unique From 8cc0d6b65fd2b15fa4a3b2432aea4396c692b1c3 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Fri, 31 Jul 2026 09:42:27 -0700 Subject: [PATCH 05/17] Fix issues with converting sparse array to pdarray --- arkouda/scipy/sparrayclass.py | 2 +- src/SparseMatrixMsg.chpl | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/arkouda/scipy/sparrayclass.py b/arkouda/scipy/sparrayclass.py index 45df14d053c..1e57833f15e 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 diff --git a/src/SparseMatrixMsg.chpl b/src/SparseMatrixMsg.chpl index 98103ee0045..486f94b0c23 100644 --- a/src/SparseMatrixMsg.chpl +++ b/src/SparseMatrixMsg.chpl @@ -56,7 +56,11 @@ 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 = if aggregatedSparseMatrixCreation then + borrowed ParSafeSparseSymEntry(SparseSymEntry_etype, 2, SparseSymEntry_matLayout) + else + borrowed SparseSymEntry(SparseSymEntry_etype, 2, SparseSymEntry_matLayout); + const e = st[msgArgs["matrix"]]: castTo; const size = e.nnz; var rows = makeDistArray(size, int), From 516bedaff0163e7a5f8430a155680c73ea822bc7 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Fri, 31 Jul 2026 09:42:58 -0700 Subject: [PATCH 06/17] Add 'to_scipy_sparse' to sparrayclass --- arkouda/scipy/sparrayclass.py | 35 ++++++++++++++++++++++++++++++ tests/scipy/sparse_test.py | 41 +++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) diff --git a/arkouda/scipy/sparrayclass.py b/arkouda/scipy/sparrayclass.py index 1e57833f15e..892fdb34406 100644 --- a/arkouda/scipy/sparrayclass.py +++ b/arkouda/scipy/sparrayclass.py @@ -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/tests/scipy/sparse_test.py b/tests/scipy/sparse_test.py index 799edea8ed1..3227356b1ac 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,43 @@ 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() + expected = sp.csr_array((vals.to_ndarray(), (rows.to_ndarray(), cols.to_ndarray())), 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) From dbbce857e3f6cf710b5001cbf249a03f5f5fa0d8 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Sat, 1 Aug 2026 16:49:27 -0700 Subject: [PATCH 07/17] Bulk-add indices since they're already sorted Also, no need to use += when the summation already occurred. Fixes a bug with the custom copy aggregation code where we need to be locking the domain while we're inserting values. Otherwise, 'find' will return bogus indices to insert the data. Finally, perform some cleanup of the original effort. --- arkouda/scipy/sparsematrix.py | 53 +------------ src/CustomCopyAggregation.chpl | 28 +++---- src/MultiTypeSymEntry.chpl | 75 ++++--------------- src/SparseMatrix.chpl | 31 ++++++-- src/SymArrayDmap.chpl | 26 +++++++ .../ge-24/ArkoudaSparseMatrixCompat.chpl | 19 +++++ 6 files changed, 94 insertions(+), 138 deletions(-) diff --git a/arkouda/scipy/sparsematrix.py b/arkouda/scipy/sparsematrix.py index 3cc5e7bc7d8..545eb587b90 100644 --- a/arkouda/scipy/sparsematrix.py +++ b/arkouda/scipy/sparsematrix.py @@ -11,66 +11,17 @@ from arkouda.numpy.dtypes import NumericDTypes, int64 from arkouda.numpy.dtypes import dtype as akdtype from arkouda.numpy.pdarrayclass import pdarray -from arkouda.numpy.pdarraysetops import concatenate from arkouda.scipy.sparrayclass import create_sparray, sparray -from arkouda.numpy.pdarraycreation import zeros -from arkouda.numpy.random import randint -from arkouda.numpy.sorting import argsort -from arkouda.pandas.groupbyclass import unique __all__ = [ - "rmat", - "random_sparse_matrix", - "sparse_matrix_matrix_mult", + "random_sparse_matrix", + "sparse_matrix_matrix_mult", "create_sparse_matrix" ] logger = get_arkouda_logger(name="sparsematrix") -@typechecked -def rmat( - scale: int, - a = 0.57, b = 0.19, c = 0.19, d = 0.05, - edge_factor = 16 -) -> tuple[pdarray,pdarray]: - p = (a,b,c,d) - n = 2 ** scale # number vertices - m = n * edge_factor # number edges - - if isinstance(p, float) and 0 <= p <= 1: - a = p - b = c = d = (1.0 - p) / 3.0 - elif isinstance(p, tuple) and all(0 <= x <= 1 for x in p) and sum(p) == 1: - a, b, c, d = p - else: - raise ValueError(f"p = {p} doesn't represent valid RMAT probability.") - - ab, cNorm, aNorm = a + b, c / (c + d), a / (a + b) - - U, V = zeros(m, dtype="int64"), zeros(m, dtype="int64") - - for s in range(0, scale): - uMask = randint(0, 1, m, dtype="float64") > ab - vMask = (randint(0, 1, m, dtype="float64") > - (cNorm * uMask + aNorm * (~uMask)) - ) - U += uMask * (2 ** s) - V += vMask * (2 ** s) - - pi = argsort(randint(0, 1, n, dtype="float64")) - U, V = pi[U], pi[V] - - Us = concatenate([U,V]) - Vs = concatenate([V,U]) - - U,V = unique([Us, Vs]) - U += 1 - V += 1 - - return (U,V) - - @typechecked def random_sparse_matrix( size: int, density: float, layout: str, dtype: Union[type, str] = int64 diff --git a/src/CustomCopyAggregation.chpl b/src/CustomCopyAggregation.chpl index 34e2c147d51..81985c7d07f 100644 --- a/src/CustomCopyAggregation.chpl +++ b/src/CustomCopyAggregation.chpl @@ -47,10 +47,10 @@ module CustomCopyAggregation { } } - /* + /* 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. + High memory usage since there are per-destination buffers. */ record CustomDstAggregator { type elemType; @@ -69,26 +69,26 @@ module CustomCopyAggregation { 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. + 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. + 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. @@ -219,14 +219,6 @@ module CustomCopyAggregation { // Process remote buffer on Locales[loc] { if !isNothing(handler) then rHandler!.flush(rBuffer, remBufferPtr, myBufferIdx); - // for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) { - // if !isNothing(handler) { - // rHandler!.flush(dstAddr, srcVal); - // } - // else { - // dstAddr.deref() = srcVal; - // } - // } } if freeData { rBuffer.markFreed(); @@ -234,4 +226,4 @@ module CustomCopyAggregation { bufferIdx = 0; } } -} \ No newline at end of file +} diff --git a/src/MultiTypeSymEntry.chpl b/src/MultiTypeSymEntry.chpl index 0ea83a491de..3c0c29f58fd 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,53 +738,6 @@ module MultiTypeSymEntry } } - use CompressedSparseLayout; - - proc getParSafeSparseDom(param layout: Layout) { - select layout { - when Layout.CSR do return new csrLayout(parSafe=true); - when Layout.CSC do return new cscLayout(parSafe=true); - } - } - - proc getParSafeDenseDom(dom, localeGrid, param layout: Layout) { - if layout == Layout.CSR { - return dom dmapped new blockDist(boundingBox=dom, - targetLocales=localeGrid, - sparseLayoutType=csrLayout(parSafe=true)); - } else { - return dom dmapped new blockDist(boundingBox=dom, - targetLocales=localeGrid, - sparseLayoutType=cscLayout(parSafe=true)); - } - } - - proc makeParSafeSparseDomain(shape: 2*int, param matLayout: Layout) { - const dom = {1..shape[0], 1..shape[1]}; // TODO: change domain to be zero based? - select MyDmap { - when Dmap.defaultRectangular { - var spsDom: sparse subdomain(dom) dmapped getParSafeSparseDom(matLayout); - return (spsDom, dom); - } - when Dmap.blockDist { - const locsPerDim = sqrt(numLocales:real): int, - grid = {0.. Date: Wed, 19 Aug 2026 17:01:46 -0700 Subject: [PATCH 08/17] Add benchmark for sparse matrix performance --- benchmarks/graph_infra/arkouda.graph | 12 ++ benchmarks/graph_infra/sparse.perfkeys | 3 + benchmarks/run_benchmarks.py | 1 + benchmarks/sparse.py | 192 +++++++++++++++++++++++++ 4 files changed, 208 insertions(+) create mode 100644 benchmarks/graph_infra/sparse.perfkeys create mode 100644 benchmarks/sparse.py diff --git a/benchmarks/graph_infra/arkouda.graph b/benchmarks/graph_infra/arkouda.graph index 1415f6d6d30..64c64a55f80 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 (N=1M, NNZ=10M) +ylabel: Time (Seconds) + +perfkeys: Average CSCxCSR time = +graphkeys: CSCxCSR time +files: sparse.dat +graphtitle: Sparse Matrix Multiplication Time (N=1M, NNZ=10M) +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..cfdd078f2bc --- /dev/null +++ b/benchmarks/sparse.py @@ -0,0 +1,192 @@ +#!/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 = csr_matrix.dot(csc, 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") + csc = create_sparse_matrix(N, cols, rows, vals, layout="CSC") + result = sparse_matrix_matrix_mult(csc, 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" + + 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" + + scipy_result = csr_matrix.dot(scipy_csc, 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) + sys.exit(0) From bec6ab1c3fb0acd1c3a2efae42fbd1aa2e19acbc Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Mon, 24 Aug 2026 10:14:09 -0700 Subject: [PATCH 09/17] Use correct name for creation of index buffer --- src/SparseMatrix.chpl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index e542771b1c5..2c761713de6 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -549,7 +549,7 @@ module SparseMatrix { inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) { const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here); var locDomVal = this.domVal.locDoms[locid]!.mySparseBlock._value; - var locIdxBuf = locDomVal.createIndexBuffer(bufSize,false,false); + var locIdxBuf = locDomVal.dsiCreateIndexBuffer(bufSize,false,false); for (dstAddr, srcVal) in rBuffer.localIter(remBufferPtr, myBufferIdx) { assert(dstAddr == nil); var (i,j,_) = srcVal; From 1d18a180935150579f52c9f7a4156099ebe0f35a Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Mon, 24 Aug 2026 12:12:26 -0700 Subject: [PATCH 10/17] Resolve some CI failures --- .chplcheckignore | 1 + .flake8 | 2 +- arkouda/scipy/sparsematrix.py | 1 + tests/scipy/sparse_test.py | 3 ++- 4 files changed, 5 insertions(+), 2 deletions(-) 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/scipy/sparsematrix.py b/arkouda/scipy/sparsematrix.py index 545eb587b90..e20187d55cf 100644 --- a/arkouda/scipy/sparsematrix.py +++ b/arkouda/scipy/sparsematrix.py @@ -22,6 +22,7 @@ logger = get_arkouda_logger(name="sparsematrix") + @typechecked def random_sparse_matrix( size: int, density: float, layout: str, dtype: Union[type, str] = int64 diff --git a/tests/scipy/sparse_test.py b/tests/scipy/sparse_test.py index 3227356b1ac..0e2ae209d2a 100644 --- a/tests/scipy/sparse_test.py +++ b/tests/scipy/sparse_test.py @@ -236,7 +236,8 @@ def test_to_scipy_sparse(self): mat = create_sparse_matrix(10, rows, cols, vals, "CSR") conv = mat.to_scipy_sparse() - expected = sp.csr_array((vals.to_ndarray(), (rows.to_ndarray(), cols.to_ndarray())), shape=(10, 10)) + 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 From 32cf84c8cd6bb0a9ac38812a9d44df502fd362c9 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Mon, 24 Aug 2026 12:20:13 -0700 Subject: [PATCH 11/17] Add 'scipy-stubs' to developer dependencies --- arkouda-env-dev.yml | 1 + docker/ci/almalinux-with-arkouda-deps/requirements.txt | 1 + pydoc/requirements.txt | 1 + pydoc/setup/REQUIREMENTS.md | 1 + pyproject.toml | 1 + 5 files changed, 5 insertions(+) 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/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", From 3b05be1f0f59a40b7573951f895f0a628ab4a4d0 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Mon, 24 Aug 2026 12:25:20 -0700 Subject: [PATCH 12/17] Update formatting --- arkouda/scipy/sparsematrix.py | 6 +--- benchmarks/sparse.py | 65 +++++++++++++++++------------------ 2 files changed, 33 insertions(+), 38 deletions(-) diff --git a/arkouda/scipy/sparsematrix.py b/arkouda/scipy/sparsematrix.py index e20187d55cf..fc6f6b38196 100644 --- a/arkouda/scipy/sparsematrix.py +++ b/arkouda/scipy/sparsematrix.py @@ -14,11 +14,7 @@ from arkouda.scipy.sparrayclass import create_sparray, sparray -__all__ = [ - "random_sparse_matrix", - "sparse_matrix_matrix_mult", - "create_sparse_matrix" - ] +__all__ = ["random_sparse_matrix", "sparse_matrix_matrix_mult", "create_sparse_matrix"] logger = get_arkouda_logger(name="sparsematrix") diff --git a/benchmarks/sparse.py b/benchmarks/sparse.py index cfdd078f2bc..652347cf341 100644 --- a/benchmarks/sparse.py +++ b/benchmarks/sparse.py @@ -5,12 +5,11 @@ 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 -) +from arkouda.scipy.sparsematrix import create_sparse_matrix, sparse_matrix_matrix_mult # TODO: Add support for 'float64' -TYPES = ("int64", ) +TYPES = ("int64",) + def compare_scipy(left, right, rtol=1e-9, atol=0.0, equal_nan=False): if left.shape != right.shape: @@ -34,16 +33,16 @@ def time_ak_sparse(N, trials, dtype, seed): cfg = ak.get_config() print("numLocales = {}, numNodes {}, N = {:,}".format(cfg["numLocales"], cfg["numNodes"], N)) - nnz = N*10 + 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]) + 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) + vals = ak.randint(1, len(rows), len(rows), seed=seed * 3) elif dtype == "float64": - vals = ak.uniform(len(rows), seed=seed*3) + 0.5 + vals = ak.uniform(len(rows), seed=seed * 3) + 0.5 csr_times = [] csc_times = [] @@ -71,16 +70,16 @@ def time_np_sparse(N, trials, dtype, seed): print(">>> numpy {} sparse".format(dtype)) print("N = {:,}".format(N)) - nnz = N*10 + 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]) + 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) + vals = ak.randint(1, len(rows), len(rows), seed=seed * 3) elif dtype == "float64": - vals = ak.uniform(len(rows), seed=seed*3) + 0.5 + vals = ak.uniform(len(rows), seed=seed * 3) + 0.5 rows = rows.to_ndarray() cols = cols.to_ndarray() @@ -92,11 +91,11 @@ def time_np_sparse(N, trials, dtype, seed): for i in range(trials): start = time.time() - csr = coo_array((vals, (rows, cols)), shape=(N,N)).tocsr() + 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 = coo_array((vals, (cols, rows)), shape=(N, N)).tocsc() csc_times.append(time.time() - start) start = time.time() @@ -111,26 +110,30 @@ def time_np_sparse(N, trials, dtype, seed): def check_correctness(dtype): seed = 1234 N = 10_000 - nnz = N*10 + 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]) + 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) + vals = ak.randint(1, len(rows), len(rows), seed=seed * 3) elif dtype == "float64": - vals = ak.uniform(len(rows), seed=seed*3) + 0.5 + vals = ak.uniform(len(rows), seed=seed * 3) + 0.5 csr = create_sparse_matrix(N, rows, cols, vals, layout="CSR") csc = create_sparse_matrix(N, cols, rows, vals, layout="CSC") result = sparse_matrix_matrix_mult(csc, csr) - scipy_csr = coo_array((vals.to_ndarray(), (rows.to_ndarray(), cols.to_ndarray())), shape=(N, N)).tocsr() + 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" - scipy_csc = coo_array((vals.to_ndarray(), (cols.to_ndarray(), rows.to_ndarray())), shape=(N, N)).tocsc() + 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" @@ -140,16 +143,11 @@ def check_correctness(dtype): def create_parser(): - parser = argparse.ArgumentParser( - description="Benchmark sparse matrix creation and multiplication." - ) + 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("-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)) ) @@ -173,6 +171,7 @@ def create_parser(): if __name__ == "__main__": import sys + args = create_parser().parse_args() if args.dtype not in TYPES: From f8b4dae3923b3e74e8ed2f3d8e8627029e4b577c Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Tue, 25 Aug 2026 14:24:42 -0700 Subject: [PATCH 13/17] Use 'dot' method directly --- benchmarks/sparse.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/benchmarks/sparse.py b/benchmarks/sparse.py index 652347cf341..f503e1ca20e 100644 --- a/benchmarks/sparse.py +++ b/benchmarks/sparse.py @@ -99,7 +99,7 @@ def time_np_sparse(N, trials, dtype, seed): csc_times.append(time.time() - start) start = time.time() - result = csr_matrix.dot(csc, csr).tocsr() + result = csc.dot(csr).tocsr() multiplication_times.append(time.time() - start) print("Average CSR time = {:.4f} seconds".format(np.mean(csr_times))) @@ -137,7 +137,7 @@ def check_correctness(dtype): ak_csc = csc.to_scipy_sparse() assert compare_scipy(scipy_csc, ak_csc), "CSC matrices do not match" - scipy_result = csr_matrix.dot(scipy_csc, scipy_csr).tocsr() + 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" @@ -188,4 +188,8 @@ def create_parser(): 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) From d088de0b443d5c43389d76f706e63e2ad880c205 Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Thu, 27 Aug 2026 10:03:34 -0700 Subject: [PATCH 14/17] Use per-locale locks rather than lock embedded in sparse array --- src/MultiTypeSymEntry.chpl | 129 ------------------ src/SparseMatrix.chpl | 72 ++++------ src/SparseMatrixMsg.chpl | 30 +--- src/SymArrayDmap.chpl | 26 ---- .../ge-24/ArkoudaSparseMatrixCompat.chpl | 19 --- 5 files changed, 31 insertions(+), 245 deletions(-) diff --git a/src/MultiTypeSymEntry.chpl b/src/MultiTypeSymEntry.chpl index 3c0c29f58fd..10499272988 100644 --- a/src/MultiTypeSymEntry.chpl +++ b/src/MultiTypeSymEntry.chpl @@ -738,135 +738,6 @@ module MultiTypeSymEntry } } - /* Symbol table entry */ - class ParSafeSparseSymEntry : GenSparseSymEntry - { - /* - generic element type array - etype is different from dtype (chapel vs numpy) - */ - type etype; - - /* - number of dimensions, to be passed back to the `GenSparseSymEntry` so that - we are able to make it visible to the Python client - */ - param dimensions: int; // TODO: should we only support 2D sparse arrays and remove this field? - - /* - the actual shape of the array, this has to live here, since GenSparseSymEntry - has to stay generic - For now, each dimension is assumed to be equal. - */ - var tupShape: dimensions*int; - - /* - layout of the sparse array: CSC or CSR - */ - param matLayout : Layout; - - /* - 'a' is the distributed sparse array - */ - // Hardcode 2D matrix for now (makeSparseArray accepts a 2-tuple shape) - var a = makeParSafeSparseArray((...tupShape), etype, matLayout); - - /* - Create a SparseSymEntry from a sparse array - */ - proc init(a: [?D] ?eltType, param matLayout) - where a.domain.parentDom.rank == 2 // Hardcode a 2D matrix for now - { - const size = D.shape[0] * D.shape[1]; - super.init(eltType, size, a.domain.getNNZ(), /*ndim*/2, layoutToStr(matLayout)); // Hardcode a 2D matrix for now - this.entryType = SymbolEntryType.SparseSymEntry; - assignableTypes.add(this.entryType); - this.etype = eltType; - this.dimensions = 2; // Hardcode a 2D matrix for now - this.tupShape = D.shape; - this.matLayout = matLayout; - this.a = a; - init this; - this.shape = tupShapeString(this.tupShape); - this.ndim = 2; - } - - /* - Formats and returns data in this entry up to the specified threshold. - Matrices with nnz less than threshold will be printed in their entirety. - Matrices with nnz greater than or equal to threshold will print the first 3 and last 3 elements - - :arg thresh: threshold for data to return - :type thresh: int - - :arg prefix: String to pre-pend to the front of the data string - :type prefix: string - - :arg suffix: String to append to the tail of the data string - :type suffix: string - - :arg baseFormat: String which represents the base format string for the data type - :type baseFormat: string - - :returns: s (string) containing the array data - */ - override proc entry__str__(thresh:int=6, prefix:string = "noprefix", suffix:string = "nosuffix", baseFormat:string = "%?"): string throws { - var s:string; - const ref sparseDom = this.a.domain, - denseDom = sparseDom.parentDom; - if this.a.domain.getNNZ() >= thresh { - var count = 0; - // Normal iteration like this is more efficient than - // dense iteration, so we prefer that for the first elements - for (_, (i, j)) in zip(1..3, sparseDom) { - const idxStr = " (%?, %?)".format(i, j); // Padding to match SciPy - s += "%<16s%?\n".format(idxStr, this.a[i,j]); - } - - s += " : :\n"; // Dot dot seperator, but vertical - - // For the last elements, we iterate in dense order - // Since sparseArrays cant be strided by -1 - // We also have to do some i,j swaps for CSC vs CSR differences - count = 0; - var backString = ""; - for (i, j) in denseDom by -1 { - var row = i, col = j; - if this.matLayout==Layout.CSC { - row = j; // Iterate in Col Major Order for CSC - col = i; // To match SciPy behavior - } - if !sparseDom.contains(row, col) then continue; - const idxStr = " (%?, %?)".format(row, col); // Padding to match SciPy - backString = "%<16s%?\n".format(idxStr, this.a[row,col]) + backString; - count += 1; - if count == 3 then break; - } - s+=backString; - } else { - for (i,j) in sparseDom { - const idxStr = " (%?, %?)".format(i, j); // Padding to match SciPy - s += "%<16s%?\n".format(idxStr, this.a[i,j]); - } - } - - if this.etype == bool { - s = s.replace("true","True"); - s = s.replace("false","False"); - } - - return s; - } - - /* - Verbose flag utility method - */ - proc deinit() { - if logLevel == LogLevel.DEBUG {writeln("deinit SparseSymEntry");try! stdout.flush();} - } - } - - class GeneratorSymEntry:AbstractSymEntry { type etype; var generator: randomStream(etype); diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index 2c761713de6..a5ab484578d 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -6,6 +6,8 @@ module SparseMatrix { use CompressedSparseLayout; use CommAggregation; use CustomCopyAggregation; + use PrivateDist; + use ChapelLocks; config param aggregatedSparseMatrixCreation = true; @@ -422,8 +424,7 @@ module SparseMatrix { // If we need to improve memory usage, could we use a SparseIndexBuffer // instead? const locInds = A.domain.parentDom.localSubdomain(); - var cBlk = if !aggregatedSparseMatrixCreation then makeSparseMat(locInds, spsData) - else makeParSafeSparseMat(locInds, spsData); + var cBlk = makeSparseMat(locInds, spsData); // Stitch the local portions back together into the global-view // @@ -504,17 +505,10 @@ module SparseMatrix { errorClass="InvalidArgumentError" ); - if !aggregatedSparseMatrixCreation { - var (SD, dense) = makeSparseDomain(shape, layout); - var A: [SD] eltType; - addElementsToSparseArray(A, SD, rows, cols, vals); - return A; - } else { - var (SD, dense) = makeParSafeSparseDomain(shape, layout); - var A: [SD] eltType; - addElementsToSparseArray(A, SD, rows, cols, vals); - return A; - } + var (SD, dense) = makeSparseDomain(shape, layout); + var A: [SD] eltType; + addElementsToSparseArray(A, SD, rows, cols, vals); + return A; } proc addElementsToSparseArray(ref A, ref SD, const ref rows,const ref cols, @@ -540,13 +534,17 @@ module SparseMatrix { class DestinationHandler { var domVal; var arrVal; + var lockVal; - proc init(domVal, arrVal) { + proc init(domVal, arrVal, lockVal) { this.domVal = domVal; this.arrVal = arrVal; + this.lockVal = lockVal; } inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) { + lockVal.dsiAccess(here.id).lock(); + const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here); var locDomVal = this.domVal.locDoms[locid]!.mySparseBlock._value; var locIdxBuf = locDomVal.dsiCreateIndexBuffer(bufSize,false,false); @@ -557,37 +555,38 @@ module SparseMatrix { } locIdxBuf.commit(); - // We need a lock around this loop because another parallel 'flush' - // might be inserting indices, which will make the returned results from - // 'find' invalid. + // 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 - // - is there a benefit to avoiding the second lock by adding a version - // of 'bulkAdd' that also handles the values? - locDomVal.lockDomain(); + // - 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; } - locDomVal.unlockDomain(); + + lockVal.dsiAccess(here.id).unlock(); } } class SourceHandler { var domVal; var arrVal; + var lockVal; type elemType = (int,int,int); - proc init(D, A) { + proc init(D, A, locks) { this.domVal = D._value; this.arrVal = A._value; + this.lockVal = locks._value; } proc sourceCopy() { - return new unmanaged DestinationHandler(domVal,arrVal); + return new unmanaged DestinationHandler(domVal,arrVal, lockVal); } proc getDestinationLocale(val: elemType) { @@ -630,8 +629,9 @@ module SparseMatrix { } } } else { + var locks : [PrivateSpace] chpl_LocalSpinlock; forall (i,j,v) in zip(rows, cols, vals) - with (var agg = new CustomDstAggregator(new shared SourceHandler(SD, A))) do + with (var agg = new CustomDstAggregator(new shared SourceHandler(SD, A, locks))) do agg.copy((i,j,v)); } @@ -707,29 +707,6 @@ module SparseMatrix { sort(inds); - for ij in inds do - CDom += ij; - - var C: [CDom] int; - for ij in inds do - try! C[ij] += spsData[ij]; // TODO: Should this really throw? - - return C; - } - - // create a new sparse matrix from a map from sparse indices to values - // - proc makeParSafeSparseMat(parentDom, spsData) { - use ArkoudaSparseMatrixCompat; - use Sort; - - var CDom: sparse subdomain(parentDom) dmapped getParSafeSparseDom(Layout.CSR); - var inds: [0.. Date: Thu, 27 Aug 2026 16:35:29 -0700 Subject: [PATCH 15/17] Check correctness earlier --- benchmarks/sparse.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/benchmarks/sparse.py b/benchmarks/sparse.py index f503e1ca20e..c927018d1f4 100644 --- a/benchmarks/sparse.py +++ b/benchmarks/sparse.py @@ -122,21 +122,20 @@ def check_correctness(dtype): vals = ak.uniform(len(rows), seed=seed * 3) + 0.5 csr = create_sparse_matrix(N, rows, cols, vals, layout="CSR") - csc = create_sparse_matrix(N, cols, rows, vals, layout="CSC") - result = sparse_matrix_matrix_mult(csc, 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" From 54a6e1a23d3636f62faefe23ae353bdebd9ce81c Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Thu, 27 Aug 2026 16:40:38 -0700 Subject: [PATCH 16/17] Wrap the PrivateDist array in a class My initial attempt at using '_value' ran into an issue in PrivateDist's implementation, where it would assume that when dsiAccess was called that 'this.locale == here'. When 'flush' was initiated from a locale other than Locale0, it would grab the wrong lock, and a data race would ensue. Some kind of locale-private variable that could be declared over an array of locales would be useful here. --- src/SparseMatrix.chpl | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index a5ab484578d..d4183ff5b42 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -528,22 +528,34 @@ 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 lockVal; + var lockObj; - proc init(domVal, arrVal, lockVal) { + proc init(domVal, arrVal, lockObj) { this.domVal = domVal; this.arrVal = arrVal; - this.lockVal = lockVal; + this.lockObj = lockObj; } inline proc flush(ref rBuffer, const ref remBufferPtr, const ref myBufferIdx) { - lockVal.dsiAccess(here.id).lock(); + ref lock = lockObj.data[here.id]; + lock.lock(); const (_, locid) = this.domVal.dist.chpl__locToLocIdx(here); var locDomVal = this.domVal.locDoms[locid]!.mySparseBlock._value; @@ -569,24 +581,24 @@ module SparseMatrix { this.arrVal.locArr[locid]!.myElems._value.data[loc] = v; } - lockVal.dsiAccess(here.id).unlock(); + lock.unlock(); } } class SourceHandler { var domVal; var arrVal; - var lockVal; + var lockObj; type elemType = (int,int,int); proc init(D, A, locks) { this.domVal = D._value; this.arrVal = A._value; - this.lockVal = locks._value; + this.lockObj = locks; } proc sourceCopy() { - return new unmanaged DestinationHandler(domVal,arrVal, lockVal); + return new unmanaged DestinationHandler(domVal,arrVal, lockObj); } proc getDestinationLocale(val: elemType) { @@ -629,10 +641,11 @@ module SparseMatrix { } } } else { - var locks : [PrivateSpace] chpl_LocalSpinlock; + var locks = new unmanaged LockHelper(); 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)); + delete locks; } } From fccd68bff2bcf65d4321a232213387c500324c1e Mon Sep 17 00:00:00 2001 From: Ben Harshbarger Date: Wed, 2 Sep 2026 10:42:19 -0700 Subject: [PATCH 17/17] Respond to reviewer feedback --- benchmarks/graph_infra/arkouda.graph | 4 ++-- src/SparseMatrix.chpl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/benchmarks/graph_infra/arkouda.graph b/benchmarks/graph_infra/arkouda.graph index 64c64a55f80..c6f78591c69 100644 --- a/benchmarks/graph_infra/arkouda.graph +++ b/benchmarks/graph_infra/arkouda.graph @@ -163,11 +163,11 @@ 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 (N=1M, NNZ=10M) +graphtitle: Sparse Matrix Creation Time ylabel: Time (Seconds) perfkeys: Average CSCxCSR time = graphkeys: CSCxCSR time files: sparse.dat -graphtitle: Sparse Matrix Multiplication Time (N=1M, NNZ=10M) +graphtitle: Sparse Matrix Multiplication Time ylabel: Time (Seconds) diff --git a/src/SparseMatrix.chpl b/src/SparseMatrix.chpl index d4183ff5b42..77367085b57 100644 --- a/src/SparseMatrix.chpl +++ b/src/SparseMatrix.chpl @@ -556,6 +556,7 @@ module SparseMatrix { 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; @@ -580,8 +581,6 @@ module SparseMatrix { var (_,loc) = locDomVal.find((i,j)); this.arrVal.locArr[locid]!.myElems._value.data[loc] = v; } - - lock.unlock(); } } @@ -642,10 +641,11 @@ module SparseMatrix { } } 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)); - delete locks; } }