Skip to content

[SYSTEMDS-3168] Add optimized transposed dense-sparse and sparse-dense matrix multiplication kernels. #2544

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
willy-rabe wants to merge 2 commits into
base: main
Choose a base branch
from
+605 −82
Open

[SYSTEMDS-3168] Add optimized transposed dense-sparse and sparse-dense matrix multiplication kernels. #2544

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
624b287
[SYSTEMDS-3168] sparse dense transpose kernels plus performance test
johannesmarold Jul 13, 2026
4ca22c1
[SYSTEMDS-3168] dense-sparse transpose kernels and component tests
willy-rabe Jul 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
339 changes: 339 additions & 0 deletions src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1544,6 +1544,175 @@ private static void matrixMultDenseDenseOutSparseVector(MatrixBlock m1, MatrixBl
}
}

public static void matrixMultDenseSparseMM(DenseBlock a, SparseBlock b, DenseBlock c,
boolean transA, boolean transB, int n, int cd, long xsp, int rl, int ru)
{
if(!transA && !transB){
// dispatcher parameter mismatch necessitates dummy wrappers
MatrixBlock m1 = new MatrixBlock(c.numRows(), cd, false);
m1.setDenseBlock(a);

MatrixBlock m2 = new MatrixBlock();
m2.sparseBlock = b;

MatrixBlock ret = new MatrixBlock();
ret.setDenseBlock(c);

matrixMultDenseSparseOutDense(m1, m2, ret, false, rl, ru);
}
else if(transA && !transB)
multDenseSparseTransA(a, b, c, n, cd, xsp, rl, ru);
else if(!transA && transB)
multDenseSparseTransB(a, b, c, n, cd, xsp, rl, ru);
else
multDenseSparseTransATransB(a, b, c, n, cd, xsp, rl, ru);
}

private static void multDenseSparseTransA(DenseBlock a, SparseBlock b, DenseBlock c,
int n, int cd, long xsp, int rl, int ru)
{
final int blocksizeJ = 1024;

for(int k = 0; k < cd; k++) {

if(b.isEmpty(k))
continue;

final int bpos = b.pos(k);
final int blen = b.size(k);
final int[] bix = b.indexes(k);
final double[] bvals = b.values(k);

final double[] arow = a.values(k);
final int apos = a.pos(k);

for(int bj = 0; bj < n; bj += blocksizeJ) {

int p1 = (bj == 0) ? bpos :
((b.posFIndexGTE(k, bj) >= 0) ?
bpos + b.posFIndexGTE(k, bj) :
bpos + blen);

int p2 =
((b.posFIndexGTE(k, bj + blocksizeJ) >= 0) ?
bpos + b.posFIndexGTE(k, bj + blocksizeJ) :
bpos + blen);

if(p1 >= p2)
continue;

for(int i = rl; i < ru; i++) {
final double aval = arow[apos + i];

if(aval == 0)
continue;

final double[] cvals = c.values(i);
final int cix = c.pos(i);

vectMultiplyAdd(
aval,
bvals,
cvals,
bix,
p1,
cix,
p2 - p1);
}
}
}
}

private static void multDenseSparseTransB(DenseBlock a, SparseBlock b, DenseBlock c,
int n, int cd, long xsp, int rl, int ru)
{
if( a.isContiguous() ) {
final double[] adata = a.values(0);
for(int i = rl; i < ru; i++) {
final int apos = i * cd;
final double[] cvals = c.values(i);
final int cix = c.pos(i);
for(int j = 0; j < n; j++) {
if(b.isEmpty(j))
continue;
final int bpos = b.pos(j);
final int blen = b.size(j);
final int[] bix = b.indexes(j);
final double[] bvals = b.values(j);
cvals[cix + j] = dotProduct(bvals, adata, bix, bpos, apos, blen);
}
}
}
else {
for(int i = rl; i < ru; i++) {
final double[] arow = a.values(i);
final int apos = a.pos(i);
final double[] cvals = c.values(i);
final int cix = c.pos(i);
for(int j = 0; j < n; j++) {
if(b.isEmpty(j))
continue;
final int bpos = b.pos(j);
final int blen = b.size(j);
final int[] bix = b.indexes(j);
final double[] bvals = b.values(j);
cvals[cix + j] = dotProduct(bvals, arow, bix, bpos, apos, blen);
}
}
}
}

