jessicapriebe commented on code in PR #2199:
URL: https://github.com/apache/systemds/pull/2199#discussion_r1937383436
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixReorg.java:
##########
@@ -4227,4 +4234,231 @@ public Object call() {
return null;
}
}
+
+ /**
+ * Transposes a dense matrix in-place using following cycles based on
Brenner's method. This
+ * method shifts cycles with a focus on less storage by using cycle
leaders based on prime factorization. The used
+ * storage is in O(n+m). Quadratic matrices should be handled outside
this method (using the trivial method) for a
+ * speedup. This method is based on: Algorithm 467, Brenner,
https://dl.acm.org/doi/pdf/10.1145/355611.362542.
+ *
+ * @param in The input matrix to be transposed.
+ * @param k The number of threads.
+ */
+ public static void transposeInPlaceDenseBrenner(MatrixBlock in, int k) {
+
+ DenseBlock denseBlock = in.getDenseBlock();
+ double[] matrix = in.getDenseBlockValues();
+
+ final int rows = in.getNumRows();
+ final int cols = in.getNumColumns();
+
+ // Brenner: rows + cols / 2 is sufficient for most cases.
+ int workSize = rows + cols;
+ int maxIndex = rows * cols - 1;
+
+ // prime factorization of the maximum index to identify cycle
structures
+ // Brenner: length 8 is sufficient up to maxIndex 2*3*5*...*19
= 9,767,520.
+ int[] primes = new int[8];
+ int[] exponents = new int[8];
+ int[] powers = new int[8];
+ int numPrimes = primeFactorization(maxIndex, primes, exponents,
powers);
+
+ int[] iExponents = new int[numPrimes];
+ int div = 1;
+
+ div:
+ while(div < maxIndex / 2) {
+
+ // number of indices divisible by div and no other
divisor of maxIndex
+ int count = eulerTotient(primes, exponents, iExponents,
numPrimes, maxIndex / div);
+ // all false
+ boolean[] moved = new boolean[workSize];
+ // starting point cycle
+ int start = div;
+
+ count:
+ do {
+ // companion of start
+ int comp = maxIndex - start;
+
+ if(start == div) {
+ // shift cycles
+ count = simultaneousCycleShift(matrix,
moved, rows, maxIndex, count, workSize, start, comp);
+ start += div;
+ }
+ else if(start < workSize && moved[start]) {
+ // already moved
+ start += div;
+ }
+ else {
+ // handle other cycle starts
+ int cycleLeader = start / div;
+ for(int ip = 0; ip < numPrimes; ip++) {
+ if(iExponents[ip] !=
exponents[ip] && cycleLeader % primes[ip] == 0) {
+ start += div;
+ continue count;
+ }
+ }
+
+ if(start < workSize) {
+ count =
simultaneousCycleShift(matrix, moved, rows, maxIndex, count, workSize, start,
comp);
+ start += div;
+ continue;
+ }
+
+ int test = start;
+ do {
+ test = prevIndexCycle(test,
rows, cols);
+ if(test < start || test > comp)
{
+ start += div;
+ continue count;
+ }
+ }
+ while(test > start && test < comp);
+
+ count = simultaneousCycleShift(matrix,
moved, rows, maxIndex, count, workSize, start, comp);
+ start += div;
+ }
+ }
+ while(count > 0);
+
+ // update cycle divisor for the next set of cycles
based on prime factors
+ for(int ip = 0; ip < numPrimes; ip++) {
+ if(iExponents[ip] != exponents[ip]) {
+ iExponents[ip]++;
+ div *= primes[ip];
+ continue div;
+ }
+ iExponents[ip] = 0;
+ div /= powers[ip];
+ }
+ break;
+ }
+
+ denseBlock.setDims(new int[] {cols, rows});
+ in.setNumColumns(rows);
+ in.setNumRows(cols);
+ }
+
+ /**
+ * Performs a simultaneous cycle shift for a cycle and its companion
cycle. This method ensures that distinct cycles
+ * or self-dual cycles are handled correctly. This method is based on:
Algorithm 2, Karlsson,
+ * https://webapps.cs.umu.se/uminf/reports/2009/011/part1.pdf and
Algorithm 467, Brenner,
+ * https://dl.acm.org/doi/pdf/10.1145/355611.362542.
+ *
+ * @param matrix The matrix whose elements are being shifted.
+ * @param moved Boolean array tracking whether an element has
already been moved.
+ * @param rows The number of rows in the matrix.
+ * @param maxIndex The maximum valid index in the matrix.
+ * @param count The number of elements left to process.
+ * @param workSize The length of moved.
+ * @param start The starting index for the cycle shift.
+ * @param comp The corresponding companion index.
+ * @return The updated count of elements remaining to shift.
+ */
+ private static int simultaneousCycleShift(double[] matrix, boolean[]
moved, int rows, int maxIndex, int count,
+ int workSize, int start, int comp) {
+
+ int orig = start;
+ double val = matrix[orig];
+ double cval = matrix[comp];
+
+ int prevOrig = prevIndexCycle(orig, rows, (maxIndex + 1) /
rows);
+ int prevComp = maxIndex - prevOrig;
+
+ while(true) {
Review Comment:
thank you for your feedback! I've addressed your comments in the latest
commit :)
regarding parallelization, unfortunately, it's not possible at this level.
it would need to be implemented at a higher level, specifically for the div
loop.
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: dev-unsubscr...@systemds.apache.org
For queries about this service, please contact Infrastructure at:
us...@infra.apache.org
