Baunsgaard commented on code in PR #2038:
URL: https://github.com/apache/systemds/pull/2038#discussion_r1668973473
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -209,6 +218,339 @@ private static MatrixBlock computeSolve(MatrixBlock in1,
MatrixBlock in2) {
return DataConverter.convertToMatrixBlock(solutionMatrix);
}
+
+ /**
+ * Computes the eigen decomposition of a symmetric matrix.
+ *
+ * @param in The input matrix to compute the eigen decomposition on.
+ * @param threads The number of threads to use for computation.
+ * @return An array of MatrixBlock objects containing the real eigen
values and eigen vectors.
+ */
+ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock
in, int threads) {
+
+ // TODO: Verify matrix is symmetric
+ final double[] mainDiag = new double[in.rlen];
+ final double[] secDiag = new double[in.rlen - 1];
+
+ final MatrixBlock householderVectors =
transformToTridiagonal(in, mainDiag, secDiag, threads);
+
+ // TODO: Consider using sparse blocks
+ final double[] hv = householderVectors.getDenseBlockValues();
+ MatrixBlock houseHolderMatrix = getQ(hv, mainDiag, secDiag);
+
+ MatrixBlock[] evResult = findEigenVectors(mainDiag, secDiag,
houseHolderMatrix, EIGEN_MAX_ITER, threads);
+
+ MatrixBlock realEigenValues = evResult[0];
+ MatrixBlock eigenVectors = evResult[1];
+
+ realEigenValues.setNonZeros(realEigenValues.rlen *
realEigenValues.rlen);
+ eigenVectors.setNonZeros(eigenVectors.rlen * eigenVectors.clen);
+
+ return new MatrixBlock[] { realEigenValues, eigenVectors };
+ }
+
+ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix,
double[] main, double[] secondary, int threads) {
+ final int m = matrix.rlen;
+
+ // MatrixBlock householderVectors = ;
+ MatrixBlock householderVectors = matrix.extractTriangular(new
MatrixBlock(m, m, false), false, true, true);
+ if (householderVectors.isInSparseFormat()) {
+ householderVectors.sparseToDense(threads);
+ }
+
+ final double[] z = new double[m];
+
+ // TODO: Consider sparse block case
+ final double[] hv = householderVectors.getDenseBlockValues();
+
+ for (int k = 0; k < m - 1; k++) {
+ final int k_kp1 = k * m + k + 1;
+ final int km = k * m;
+
+ // zero-out a row and a column simultaneously
+ main[k] = hv[km + k];
+
+ // TODO: may or may not improve, seems it does not
+ // double xNormSqr = householderVectors.slice(k, k, k +
1, m - 1).sumSq();
+ // double xNormSqr = LibMatrixMult.dotProduct(hv, hv,
k_kp1, k_kp1, m - (k + 1));
+ double xNormSqr = 0;
+ for (int j = k + 1; j < m; ++j) {
Review Comment:
okay, then we do no change it.
--
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