private static void multDenseSparseTransATransB(DenseBlock a, SparseBlock b, DenseBlock c,
int n, int cd, long xsp, int rl, int ru)
{
final int m = a.numCols();
if( a.isContiguous() && c.isContiguous() ) {
final double[] adata = a.values(0);
final double[] cvals = c.values(0);
for(int j = 0; j < n; j++) {
if(b.isEmpty(j))
continue;
final int bpos = b.pos(j);
final int blen = b.size(j);
final int[] bix = b.indexes(j);
final double[] bvals = b.values(j);
for(int p = bpos; p < bpos + blen; p++) {
final int k = bix[p];
final double bval = bvals[p];
if (bval == 0)
continue;
final int apos = k * m;
for(int i = rl; i < ru; i++) {
cvals[i * n + j] += bval * adata[apos + i];
}
}
}
}
else {
for(int j = 0; j < n; j++) {
if(b.isEmpty(j))
continue;
final int bpos = b.pos(j);
final int blen = b.size(j);
final int[] bix = b.indexes(j);
final double[] bvals = b.values(j);
for(int p = bpos; p < bpos + blen; p++) {
final int k = bix[p];
final double bval = bvals[p];
if (bval == 0)
continue;
final double[] arow = a.values(k);
final int apos = a.pos(k);
for(int i = rl; i < ru; i++) {
final double[] cvals = c.values(i);
final int cix = c.pos(i);
cvals[cix + j] += bval * arow[apos + i];
}
}
}
}
}

private static void matrixMultDenseSparseOutSparse(MatrixBlock m1, MatrixBlock m2, MatrixBlock ret, boolean pm2,
int rl, int ru) {
final DenseBlock a = m1.getDenseBlock();
Expand Down Expand Up @@ -1753,6 +1922,176 @@ private static void matrixMultDenseSparseOutDenseVector(MatrixBlock m1, MatrixBl
}
}

public static void matrixMultSparseDenseMM(SparseBlock a, DenseBlock b, DenseBlock c,
boolean transA, boolean transB, int n, int cd, long xsp, int rl, int ru) {
if(!transA && !transB)
matrixMultSparseDenseMM(a, b, c, n, cd, xsp, rl, ru);
else if(transA && !transB)
multSparseDenseTransA(a, b, c, n, cd, xsp, rl, ru);
else if(!transA && transB)
multSparseDenseTransB(a, b, c, n, cd, xsp, rl, ru);
else
multSparseDenseTransATransB(a, b, c, n, cd, xsp, rl, ru);
}

private static void multSparseDenseTransA(SparseBlock a, DenseBlock b, DenseBlock c, int n, int cd, long xsp, int rl, int ru) {
final int blocksizeK = (int) (8L*xsp);
final int blocksizeJ = 1024;

for(int bk = 0; bk < cd; bk += blocksizeK) {
final int bkmin = Math.min(cd, bk + blocksizeK);

for(int bj = 0; bj < n; bj += blocksizeJ) {
final int bjlen = Math.min(n, bj + blocksizeJ) - bj;
final boolean contiguous = b.isContiguous(bk, bkmin - 1);
final double[] bvals = contiguous ? b.values(bk) : null;

for(int k = bk; k < bkmin; k++) {
if(a.isEmpty(k))
continue;

final int apos = a.pos(k);
final int alen = a.size(k);
final int[] aix = a.indexes(k);
final double[] avals = a.values(k);

int p1 = (rl == 0) ? 0 : a.posFIndexGTE(k, rl);
p1 = (p1 >= 0) ? apos + p1 : apos + alen;

int p2 = a.posFIndexGTE(k, ru);
p2 = (p2 >= 0) ? apos + p2 : apos + alen;

if(p1 >= p2)
continue;

if(contiguous) {
final int bpos = b.pos(k, bj);

for(int p = p1; p < p2; p++) {
final double aval = avals[p];
if(aval != 0) {
final int row = aix[p];
final double[] cvals = c.values(row);
final int cix = c.pos(row, bj);
vectMultiplyAdd(aval, bvals, cvals, bpos, cix, bjlen);
}
}
}
else {
final double[] kbvals = b.values(k);
final int bix = b.pos(k, bj);

for(int p = p1; p < p2; p++) {
final double aval = avals[p];
if(aval != 0) {
final int row = aix[p];
final double[] cvals = c.values(row);
final int cix = c.pos(row, bj);
vectMultiplyAdd(aval, kbvals, cvals, bix, cix, bjlen);
}
}
}
}
}
}
}

