Baunsgaard commented on code in PR #2038:
URL: https://github.com/apache/systemds/pull/2038#discussion_r1664306304
##########
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:
this looks like a square row sum,
this is an operation you can perform directly.
##########
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];
Review Comment:
This temporary allocation seems like the only thing that limits you from
parallelizing the entire for loop.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -74,6 +79,8 @@ public class LibCommonsMath
private static final Log LOG =
LogFactory.getLog(LibCommonsMath.class.getName());
private static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-6;
private static final double EIGEN_LAMBDA = 1e-8;
+ private static final double EIGEN_EPS = Precision.EPSILON;
Review Comment:
Define this value directly without using commons math lib.
##########
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) {
+ final double c = hv[k * m + j];
+ xNormSqr += c * c;
+ }
+
+ final double a = (hv[k_kp1] > 0) ?
-FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
+
+ secondary[k] = a;
+
+ if (a != 0.0) {
+ // apply Householder transform from left and right
simultaneously
+ hv[k_kp1] -= a;
+
+ final double beta = -1 / (a * hv[k_kp1]);
+ double gamma = 0;
+
+ // compute a = beta A v, where v is the Householder
vector
+ // this loop is written in such a way
+ // 1) only the upper triangular part of the matrix is
accessed
+ // 2) access is cache-friendly for a matrix stored in
rows
+ Arrays.fill(z, k + 1, m, 0);
+ for (int i = k + 1; i < m; ++i) {
+ final double hKI = hv[km + i];
+ double zI = hv[i * m + i] * hKI;
+ for (int j = i + 1; j < m; ++j) {
+ final double hIJ = hv[i * m + j];
+ zI += hIJ * hv[k * m + j];
+ z[j] += hIJ * hKI;
+ }
+ z[i] = beta * (z[i] + zI);
+
+ gamma += z[i] * hv[km + i];
+ }
+
+ gamma *= beta / 2;
+
+ // compute z = z - gamma v
+ // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1,
k + 1, m - (k + 1));
+ for (int i = k + 1; i < m; ++i) {
+ z[i] -= gamma * hv[km + i];
+ }
+
+ // update matrix: A = A - v zT - z vT
+ // only the upper triangular part of the matrix is
updated
+ for (int i = k + 1; i < m; ++i) {
+ final double hki = hv[km + i];
+ for (int j = i; j < m; ++j) {
+ final double hkj = hv[km + j];
+
+ hv[i * m + j] -= (hki * z[j] + z[i] * hkj);
+ }
+ }
+ }
+ }
+
+ main[m - 1] = hv[(m - 1) * m + m - 1];
+
+ return householderVectors;
+ }
+
+
+ /**
+ * Computes the orthogonal matrix Q using Householder transforms.
+ * The matrix Q is built by applying Householder transforms to the
input vectors.
+ *
+ * @param hv The input vector containing the Householder vectors.
+ * @param main The main diagonal of the matrix.
+ * @param secondary The secondary diagonal of the matrix.
+ * @return The orthogonal matrix Q.
+ */
+ public static MatrixBlock getQ(final double[] hv, double[] main,
double[] secondary) {
+ final int m = main.length;
+
+ // orthogonal matrix, most operations are column major (TODO)
+ DenseBlock qaB = DenseBlockFactory.createDenseBlock(m, m);
+ double[] qaV = qaB.valuesAt(0);
+
+ // build up first part of the matrix by applying Householder
transforms
+ for (int k = m - 1; k >= 1; --k) {
Review Comment:
is there a reason why it should be done in opposite order?
##########
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) {
+ final double c = hv[k * m + j];
Review Comment:
please fix all the indentations, you can optionally use (i recommend this)
our code formatting settings file and apply it on the code parts that you have
added.
##########
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
Review Comment:
remove todos that are already verified.
##########
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) {
Review Comment:
i would like to know where the time is spend inside this function call.
is it in transformToTridiagonal, getQ or findEigenVectors?
##########
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) {
+ final double c = hv[k * m + j];
+ xNormSqr += c * c;
+ }
+
+ final double a = (hv[k_kp1] > 0) ?
-FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
+
+ secondary[k] = a;
+
+ if (a != 0.0) {
+ // apply Householder transform from left and right
simultaneously
+ hv[k_kp1] -= a;
+
+ final double beta = -1 / (a * hv[k_kp1]);
+ double gamma = 0;
+
+ // compute a = beta A v, where v is the Householder
vector
+ // this loop is written in such a way
+ // 1) only the upper triangular part of the matrix is
accessed
+ // 2) access is cache-friendly for a matrix stored in
rows
+ Arrays.fill(z, k + 1, m, 0);
+ for (int i = k + 1; i < m; ++i) {
+ final double hKI = hv[km + i];
+ double zI = hv[i * m + i] * hKI;
+ for (int j = i + 1; j < m; ++j) {
+ final double hIJ = hv[i * m + j];
+ zI += hIJ * hv[k * m + j];
+ z[j] += hIJ * hKI;
+ }
+ z[i] = beta * (z[i] + zI);
+
+ gamma += z[i] * hv[km + i];
+ }
+
+ gamma *= beta / 2;
+
+ // compute z = z - gamma v
+ // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1,
k + 1, m - (k + 1));
+ for (int i = k + 1; i < m; ++i) {
+ z[i] -= gamma * hv[km + i];
+ }
+
+ // update matrix: A = A - v zT - z vT
+ // only the upper triangular part of the matrix is
updated
+ for (int i = k + 1; i < m; ++i) {
+ final double hki = hv[km + i];
+ for (int j = i; j < m; ++j) {
+ final double hkj = hv[km + j];
+
+ hv[i * m + j] -= (hki * z[j] + z[i] * hkj);
+ }
+ }
+ }
+ }
+
+ main[m - 1] = hv[(m - 1) * m + m - 1];
+
+ return householderVectors;
+ }
+
+
+ /**
+ * Computes the orthogonal matrix Q using Householder transforms.
+ * The matrix Q is built by applying Householder transforms to the
input vectors.
+ *
+ * @param hv The input vector containing the Householder vectors.
+ * @param main The main diagonal of the matrix.
+ * @param secondary The secondary diagonal of the matrix.
+ * @return The orthogonal matrix Q.
+ */
+ public static MatrixBlock getQ(final double[] hv, double[] main,
double[] secondary) {
+ final int m = main.length;
+
+ // orthogonal matrix, most operations are column major (TODO)
+ DenseBlock qaB = DenseBlockFactory.createDenseBlock(m, m);
+ double[] qaV = qaB.valuesAt(0);
+
+ // build up first part of the matrix by applying Householder
transforms
+ for (int k = m - 1; k >= 1; --k) {
+ final int km = k * m;
+ final int km1m = (k - 1) * m;
+
+ qaV[km + k] = 1.0;
+ if (hv[km1m + k] != 0.0) {
+ final double inv = 1.0 / (secondary[k - 1] * hv[km1m +
k]);
+
+ double beta = 1.0 / secondary[k - 1];
+
+ qaV[km + k] = 1 + beta * hv[km1m + k];
+
+ // TODO: may speedup vector operations
+ for (int i = k + 1; i < m; ++i) {
+ qaV[i * m + k] = beta * hv[km1m + i];
+ }
+
+ for (int j = k + 1; j < m; ++j) {
+ beta = 0;
+ for (int i = k + 1; i < m; ++i) {
+ beta += qaV[m * i + j] * hv[km1m + i];
+ }
+ beta *= inv;
+ qaV[m * k + j] = hv[km1m + k] * beta;
+
+ for (int i = k + 1; i < m; ++i) {
+ qaV[m * i + j] += beta * hv[km1m + i];
+ }
+ }
+ }
+ }
+
+ qaV[0] = 1.0;
+ MatrixBlock res = new MatrixBlock(m, m, qaB);
+
+ return res;
+ }
+
+
+ /**
+ * Finds the eigen vectors corresponding to the given eigen values
using the Householder transformation.
+ * (Dubrulle et al., 1971).
+ *
+ * @param main The main diagonal of the tridiagonal matrix.
+ * @param secondary The secondary diagonal of the tridiagonal matrix.
+ * @param hhMatrix The Householder matrix.
+ * @param maxIter The maximum number of iterations for convergence.
+ * @param threads The number of threads to use for computation.
+ * @return An array of two MatrixBlock objects: eigen values and eigen
vectors.
+ * @throws MaxCountExceededException If the maximum number of
iterations is exceeded and convergence fails.
+ */
+ private static MatrixBlock[] findEigenVectors(double[] main, double[]
secondary, final MatrixBlock hhMatrix,
+ final int maxIter, int threads) {
+
+ DenseBlock hhDense = hhMatrix.getDenseBlock();
+
+ final int n = hhMatrix.rlen;
+ MatrixBlock eigenValues = new MatrixBlock(n, 1, main);
+
+ double[] ev = eigenValues.denseBlock.valuesAt(0);
+ final double[] e = new double[n];
+
+ System.arraycopy(secondary, 0, e, 0, n - 1);
+ e[n - 1] = 0;
+
+ // Determine the largest main and secondary value in absolute
term.
+ double maxAbsoluteValue = 0;
+ for (int i = 0; i < n; i++) {
+ maxAbsoluteValue = FastMath.max(maxAbsoluteValue,
FastMath.abs(ev[i]));
+ maxAbsoluteValue = FastMath.max(maxAbsoluteValue,
FastMath.abs(e[i]));
+ }
+ // Make null any main and secondary value too small to be
significant
+ if (maxAbsoluteValue != 0) {
+ for (int i = 0; i < n; i++) {
+ if (FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue
&& ev[i] != 0.0) {
+ ev[i] = 0;
+ }
+ if (FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue
&& e[i] != 0.0) {
+ e[i] = 0;
+ }
+ }
+ }
+
+ for (int j = 0; j < n; j++) {
+ int its = 0;
+ int m;
+ do {
Review Comment:
change the do while loops to while do. it is easier to read.
--
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
