From 624b287a33cea155496fc17766c57c81e13f48ea Mon Sep 17 00:00:00 2001 From: johannesmarold Date: Mon, 13 Jul 2026 19:25:04 +0200 Subject: [PATCH 1/2] [SYSTEMDS-3168] sparse dense transpose kernels plus performance test --- .../runtime/matrix/data/LibMatrixMult.java | 170 +++++++++++++++ .../MatrixMultTransposedPerformanceTest.java | 203 +++++++++++++----- 2 files changed, 324 insertions(+), 49 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java index 9763e4ea57d..fe1b5dea821 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java @@ -1753,6 +1753,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(); diff --git a/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedPerformanceTest.java b/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedPerformanceTest.java index 94d8a5cfbbd..cf1dcd2fb0e 100644 --- a/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedPerformanceTest.java +++ b/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedPerformanceTest.java @@ -25,83 +25,188 @@ import org.junit.Test; public class MatrixMultTransposedPerformanceTest { - // could be adjusted, as it takes a lot of runtime with higher dimensions - private final int m = 200; - private final int n = 200; - private final int k = 200; @Test - public void testPerf1NoTransATransB() { - System.out.println("Case: C = A %*% t(B)"); - runTest(false, true); - System.out.println(); + public void testPerfDenseDense() throws Exception { + System.out.println("=============================================================="); + System.out.println("BENCHMARK: Dense-Dense Transposed Kernels"); + System.out.println("=============================================================="); + runBenchmarkSuite(1.0, 1.0); } @Test - public void testPerf2TransANoTransB() { - System.out.println("Case: C = t(A) %*% B"); - runTest(true, false); - System.out.println(); + public void testPerfSparseDense() throws Exception { + System.out.println("=============================================================="); + System.out.println("BENCHMARK: Sparse-Dense Transposed Kernels"); + System.out.println("=============================================================="); + for (double sp : new double[]{0.01, 0.05}) { + System.out.println(">>> Sparsity A: " + sp); + runBenchmarkSuite(sp, 1.0); + } } @Test - public void testPerf3TransATransB() { - System.out.println("Case: C = t(A) %*% t(B)"); - runTest(true, true); + public void testPerfDenseSparse() throws Exception { + System.out.println("=============================================================="); + System.out.println("BENCHMARK: Dense-Sparse Transposed Kernels"); + System.out.println("=============================================================="); + for (double sp : new double[]{0.01, 0.05}) { + System.out.println(">>> Sparsity B: " + sp); + runBenchmarkSuite(1.0, sp); + } + } + + private void runBenchmarkSuite(double spA, double spB) throws Exception { + boolean[][] transConfigs = { + {true, false}, + {false, true}, + {true, true} + }; + + int[] sizes = {200, 500}; + + for (boolean[] tc : transConfigs) { + boolean tA = tc[0]; + boolean tB = tc[1]; + String exprName = (tA && tB) ? "t(A) %*% t(B)" : tA ? "t(A) %*% B" : "A %*% t(B)"; + + System.out.printf("--- Case: C = %s ---%n", exprName); + + for (int size : sizes) { + System.out.printf("Size: %d%n", size); + runTest(spA, spB, tA, tB, size, size, size); + } + System.out.println(); + } } - private void runTest(boolean tA, boolean tB) { - int REP = 100; + private void runTest(double spA, double spB, boolean tA, boolean tB, int m, int n, int k) throws Exception { + final int REP = 100; + final int WARMUP = 50; - // setup Dimensions int rowsA = tA ? k : m; int colsA = tA ? m : k; int rowsB = tB ? n : k; int colsB = tB ? k : n; - // generate random matrices - MatrixBlock A = MatrixBlock.randOperations(rowsA, colsA, 1.0, -1, 1, "uniform", 7); - MatrixBlock B = MatrixBlock.randOperations(rowsB, colsB, 1.0, -1, 1, "uniform", 3); - MatrixBlock C = new MatrixBlock(m, n, false); - C.allocateDenseBlock(); + MatrixBlock ma = generateInput(rowsA, colsA, spA, 7); + MatrixBlock mb = generateInput(rowsB, colsB, spB, 3); + MatrixBlock mc = new MatrixBlock(m, n, false); + mc.allocateDenseBlock(); - for(int i=0; i<50; i++) { - runOldMethod(A, B, tA, tB); - runNewKernel(A, B, C, tA, tB); + for (int i = 0; i < WARMUP; i++) { + runOldMethod(ma, mb, tA, tB); + runNewKernel(ma, mb, mc, tA, tB); } - // Measure Old Method - long startTimeOld = System.nanoTime(); - for(int i = 0; i < REP; i++) { - runOldMethod(A, B, tA, tB); + long startOld = System.nanoTime(); + for (int i = 0; i < REP; i++) { + runOldMethod(ma, mb, tA, tB); } - double avgTimeOld = (System.nanoTime() - startTimeOld) / 1e6 / REP; + double avgOld = (System.nanoTime() - startOld) / 1e6 / REP; - // Measure New Kernel - double startTimeNew = System.nanoTime(); - for(int i = 0; i < REP; i++) { - runNewKernel(A, B, C, tA, tB); + long startNew = System.nanoTime(); + for (int i = 0; i < REP; i++) { + runNewKernel(ma, mb, mc, tA, tB); } - double avgTimeNew = (System.nanoTime() - startTimeNew) / 1e6 / REP; + double avgNew = (System.nanoTime() - startNew) / 1e6 / REP; - // print results comparison - System.out.printf("Old Method: %.3f ms | New Kernel: %.3f ms%n", avgTimeOld, avgTimeNew); + System.out.printf(" Old Method: %.3f ms | New Kernel: %.3f ms | Speedup: %.2fx%n", + avgOld, avgNew, avgOld / Math.max(1e-9, avgNew)); } - private void runNewKernel(MatrixBlock A, MatrixBlock B, MatrixBlock C, boolean tA, boolean tB) { - C.reset(); - LibMatrixMult.matrixMultDenseDenseMM(A.getDenseBlock(), B.getDenseBlock(), C.getDenseBlock(), tA, tB, m, k, 0, m, 0, n); + private MatrixBlock generateInput(int rows, int cols, double sparsity, long seed) { + MatrixBlock mb = MatrixBlock.randOperations(rows, cols, sparsity, -1, 1, "uniform", seed); + mb.examSparsity(); + if (sparsity < 1.0) { + if (!mb.isInSparseFormat()) + mb.denseToSparse(true); + if (mb.getSparseBlock() == null) + mb.allocateSparseRowsBlock(); + } + return mb; } - private void runOldMethod(MatrixBlock A, MatrixBlock B, boolean tA, boolean tB) { - // do transpose if needed - MatrixBlock A_in = tA ? LibMatrixReorg.transpose(A) : A; - MatrixBlock B_in = tB ? LibMatrixReorg.transpose(B) : B; + private void runNewKernel(MatrixBlock ma, MatrixBlock mb, MatrixBlock mc, boolean tA, boolean tB) throws Exception { + mc.reset(); + mc.allocateDenseBlock(); + + boolean sparseA = ma.isInSparseFormat(); + boolean sparseB = mb.isInSparseFormat(); - MatrixBlock C = new MatrixBlock(m, n, false); - C.allocateDenseBlock(); + int m = tA ? ma.getNumColumns() : ma.getNumRows(); + int n = tB ? mb.getNumRows() : mb.getNumColumns(); + int k = tA ? ma.getNumRows() : ma.getNumColumns(); - LibMatrixMult.matrixMultDenseDenseMM(A_in.getDenseBlock(), B_in.getDenseBlock(), C.getDenseBlock(), false, - false, m, k, 0, m, 0, n); + if (!sparseA && !sparseB) { + LibMatrixMult.matrixMultDenseDenseMM( + ma.getDenseBlock(), + mb.getDenseBlock(), + mc.getDenseBlock(), + tA, tB, n, k, 0, m, 0, n + ); + } + else if (sparseA && !sparseB) { + int cd = tB ? mb.getNumColumns() : mb.getNumRows(); + long xsp = (long) m * cd / Math.max(1L, ma.getNonZeros()); + LibMatrixMult.matrixMultSparseDenseMM( + ma.getSparseBlock(), + mb.getDenseBlock(), + mc.getDenseBlock(), + tA, tB, n, cd, xsp, 0, m + ); + } + else if (!sparseA && sparseB) { + long xsp = (long) ma.getNumRows() * ma.getNumColumns() / Math.max(1L, ma.getNonZeros()); + LibMatrixMult.matrixMultDenseSparseMM( + ma.getDenseBlock(), + mb.getSparseBlock(), + mc.getDenseBlock(), + tA, tB, n, k, xsp, 0, m + ); + } + mc.recomputeNonZeros(); + } + + private void runOldMethod(MatrixBlock ma, MatrixBlock mb, boolean tA, boolean tB) throws Exception { + MatrixBlock A_in = tA ? LibMatrixReorg.transpose(ma) : ma; + MatrixBlock B_in = tB ? LibMatrixReorg.transpose(mb) : mb; + + boolean sparseA = A_in.isInSparseFormat(); + boolean sparseB = B_in.isInSparseFormat(); + + int m = A_in.getNumRows(); + int n = B_in.getNumColumns(); + int k = A_in.getNumColumns(); + + MatrixBlock mc = new MatrixBlock(m, n, false); + mc.allocateDenseBlock(); + + if (!sparseA && !sparseB) { + LibMatrixMult.matrixMultDenseDenseMM( + A_in.getDenseBlock(), + B_in.getDenseBlock(), + mc.getDenseBlock(), + false, false, n, k, 0, m, 0, n + ); + } + else if (sparseA && !sparseB) { + long xsp = (long) m * k / Math.max(1L, A_in.getNonZeros()); + LibMatrixMult.matrixMultSparseDenseMM( + A_in.getSparseBlock(), + B_in.getDenseBlock(), + mc.getDenseBlock(), + false, false, n, k, xsp, 0, m + ); + } + else if (!sparseA && sparseB) { + long xsp = (long) A_in.getNumRows() * A_in.getNumColumns() / Math.max(1L, A_in.getNonZeros()); + LibMatrixMult.matrixMultDenseSparseMM( + A_in.getDenseBlock(), + B_in.getSparseBlock(), + mc.getDenseBlock(), + false, false, n, k, xsp, 0, m + ); + } } } From 4ca22c1fb48074cc980c5c9ab7fbcdcd052d4a5f Mon Sep 17 00:00:00 2001 From: Willy Rabe Date: Tue, 14 Jul 2026 16:34:06 +0200 Subject: [PATCH 2/2] [SYSTEMDS-3168] dense-sparse transpose kernels and component tests --- .../runtime/matrix/data/LibMatrixMult.java | 169 ++++++++++++++++++ .../matrixmult/MatrixMultTransposedTest.java | 145 +++++++++++---- 2 files changed, 281 insertions(+), 33 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java index fe1b5dea821..d411eabac00 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixMult.java @@ -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(); diff --git a/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedTest.java b/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedTest.java index 2b9f09ba13b..d6ba0b89f81 100644 --- a/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedTest.java +++ b/src/test/java/org/apache/sysds/test/component/matrixmult/MatrixMultTransposedTest.java @@ -19,73 +19,152 @@ package org.apache.sysds.test.component.matrixmult; -import org.apache.sysds.runtime.data.DenseBlock; +import java.util.Random; + import org.apache.sysds.runtime.matrix.data.LibMatrixMult; import org.apache.sysds.runtime.matrix.data.LibMatrixReorg; import org.apache.sysds.runtime.matrix.data.MatrixBlock; import org.apache.sysds.test.TestUtils; import org.junit.Test; -import java.util.Random; - public class MatrixMultTransposedTest { - // run multiple random scenarios @Test - public void testCaseNoTransATransB() { - for(int i=0; i<10; i++) { - runTest(false, true); - } + public void testDenseDenseTransA() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, false, true, false); } @Test - public void testCaseTransANoTransB() { - for(int i=0; i<10; i++) { - runTest(true, false); - } + public void testDenseDenseTransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, false, false, true); } @Test - public void testCaseTransATransB() { - for(int i=0; i<10; i++) { - runTest(true, true); - } + public void testDenseDenseTransATransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, false, true, true); } - private void runTest(boolean tA, boolean tB) { - Random rand = new Random(); + @Test + public void testSparseDenseTransA() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(true, false, true, false); + } + + @Test + public void testSparseDenseTransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(true, false, false, true); + } + + @Test + public void testSparseDenseTransATransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(true, false, true, true); + } + + @Test + public void testDenseSparseTransA() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, true, true, false); + } + + @Test + public void testDenseSparseTransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, true, false, true); + } - // generate random dimensions between 1 and 300 + @Test + public void testDenseSparseTransATransB() throws Exception { + for(int i = 0; i < 10; i++) + runRandomTest(false, true, true, true); + } + + private void runRandomTest(boolean sparseA, boolean sparseB, boolean tA, boolean tB) throws Exception { + Random rand = new Random(); int m = rand.nextInt(300) + 1; int n = rand.nextInt(300) + 1; int k = rand.nextInt(300) + 1; + double spA = sparseA ? 0.05 : 1.0; + double spB = sparseB ? 0.05 : 1.0; + + runTest(spA, spB, tA, tB, m, n, k); + } + private void runTest(double spA, double spB, boolean tA, boolean tB, int m, int n, int k) throws Exception { int rowsA = tA ? k : m; int colsA = tA ? m : k; int rowsB = tB ? n : k; int colsB = tB ? k : n; - MatrixBlock ma = MatrixBlock.randOperations(rowsA, colsA, 1.0, -1, 1, "uniform", 7); - MatrixBlock mb = MatrixBlock.randOperations(rowsB, colsB, 1.0, -1, 1, "uniform", 3); - + MatrixBlock ma = generateInput(rowsA, colsA, spA, 7); + MatrixBlock mb = generateInput(rowsB, colsB, spB, 3); MatrixBlock mc = new MatrixBlock(m, n, false); mc.allocateDenseBlock(); - DenseBlock a = ma.getDenseBlock(); - DenseBlock b = mb.getDenseBlock(); - DenseBlock c = mc.getDenseBlock(); + runNewKernel(ma, mb, mc, tA, tB); - LibMatrixMult.matrixMultDenseDenseMM(a, b, c, tA, tB, n, k, 0, m, 0, n); + MatrixBlock A_in = tA ? LibMatrixReorg.transpose(ma) : ma; + MatrixBlock B_in = tB ? LibMatrixReorg.transpose(mb) : mb; + MatrixBlock expected = LibMatrixMult.matrixMult(A_in, B_in); - mc.recomputeNonZeros(); + TestUtils.compareMatrices(expected, mc, 1e-8); + } - // calc true result with existing methods - MatrixBlock ma_in = tA ? LibMatrixReorg.transpose(ma) : ma; - MatrixBlock mb_in = tB ? LibMatrixReorg.transpose(mb) : mb; - MatrixBlock expected = LibMatrixMult.matrixMult(ma_in, mb_in); + private MatrixBlock generateInput(int rows, int cols, double sparsity, long seed) { + MatrixBlock mb = MatrixBlock.randOperations(rows, cols, sparsity, -1, 1, "uniform", seed); + mb.examSparsity(); + if (sparsity < 1.0) { + if (!mb.isInSparseFormat()) + mb.denseToSparse(true); + if (mb.getSparseBlock() == null) + mb.allocateSparseRowsBlock(); + } + return mb; + } - // compare results - TestUtils.compareMatrices(expected, mc, 1e-8); + private void runNewKernel(MatrixBlock ma, MatrixBlock mb, MatrixBlock mc, boolean tA, boolean tB) throws Exception { + mc.reset(); + mc.allocateDenseBlock(); + + boolean sparseA = ma.isInSparseFormat(); + boolean sparseB = mb.isInSparseFormat(); + + int m = tA ? ma.getNumColumns() : ma.getNumRows(); + int n = tB ? mb.getNumRows() : mb.getNumColumns(); + int k = tA ? ma.getNumRows() : ma.getNumColumns(); + + if (!sparseA && !sparseB) { + LibMatrixMult.matrixMultDenseDenseMM( + ma.getDenseBlock(), + mb.getDenseBlock(), + mc.getDenseBlock(), + tA, tB, n, k, 0, m, 0, n + ); + } + else if (sparseA && !sparseB) { + int cd = tB ? mb.getNumColumns() : mb.getNumRows(); + long xsp = (long) m * cd / Math.max(1L, ma.getNonZeros()); + LibMatrixMult.matrixMultSparseDenseMM( + ma.getSparseBlock(), + mb.getDenseBlock(), + mc.getDenseBlock(), + tA, tB, n, cd, xsp, 0, m + ); + } + else if (!sparseA && sparseB) { + long xsp = (long) ma.getNumRows() * ma.getNumColumns() / Math.max(1L, ma.getNonZeros()); + LibMatrixMult.matrixMultDenseSparseMM( + ma.getDenseBlock(), + mb.getSparseBlock(), + mc.getDenseBlock(), + tA, tB, n, k, xsp, 0, m + ); + } + mc.recomputeNonZeros(); } }