Big THANKS to Chris and Luc for this! Its great to see those bugs resolved.
Phil On 7/20/11 5:15 AM, l...@apache.org wrote: > Author: luc > Date: Wed Jul 20 12:15:00 2011 > New Revision: 1148714 > > URL: http://svn.apache.org/viewvc?rev=1148714&view=rev > Log: > Rewritten SVD implementation based on JAMA code. > > JIRA: MATH-327, MATH-383, MATH-465, MATH-583, MATH-611 > > Added: > > commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/ > > commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/matrix1.csv > > commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/matrix2.csv > Modified: > commons/proper/math/trunk/pom.xml > > commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java > commons/proper/math/trunk/src/site/xdoc/changes.xml > > commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java > > Modified: commons/proper/math/trunk/pom.xml > URL: > http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1148714&r1=1148713&r2=1148714&view=diff > ============================================================================== > --- commons/proper/math/trunk/pom.xml (original) > +++ commons/proper/math/trunk/pom.xml Wed Jul 20 12:15:00 2011 > @@ -196,6 +196,9 @@ > <name>Thomas Neidhart</name> > </contributor> > <contributor> > + <name>Christopher Nix</name> > + </contributor> > + <contributor> > <name>Fredrik Norin</name> > </contributor> > <contributor> > > Modified: > commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java > URL: > http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java?rev=1148714&r1=1148713&r2=1148714&view=diff > ============================================================================== > --- > commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java > (original) > +++ > commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java > Wed Jul 20 12:15:00 2011 > @@ -14,7 +14,6 @@ > * See the License for the specific language governing permissions and > * limitations under the License. > */ > - > package org.apache.commons.math.linear; > > import org.apache.commons.math.exception.NumberIsTooLargeException; > @@ -35,23 +34,38 @@ import org.apache.commons.math.util.Fast > * @since 2.0 > */ > public class SingularValueDecompositionImpl implements > SingularValueDecomposition { > - /** Number of rows of the initial matrix. */ > + > + /** Relative threshold for small singular values. */ > + private static final double EPS = 0x1.0p-52; > + > + /** Absolute threshold for small singular values. */ > + private static final double TINY = 0x1.0p-966; > + > + /** Computed singular values. */ > + private double[] singularValues; > + > + /** Row dimension. */ > private int m; > - /** Number of columns of the initial matrix. */ > + > + /** Column dimension. */ > private int n; > - /** Eigen decomposition of the tridiagonal matrix. */ > - private EigenDecomposition eigenDecomposition; > - /** Singular values. */ > - private double[] singularValues; > - /** Cached value of U. */ > + > + /** Indicator for transposed matrix. */ > + private boolean transposed; > + > + /** Cached value of U matrix. */ > private RealMatrix cachedU; > - /** Cached value of U<sup>T</sup>. */ > + > + /** Cached value of transposed U matrix. */ > private RealMatrix cachedUt; > - /** Cached value of S. */ > + > + /** Cached value of S (diagonal) matrix. */ > private RealMatrix cachedS; > - /** Cached value of V. */ > + > + /** Cached value of V matrix. */ > private RealMatrix cachedV; > - /** Cached value of V<sup>T</sup>. */ > + > + /** Cached value of transposed V matrix. */ > private RealMatrix cachedVt; > > /** > @@ -60,80 +74,397 @@ public class SingularValueDecompositionI > * @param matrix Matrix to decompose. > */ > public SingularValueDecompositionImpl(final RealMatrix matrix) { > + > + // Derived from LINPACK code. > + // Initialize. > + double[][] A; > m = matrix.getRowDimension(); > n = matrix.getColumnDimension(); > - > - cachedU = null; > - cachedS = null; > - cachedV = null; > - cachedVt = null; > - > - double[][] localcopy = matrix.getData(); > - double[][] matATA = new double[n][n]; > - // > - // create A^T*A > - // > - for (int i = 0; i < n; i++) { > - for (int j = i; j < n; j++) { > - matATA[i][j] = 0.0; > - for (int k = 0; k < m; k++) { > - matATA[i][j] += localcopy[k][i] * localcopy[k][j]; > - } > - matATA[j][i] = matATA[i][j]; > + if (matrix.getRowDimension() < matrix.getColumnDimension()) { > + transposed = true; > + A = matrix.transpose().getData(); > + m = matrix.getColumnDimension(); > + n = matrix.getRowDimension(); > + } else { > + transposed = false; > + A = matrix.getData(); > + m = matrix.getRowDimension(); > + n = matrix.getColumnDimension(); > + } > + int nu = FastMath.min(m, n); > + singularValues = new double[FastMath.min(m + 1, n)]; > + double[][] U = new double[m][nu]; > + double[][] V = new double[n][n]; > + double[] e = new double[n]; > + double[] work = new double[m]; > + boolean wantu = true; > + boolean wantv = true; > + // Reduce A to bidiagonal form, storing the diagonal elements > + // in s and the super-diagonal elements in e. > + int nct = FastMath.min(m - 1, n); > + int nrt = FastMath.max(0, FastMath.min(n - 2, m)); > + for (int k = 0; k < FastMath.max(nct, nrt); k++) { > + if (k < nct) { > + // Compute the transformation for the k-th column and > + // place the k-th diagonal in s[k]. > + // Compute 2-norm of k-th column without under/overflow. > + singularValues[k] = 0; > + for (int i = k; i < m; i++) { > + singularValues[k] = FastMath.hypot(singularValues[k], > A[i][k]); > + } > + if (singularValues[k] != 0.0) { > + if (A[k][k] < 0.0) { > + singularValues[k] = -singularValues[k]; > + } > + for (int i = k; i < m; i++) { > + A[i][k] /= singularValues[k]; > + } > + A[k][k] += 1.0; > + } > + singularValues[k] = -singularValues[k]; > + } > + for (int j = k + 1; j < n; j++) { > + if ((k < nct) & (singularValues[k] != 0.0)) { > + // Apply the transformation. > + double t = 0; > + for (int i = k; i < m; i++) { > + t += A[i][k] * A[i][j]; > + } > + t = -t / A[k][k]; > + for (int i = k; i < m; i++) { > + A[i][j] += t * A[i][k]; > + } > + } > + // Place the k-th row of A into e for the > + // subsequent calculation of the row transformation. > + e[j] = A[k][j]; > + } > + if (wantu & (k < nct)) { > + // Place the transformation in U for subsequent back > + // multiplication. > + for (int i = k; i < m; i++) { > + U[i][k] = A[i][k]; > + } > + } > + if (k < nrt) { > + // Compute the k-th row transformation and place the > + // k-th super-diagonal in e[k]. > + // Compute 2-norm without under/overflow. > + e[k] = 0; > + for (int i = k + 1; i < n; i++) { > + e[k] = FastMath.hypot(e[k], e[i]); > + } > + if (e[k] != 0.0) { > + if (e[k + 1] < 0.0) { > + e[k] = -e[k]; > + } > + for (int i = k + 1; i < n; i++) { > + e[i] /= e[k]; > + } > + e[k + 1] += 1.0; > + } > + e[k] = -e[k]; > + if ((k + 1 < m) & (e[k] != 0.0)) { > + // Apply the transformation. > + for (int i = k + 1; i < m; i++) { > + work[i] = 0.0; > + } > + for (int j = k + 1; j < n; j++) { > + for (int i = k + 1; i < m; i++) { > + work[i] += e[j] * A[i][j]; > + } > + } > + for (int j = k + 1; j < n; j++) { > + double t = -e[j] / e[k + 1]; > + for (int i = k + 1; i < m; i++) { > + A[i][j] += t * work[i]; > + } > + } > + } > + if (wantv) { > + // Place the transformation in V for subsequent > + // back multiplication. > + for (int i = k + 1; i < n; i++) { > + V[i][k] = e[i]; > + } > + } > + } > + } > + // Set up the final bidiagonal matrix or order p. > + int p = FastMath.min(n, m + 1); > + if (nct < n) { > + singularValues[nct] = A[nct][nct]; > + } > + if (m < p) { > + singularValues[p - 1] = 0.0; > + } > + if (nrt + 1 < p) { > + e[nrt] = A[nrt][p - 1]; > + } > + e[p - 1] = 0.0; > + // If required, generate U. > + if (wantu) { > + for (int j = nct; j < nu; j++) { > + for (int i = 0; i < m; i++) { > + U[i][j] = 0.0; > + } > + U[j][j] = 1.0; > + } > + for (int k = nct - 1; k >= 0; k--) { > + if (singularValues[k] != 0.0) { > + for (int j = k + 1; j < nu; j++) { > + double t = 0; > + for (int i = k; i < m; i++) { > + t += U[i][k] * U[i][j]; > + } > + t = -t / U[k][k]; > + for (int i = k; i < m; i++) { > + U[i][j] += t * U[i][k]; > + } > + } > + for (int i = k; i < m; i++) { > + U[i][k] = -U[i][k]; > + } > + U[k][k] = 1.0 + U[k][k]; > + for (int i = 0; i < k - 1; i++) { > + U[i][k] = 0.0; > + } > + } else { > + for (int i = 0; i < m; i++) { > + U[i][k] = 0.0; > + } > + U[k][k] = 1.0; > + } > + } > + } > + // If required, generate V. > + if (wantv) { > + for (int k = n - 1; k >= 0; k--) { > + if ((k < nrt) & (e[k] != 0.0)) { > + for (int j = k + 1; j < nu; j++) { > + double t = 0; > + for (int i = k + 1; i < n; i++) { > + t += V[i][k] * V[i][j]; > + } > + t = -t / V[k + 1][k]; > + for (int i = k + 1; i < n; i++) { > + V[i][j] += t * V[i][k]; > + } > + } > + } > + for (int i = 0; i < n; i++) { > + V[i][k] = 0.0; > + } > + V[k][k] = 1.0; > + } > + } > + // Main iteration loop for the singular values. > + int pp = p - 1; > + int iter = 0; > + while (p > 0) { > + int k, kase; > + // Here is where a test for too many iterations would go. > + // This section of the program inspects for > + // negligible elements in the s and e arrays. On > + // completion the variables kase and k are set as follows. > + // kase = 1 if s(p) and e[k-1] are negligible and k<p > + // kase = 2 if s(k) is negligible and k<p > + // kase = 3 if e[k-1] is negligible, k<p, and > + // s(k), ..., s(p) are not negligible (qr step). > + // kase = 4 if e(p-1) is negligible (convergence). > + for (k = p - 2; k >= -1; k--) { > + if (k == -1) { > + break; > + } > + final double threshold = > + TINY + EPS * (FastMath.abs(singularValues[k]) + > FastMath.abs(singularValues[k + 1])); > + if (FastMath.abs(e[k]) <= threshold) { > + e[k] = 0.0; > + break; > + } > + } > + if (k == p - 2) { > + kase = 4; > + } else { > + int ks; > + for (ks = p - 1; ks >= k; ks--) { > + if (ks == k) { > + break; > + } > + double t = (ks != p ? FastMath.abs(e[ks]) : 0.0) + > + (ks != k + 1 ? FastMath.abs(e[ks - 1]) : 0.0); > + if (FastMath.abs(singularValues[ks]) <= TINY + EPS * t) { > + singularValues[ks] = 0.0; > + break; > + } > + } > + if (ks == k) { > + kase = 3; > + } else if (ks == p - 1) { > + kase = 1; > + } else { > + kase = 2; > + k = ks; > + } > + } > + k++; > + // Perform the task indicated by kase. > + switch (kase) { > + // Deflate negligible s(p). > + case 1: { > + double f = e[p - 2]; > + e[p - 2] = 0.0; > + for (int j = p - 2; j >= k; j--) { > + double t = FastMath.hypot(singularValues[j], f); > + double cs = singularValues[j] / t; > + double sn = f / t; > + singularValues[j] = t; > + if (j != k) { > + f = -sn * e[j - 1]; > + e[j - 1] = cs * e[j - 1]; > + } > + if (wantv) { > + for (int i = 0; i < n; i++) { > + t = cs * V[i][j] + sn * V[i][p - 1]; > + V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - > 1]; > + V[i][j] = t; > + } > + } > + } > + } > + break; > + // Split at negligible s(k). > + case 2: { > + double f = e[k - 1]; > + e[k - 1] = 0.0; > + for (int j = k; j < p; j++) { > + double t = FastMath.hypot(singularValues[j], f); > + double cs = singularValues[j] / t; > + double sn = f / t; > + singularValues[j] = t; > + f = -sn * e[j]; > + e[j] = cs * e[j]; > + if (wantu) { > + for (int i = 0; i < m; i++) { > + t = cs * U[i][j] + sn * U[i][k - 1]; > + U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - > 1]; > + U[i][j] = t; > + } > + } > + } > + } > + break; > + // Perform one qr step. > + case 3: { > + // Calculate the shift. > + double scale = > FastMath.max(FastMath.max(FastMath.max(FastMath.max( > + FastMath.abs(singularValues[p - 1]), > FastMath.abs(singularValues[p - 2])), FastMath.abs(e[p - 2])), > + FastMath.abs(singularValues[k])), > FastMath.abs(e[k])); > + double sp = singularValues[p - 1] / scale; > + double spm1 = singularValues[p - 2] / scale; > + double epm1 = e[p - 2] / scale; > + double sk = singularValues[k] / scale; > + double ek = e[k] / scale; > + double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / > 2.0; > + double c = (sp * epm1) * (sp * epm1); > + double shift = 0.0; > + if ((b != 0.0) | (c != 0.0)) { > + shift = FastMath.sqrt(b * b + c); > + if (b < 0.0) { > + shift = -shift; > + } > + shift = c / (b + shift); > + } > + double f = (sk + sp) * (sk - sp) + shift; > + double g = sk * ek; > + // Chase zeros. > + for (int j = k; j < p - 1; j++) { > + double t = FastMath.hypot(f, g); > + double cs = f / t; > + double sn = g / t; > + if (j != k) { > + e[j - 1] = t; > + } > + f = cs * singularValues[j] + sn * e[j]; > + e[j] = cs * e[j] - sn * singularValues[j]; > + g = sn * singularValues[j + 1]; > + singularValues[j + 1] = cs * singularValues[j + 1]; > + if (wantv) { > + for (int i = 0; i < n; i++) { > + t = cs * V[i][j] + sn * V[i][j + 1]; > + V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + > 1]; > + V[i][j] = t; > + } > + } > + t = FastMath.hypot(f, g); > + cs = f / t; > + sn = g / t; > + singularValues[j] = t; > + f = cs * e[j] + sn * singularValues[j + 1]; > + singularValues[j + 1] = -sn * e[j] + cs * > singularValues[j + 1]; > + g = sn * e[j + 1]; > + e[j + 1] = cs * e[j + 1]; > + if (wantu && (j < m - 1)) { > + for (int i = 0; i < m; i++) { > + t = cs * U[i][j] + sn * U[i][j + 1]; > + U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + > 1]; > + U[i][j] = t; > + } > + } > + } > + e[p - 2] = f; > + iter = iter + 1; > + } > + break; > + // Convergence. > + default: { > + // Make the singular values positive. > + if (singularValues[k] <= 0.0) { > + singularValues[k] = singularValues[k] < 0.0 ? > -singularValues[k] : 0.0; > + if (wantv) { > + for (int i = 0; i <= pp; i++) { > + V[i][k] = -V[i][k]; > + } > + } > + } > + // Order the singular values. > + while (k < pp) { > + if (singularValues[k] >= singularValues[k + 1]) { > + break; > + } > + double t = singularValues[k]; > + singularValues[k] = singularValues[k + 1]; > + singularValues[k + 1] = t; > + if (wantv && (k < n - 1)) { > + for (int i = 0; i < n; i++) { > + t = V[i][k + 1]; > + V[i][k + 1] = V[i][k]; > + V[i][k] = t; > + } > + } > + if (wantu && (k < m - 1)) { > + for (int i = 0; i < m; i++) { > + t = U[i][k + 1]; > + U[i][k + 1] = U[i][k]; > + U[i][k] = t; > + } > + } > + k++; > + } > + iter = 0; > + p--; > + } > + break; > } > } > > - double[][] matAAT = new double[m][m]; > - // > - // create A*A^T > - // > - for (int i = 0; i < m; i++) { > - for (int j = i; j < m; j++) { > - matAAT[i][j] = 0.0; > - for (int k = 0; k < n; k++) { > - matAAT[i][j] += localcopy[i][k] * localcopy[j][k]; > - } > - matAAT[j][i] = matAAT[i][j]; > - } > - } > - int p; > - if (m >= n) { > - p = n; > - // compute eigen decomposition of A^T*A > - eigenDecomposition > - = new EigenDecompositionImpl(new > Array2DRowRealMatrix(matATA), 1); > - singularValues = eigenDecomposition.getRealEigenvalues(); > - cachedV = eigenDecomposition.getV(); > - // compute eigen decomposition of A*A^T > - eigenDecomposition > - = new EigenDecompositionImpl(new > Array2DRowRealMatrix(matAAT), 1); > - cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p > - 1); > + if (!transposed) { > + cachedU = MatrixUtils.createRealMatrix(U); > + cachedV = MatrixUtils.createRealMatrix(V); > } else { > - p = m; > - // compute eigen decomposition of A*A^T > - eigenDecomposition > - = new EigenDecompositionImpl(new > Array2DRowRealMatrix(matAAT), 1); > - singularValues = eigenDecomposition.getRealEigenvalues(); > - cachedU = eigenDecomposition.getV(); > - > - // compute eigen decomposition of A^T*A > - eigenDecomposition > - = new EigenDecompositionImpl(new > Array2DRowRealMatrix(matATA), 1); > - cachedV = eigenDecomposition.getV().getSubMatrix(0, n - 1 , 0, p > - 1); > - } > - for (int i = 0; i < p; i++) { > - singularValues[i] = > FastMath.sqrt(FastMath.abs(singularValues[i])); > - } > - // Up to this point, U and V are computed independently of each > other. > - // There still a sign indetermination of each column of, say, U. > - // The sign is set such that A.V_i=sigma_i.U_i (i<=p) > - // The right sign corresponds to a positive dot product of A.V_i and > U_i > - for (int i = 0; i < p; i++) { > - RealVector tmp = cachedU.getColumnVector(i); > - double > product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp); > - if (product < 0) { > - cachedU.setColumnVector(i, tmp.mapMultiply(-1)); > - } > + cachedU = MatrixUtils.createRealMatrix(V); > + cachedV = MatrixUtils.createRealMatrix(U); > + > } > } > > @@ -193,7 +524,7 @@ public class SingularValueDecompositionI > > if (dimension == 0) { > throw new > NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE, > - minSingularValue, > singularValues[0], true); > + minSingularValue, singularValues[0], true); > } > > final double[][] data = new double[dimension][p]; > @@ -217,19 +548,19 @@ public class SingularValueDecompositionI > > /** {@inheritDoc} */ > public double getConditionNumber() { > - return singularValues[0] / singularValues[singularValues.length - 1]; > + return singularValues[0] / singularValues[FastMath.min(m, n) - 1]; > } > > /** {@inheritDoc} */ > public int getRank() { > - final double threshold = FastMath.max(m, n) * > FastMath.ulp(singularValues[0]); > - > - for (int i = singularValues.length - 1; i >= 0; --i) { > - if (singularValues[i] > threshold) { > - return i + 1; > + double tol = FastMath.max(m, n) * singularValues[0] * EPS; > + int r = 0; > + for (int i = 0; i < singularValues.length; i++) { > + if (singularValues[i] > tol) { > + r++; > } > } > - return 0; > + return r; > } > > /** {@inheritDoc} */ > @@ -253,14 +584,14 @@ public class SingularValueDecompositionI > * @param nonSingular Singularity indicator. > */ > private Solver(final double[] singularValues, final RealMatrix uT, > - final RealMatrix v, final boolean nonSingular) { > + final RealMatrix v, final boolean nonSingular) { > double[][] suT = uT.getData(); > for (int i = 0; i < singularValues.length; ++i) { > final double a; > if (singularValues[i] > 0) { > - a = 1 / singularValues[i]; > + a = 1 / singularValues[i]; > } else { > - a = 0; > + a = 0; > } > final double[] suTi = suT[i]; > for (int j = 0; j < suTi.length; ++j) { > > Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml > URL: > http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1148714&r1=1148713&r2=1148714&view=diff > ============================================================================== > --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) > +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Wed Jul 20 12:15:00 > 2011 > @@ -52,6 +52,9 @@ The <action> type attribute can be add,u > If the output is not quite correct, check for invisible trailing spaces! > --> > <release version="3.0" date="TBD" description="TBD"> > + <action dev="luc" type="fix" > issue="MATH-327,MATH-383,MATH-465,MATH-583,MATH-611" due-to="Christopher Nix"> > + Rewritten SVD decomposition based on JAMA code. > + </action> > <action dev="psteitz" type="fix" issue="MATH-618" due-to="Arne Plose"> > Complex add javadoc says that if either addend has NaN parts, the > result > should be Complex.NaN. Prior to the fix for this issue, NaNs were > propagated > > Modified: > commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java > URL: > http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java?rev=1148714&r1=1148713&r2=1148714&view=diff > ============================================================================== > --- > commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java > (original) > +++ > commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java > Wed Jul 20 12:15:00 2011 > @@ -17,6 +17,10 @@ > > package org.apache.commons.math.linear; > > +import java.io.BufferedReader; > +import java.io.DataInputStream; > +import java.io.IOException; > +import java.io.InputStreamReader; > import java.util.Random; > > import org.junit.Assert; > @@ -163,7 +167,8 @@ public class SingularValueDecompositionI > } > > /** test matrices values */ > - @Test > + // This test is useless since whereas the columns of U and V are linked > + // together, the actual triplet (U,S,V) is not uniquely defined. > public void testMatricesValues1() { > SingularValueDecomposition svd = > new > SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)); > @@ -234,6 +239,56 @@ public class SingularValueDecompositionI > > } > > + /** test MATH-465 */ > + @Test > + public void testRank() { > + double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } }; > + RealMatrix m = new Array2DRowRealMatrix(d); > + SingularValueDecomposition svd = new > SingularValueDecompositionImpl(m); > + Assert.assertEquals(2, svd.getRank()); > + } > + > + /** test MATH-583 */ > + @Test > + public void testStability1() { > + RealMatrix m = new Array2DRowRealMatrix(201, 201); > + loadRealMatrix(m,"matrix1.csv"); > + try { > + new SingularValueDecompositionImpl(m); > + } catch (Exception e) { > + Assert.fail("Exception whilst constructing SVD"); > + } > + } > + > + /** test MATH-327 */ > + @Test > + public void testStability2() { > + RealMatrix m = new Array2DRowRealMatrix(7, 168); > + loadRealMatrix(m,"matrix2.csv"); > + try { > + new SingularValueDecompositionImpl(m); > + } catch (Throwable e) { > + Assert.fail("Exception whilst constructing SVD"); > + } > + } > + > + private void loadRealMatrix(RealMatrix m, String resourceName) { > + try { > + DataInputStream in = new > DataInputStream(getClass().getResourceAsStream(resourceName)); > + BufferedReader br = new BufferedReader(new > InputStreamReader(in)); > + String strLine; > + int row = 0; > + while ((strLine = br.readLine()) != null) { > + int col = 0; > + for (String entry : strLine.split(",")) { > + m.setEntry(row, col++, Double.parseDouble(entry)); > + } > + row++; > + } > + in.close(); > + } catch (IOException e) {} > + } > + > /** test condition number */ > @Test > public void testConditionNumber() { > > > --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