private static void multSparseDenseTransB(SparseBlock a, DenseBlock b, DenseBlock c, int n, int cd, long xsp, int rl, int ru) {
final int blocksizeK = 24;
final int blocksizeJ = Math.min(1024, (int)(8L*xsp));
final double[] bufB = new double[blocksizeK * blocksizeJ];
final int lenI = ru - rl;
final int[] p1s = new int[lenI];
final int[] p2s = new int[lenI];

for( int bk = 0; bk < cd; bk += blocksizeK ) {
final int bkmin = Math.min(cd, bk + blocksizeK);
final int bklen = bkmin - bk;

for( int i = rl; i < ru; i++ ) {
int idx = i - rl;
if( a.isEmpty(i) ) {
p1s[idx] = 0;
p2s[idx] = 0;
continue;
}
int p1 = a.posFIndexGTE(i, bk);
if( p1 < 0 ) {
p1s[idx] = 0;
p2s[idx] = 0;
continue;
}
int p2 = a.posFIndexGTE(i, bkmin);
p1s[idx] = a.pos(i) + p1;
p2s[idx] = (p2 >= 0) ? a.pos(i) + p2 : a.pos(i) + a.size(i);
}

for( int bj = 0; bj < n; bj += blocksizeJ ) {
final int bjmin = Math.min(n, bj + blocksizeJ);
final int bjlen = bjmin - bj;

for( int j = bj; j < bjmin; j++ ) {
final double[] bvals = b.values(j);
final int bpos = b.pos(j);
final int joff = j - bj;

for( int k = 0; k < bklen; k++ )
bufB[k * bjlen + joff] = bvals[bpos + bk + k];
}


for( int i = rl; i < ru; i++ ) {
int idx = i - rl;
int p1 = p1s[idx];
int p2 = p2s[idx];

if( p1 >= p2 )
continue;

final int[] aix = a.indexes(i);
final double[] avals = a.values(i);
final double[] cvals = c.values(i);
final int cix = c.pos(i, bj);

for( int p = p1; p < p2; p++ ) {
final int k = aix[p];
final double aval = avals[p];
final int koff = (k - bk) * bjlen;
vectMultiplyAdd(aval, bufB, cvals, koff, cix, bjlen);
}
}
}
}
}

private static void multSparseDenseTransATransB(SparseBlock a, DenseBlock b, DenseBlock c, int n, int cd, long xsp, int rl, int ru) {
final int blocksizeK = 24;
final int blocksizeJ = Math.min(1024, (int)(8L*xsp));

final DenseBlock tB = new DenseBlockFP64(new int[] {cd, n});
final double[] tBvals = tB.values(0);

for( int bj = 0; bj < n; bj += blocksizeJ ) {
final int bjmin = Math.min(n, bj + blocksizeJ);

for( int bk = 0; bk < cd; bk += blocksizeK ) {
final int bkmin = Math.min(cd, bk + blocksizeK);
final int bklen = bkmin - bk;

for( int j = bj; j < bjmin; j++ ) {
final double[] bvals = b.values(j);
final int bpos = b.pos(j);
final int joff = j - bj;

for( int k = 0; k < bklen; k++ )
tBvals[(bk + k) * n + bj + joff] = bvals[bpos + bk + k];
}
}
}

multSparseDenseTransA(a, tB, c, n, cd, xsp, rl, ru);
}

private static void matrixMultSparseDense(MatrixBlock m1, MatrixBlock m2, MatrixBlock ret, boolean pm2, int rl, int ru) {
SparseBlock a = m1.sparseBlock;
DenseBlock b = m2.getDenseBlock();
Expand Down
Loading