Author: desruisseaux
Date: Wed Oct 2 20:10:31 2013
New Revision: 1528611
URL: http://svn.apache.org/r1528611
Log:
Added a STRICTFP static final flag for verification purpose only.
First piece of double-double arithmetic in matrix inversion.
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -400,22 +400,6 @@ class GeneralMatrix extends MatrixSIS {
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS inverse() throws NoninvertibleMatrixException {
- return Solver.inverse(this);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS solve(final Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException {
- throw new UnsupportedOperationException(); // TODO
- }
-
- /**
* Sets this matrix to the product of the given matrices: {@code this = A
× B}.
* The matrix sizes much match - this is not verified unless assertions
are enabled.
*/
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -265,27 +265,6 @@ public final class Matrix2 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS inverse() throws NoninvertibleMatrixException {
- final double det = m00*m11 - m01*m10;
- if (det == 0) {
- throw new NoninvertibleMatrixException();
- }
- return new Matrix2(m11 / det, -m01 / det,
- -m10 / det, m00 / det);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS solve(final Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException {
- throw new UnsupportedOperationException(); // TODO
- }
-
- /**
* Returns {@code true} if the specified object is of type {@code Matrix2}
and
* all of the data members are equal to the corresponding data members in
this matrix.
*
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -307,22 +307,6 @@ public final class Matrix3 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS inverse() throws NoninvertibleMatrixException {
- return Solver.inverse(this);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS solve(final Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException {
- throw new UnsupportedOperationException(); // TODO
- }
-
- /**
* Returns {@code true} if the specified object is of type {@code Matrix3}
and
* all of the data members are equal to the corresponding data members in
this matrix.
*
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -358,22 +358,6 @@ public final class Matrix4 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS inverse() throws NoninvertibleMatrixException {
- return Solver.inverse(this);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS solve(final Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException {
- throw new UnsupportedOperationException(); // TODO
- }
-
- /**
* Returns {@code true} if the specified object is of type {@code Matrix4}
and
* all of the data members are equal to the corresponding data members in
this matrix.
*
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
(original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
Wed Oct 2 20:10:31 2013
@@ -230,7 +230,9 @@ public abstract class MatrixSIS implemen
* to the number of columns in this matrix.
* @throws NoninvertibleMatrixException if this matrix is not invertible.
*/
- public abstract MatrixSIS solve(Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException;
+ public MatrixSIS solve(final Matrix matrix) throws
MismatchedMatrixSizeException, NoninvertibleMatrixException {
+ throw new UnsupportedOperationException(); // TODO
+ }
/**
* Returns the inverse of this matrix.
@@ -240,7 +242,9 @@ public abstract class MatrixSIS implemen
*
* @see java.awt.geom.AffineTransform#createInverse()
*/
- public abstract MatrixSIS inverse() throws NoninvertibleMatrixException;
+ public MatrixSIS inverse() throws NoninvertibleMatrixException {
+ return Solver.inverse(this);
+ }
/**
* Compares the given matrices for equality, using the given absolute
tolerance threshold.
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -17,6 +17,7 @@
package org.apache.sis.referencing.operation.matrix;
import org.opengis.referencing.operation.Matrix;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.util.ArraysExt;
@@ -100,7 +101,8 @@ final class Solver implements Matrix {
* (this is <strong>not</strong> verified by this method).
*/
static MatrixSIS inverse(final MatrixSIS X) throws
NoninvertibleMatrixException {
- return solve(X, IDENTITY, X.getNumRow());
+ final int size = X.getNumRow();
+ return solve(X, IDENTITY, size, size);
}
/**
@@ -121,28 +123,34 @@ final class Solver implements Matrix {
* }
*
* @param X The matrix to invert.
+ * @param size The value of {@code X.getNumRow()}, {@code X.getNumCol()}
and {@code Y.getNumRow()}.
* @param innerSize The value of {@code Y.getNumCol()}.
* @throws NoninvertibleMatrixException If the {@code X} matrix is
singular.
*/
- static MatrixSIS solve(final MatrixSIS X, final Matrix Y, final int
innerSize)
+ static MatrixSIS solve(final MatrixSIS X, final Matrix Y, final int size,
final int innerSize)
throws NoninvertibleMatrixException
{
- final int size = X.getNumRow();
/*
* Use a "left-looking", dot-product, Crout/Doolittle algorithm.
*/
- final double[] LU = X.getElements();
+ final int errorLU = size * size;
+ final double[] LU = GeneralMatrix.getExtendedElements(X, size, size);
+ assert errorLU == GeneralMatrix.indexOfErrors(size, size, LU);
final int[] pivot = new int[size];
for (int j=0; j<size; j++) {
pivot[j] = j;
}
- final double[] column = new double[size];
+ final double[] column = new double[size * 2];
+ final DoubleDouble dot = new DoubleDouble();
+ final DoubleDouble sum = new DoubleDouble();
for (int i=0; i<size; i++) {
/*
* Make a copy of the i-th column.
*/
for (int j=0; j<size; j++) {
- column[j] = LU[j*size + i];
+ final int k = j*size + i;
+ column[j] = LU[k];
+ column[j + size] = LU[k + errorLU];
}
/*
* Apply previous transformations.
@@ -150,11 +158,16 @@ final class Solver implements Matrix {
for (int j=0; j<size; j++) {
final int rowOffset = j*size;
final int kmax = Math.min(j,i);
- double s = 0.0;
+ sum.clear();
for (int k=0; k<kmax; k++) {
- s += LU[rowOffset + k] * column[k];
- }
- LU[rowOffset + i] = (column[j] -= s);
+ dot.value = LU[rowOffset + k];
+ dot.error = LU[rowOffset + k + errorLU];
+ dot.multiply(column[k], column[k + size]);
+ sum.add(dot);
+ }
+ sum.add(-column[j], -column[j + size]);
+ LU[rowOffset + i] = column[j] = -sum.value;
+ LU[rowOffset + i + errorLU] = column[j + size] = -sum.error;
}
/*
* Find pivot and exchange if necessary.
Modified:
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -18,6 +18,7 @@ package org.apache.sis.referencing.opera
import java.util.Random;
import Jama.Matrix;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.test.TestCase;
import org.apache.sis.test.TestUtilities;
import org.apache.sis.test.DependsOnMethod;
@@ -58,7 +59,7 @@ public abstract strictfp class MatrixTes
* The matrix elements used in this class varies between 0 and 100,
* and the {@code Math.ulp(100.0)} value is approximatively 1.4E-14.
*/
- static final double TOLERANCE = 1E-10;
+ static final double TOLERANCE = DoubleDouble.STRICTFP ? STRICT : 1E-8;
/**
* Number of random matrices to try in arithmetic operation tests.
@@ -104,9 +105,10 @@ public abstract strictfp class MatrixTes
final int numCol = actual.getNumCol();
assertEquals("numRow", expected.getRowDimension(), numRow);
assertEquals("numCol", expected.getColumnDimension(), numCol);
+ final String name = actual.getClass().getSimpleName();
for (int j=0; j<numRow; j++) {
for (int i=0; i<numCol; i++) {
- assertEquals(expected.get(j,i), actual.getElement(j,i),
tolerance);
+ assertEquals(name, expected.get(j,i), actual.getElement(j,i),
tolerance);
}
}
}
@@ -273,7 +275,7 @@ public abstract strictfp class MatrixTes
m += e*e;
}
m = StrictMath.sqrt(m);
- assertEquals(1, m, TOLERANCE);
+ assertEquals(1, m, 1E-12);
}
}
Modified:
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -80,7 +80,7 @@ public final strictfp class SolverTest e
out.println(e); // "Matrix is singular."
continue;
}
- final MatrixSIS U = Solver.solve(matrix, matrixArg,
matrixArg.getNumCol());
+ final MatrixSIS U = Solver.solve(matrix, matrixArg,
matrixArg.getNumRow(), matrixArg.getNumCol());
MatrixTestCase.assertMatrixEquals(jama, U,
MatrixTestCase.TOLERANCE);
}
}
Modified:
sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -62,6 +62,23 @@ import org.apache.sis.math.MathFunctions
*/
public final class DoubleDouble extends Number {
/**
+ * {@code true} for disabling the extended precision. This variable should
always be {@code false},
+ * except for testing purpose. If set to {@code true}, then all
double-double arithmetic operations
+ * are immediately followed by a clearing of {@link DoubleDouble#error}.
The result should then be
+ * identical to computation performed using the normal {@code double}
arithmetic.
+ *
+ * <p>Since this flag is static final, all expressions of the form {@code
if (STRICTFP)} should be
+ * omitted by the compiler from the class files in normal operations.</p>
+ *
+ * <p>Setting this flag to {@code true} causes some JUnit tests to fail.
This is normal. The main
+ * purpose of this flag is to allow {@link
org.apache.sis.referencing.operation.matrix.MatrixTestCase}
+ * to perform strict comparisons of matrix operation results with JAMA,
which is taken as the reference
+ * implementation. Since JAMA uses {@code double} arithmetic, SIS needs to
disable {@code double-double}
+ * arithmetic if the results are to be compared for strict equality.</p>
+ */
+ public static final boolean STRICTFP = false;
+
+ /**
* For cross-version compatibility.
*/
private static final long serialVersionUID = -7602414219228638550L;
@@ -201,6 +218,7 @@ public final class DoubleDouble extends
* the given value is assumed to be the most accurate available
representation.
*/
public static double errorForWellKnownValue(final double value) {
+ if (STRICTFP) return 0;
final int i = Arrays.binarySearch(VALUES, Math.abs(value));
return (i >= 0) ? MathFunctions.xorSign(ERRORS[i], value) : 0;
}
@@ -219,6 +237,7 @@ public final class DoubleDouble extends
*/
final void normalize() {
error += (value - (value += error));
+ if (STRICTFP) error = 0;
}
/**
@@ -233,6 +252,7 @@ public final class DoubleDouble extends
public void setToQuickSum(final double a, final double b) {
value = a + b;
error = b - (value - a);
+ if (STRICTFP) error = 0;
}
/**
@@ -247,6 +267,7 @@ public final class DoubleDouble extends
value = a + b;
final double v = value - a;
error = (a - (value - v)) + (b - v);
+ if (STRICTFP) error = 0;
}
/**
@@ -267,6 +288,7 @@ public final class DoubleDouble extends
final double bhi = t - (t - b);
final double blo = b - bhi;
error = ((ahi*bhi - value) + ahi*blo + alo*bhi) + alo*blo;
+ if (STRICTFP) error = 0;
}
/**
Modified:
sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java?rev=1528611&r1=1528610&r2=1528611&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] Wed Oct 2 20:10:31 2013
@@ -33,6 +33,7 @@ import static java.lang.StrictMath.*;
/**
* Tests {@link DoubleDouble} using {@link BigDecimal} as the references.
+ * Those tests need {@link DoubleDouble#STRICTFP} to be set to {@code false}.
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.4
