diff --git a/projects/Math/76/org/apache/commons/math/linear/SingularValueDecompositionImpl.java b/projects/Math/76/org/apache/commons/math/linear/SingularValueDecompositionImpl.java index 1436881..e418c08 100644 --- a/projects/Math/76/org/apache/commons/math/linear/SingularValueDecompositionImpl.java +++ b/projects/Math/76/org/apache/commons/math/linear/SingularValueDecompositionImpl.java @@ -159,24 +159,27 @@ public RealMatrix getU() if (m >= n) { // the tridiagonal matrix is Bt.B, where B is upper bidiagonal final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1); + eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1); final double[][] eData = e.getData(); final double[][] wData = new double[m][p]; double[] ei1 = eData[0]; - for (int i = 0; i < p - 1; ++i) { + for (int i = 0; i < p; ++i) { // compute W = B.E.S^(-1) where E is the eigenvectors matrix final double mi = mainBidiagonal[i]; final double[] ei0 = ei1; final double[] wi = wData[i]; + if (i < n - 1) { ei1 = eData[i + 1]; final double si = secondaryBidiagonal[i]; for (int j = 0; j < p; ++j) { wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; } - } + } else { for (int j = 0; j < p; ++j) { - wData[p - 1][j] = ei1[j] * mainBidiagonal[p - 1] / singularValues[j]; + wi[j] = mi * ei0[j] / singularValues[j]; } + } + } for (int i = p; i < m; ++i) { wData[i] = new double[p]; @@ -245,23 +248,26 @@ public RealMatrix getV() // the tridiagonal matrix is B.Bt, where B is lower bidiagonal // compute W = Bt.E.S^(-1) where E is the eigenvectors matrix final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1); + eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); final double[][] eData = e.getData(); final double[][] wData = new double[n][p]; double[] ei1 = eData[0]; - for (int i = 0; i < p - 1; ++i) { + for (int i = 0; i < p; ++i) { final double mi = mainBidiagonal[i]; final double[] ei0 = ei1; final double[] wi = wData[i]; + if (i < m - 1) { ei1 = eData[i + 1]; final double si = secondaryBidiagonal[i]; for (int j = 0; j < p; ++j) { wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; } - } + } else { for (int j = 0; j < p; ++j) { - wData[p - 1][j] = ei1[j] * mainBidiagonal[p - 1] / singularValues[j]; + wi[j] = mi * ei0[j] / singularValues[j]; } + } + } for (int i = p; i < n; ++i) { wData[i] = new double[p]; }