Author: desruisseaux
Date: Wed Oct 2 15:13:58 2013
New Revision: 1528514
URL: http://svn.apache.org/r1528514
Log:
Matrix multiplications now use double-double arithmetic.
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/Matrix1.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/NonSquareMatrix.java
sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -122,11 +122,7 @@ class GeneralMatrix extends MatrixSIS {
this.numRow = (short) numRow;
this.numCol = (short) numCol;
elements = new double[numRow * numCol * precision];
- for (int k=0,j=0; j<numRow; j++) {
- for (int i=0; i<numCol; i++) {
- elements[k++] = matrix.getElement(j, i);
- }
- }
+ getElements(matrix, numRow, numCol, elements);
inferErrors();
}
@@ -140,6 +136,23 @@ class GeneralMatrix extends MatrixSIS {
}
/**
+ * Copies the elements of the given matrix in the given array.
+ * This method ignores the error terms, if any.
+ *
+ * @param matrix The matrix to copy.
+ * @param numRow The number of rows to copy (usually {@code
matrix.getNumRow()}).
+ * @param numCol The number of columns to copy (usually {@code
matrix.getNumCol()}).
+ * @param elements Where to copy the elements.
+ */
+ private static void getElements(final Matrix matrix, final int numRow,
final int numCol, final double[] elements) {
+ for (int k=0,j=0; j<numRow; j++) {
+ for (int i=0; i<numCol; i++) {
+ elements[k++] = matrix.getElement(j, i);
+ }
+ }
+ }
+
+ /**
* Infers all {@link DoubleDouble#error} with a default values inferred
from {@link DoubleDouble#value}.
* For example if a matrix element is exactly 3.141592653589793, there is
good chances that the user's
* intend was to specify the {@link Math#PI} value, in which case this
method will infer that we would
@@ -153,6 +166,16 @@ class GeneralMatrix extends MatrixSIS {
}
/**
+ * Returns the index of the first {@link DoubleDouble#error} value in the
{@link #elements} array,
+ * or 0 if none. This method returns a non-zero value only if the matrix
has been created in extended
+ * precision mode.
+ */
+ static int indexOfErrors(final int numRow, final int numCol, final
double[] elements) {
+ assert elements.length % (numRow * numCol) == 0;
+ return (numRow * numCol) % elements.length; // A % B is for getting 0
without branching if A == B.
+ }
+
+ /**
* Ensures that the given matrix size is valid for this {@code
GeneralMatrix} implementation.
*/
private static void ensureValidSize(final int numRow, final int numCol) {
@@ -214,6 +237,31 @@ class GeneralMatrix extends MatrixSIS {
}
/**
+ * Returns all elements of the given matrix, possibly including the error
terms for extended-precision arithmetic.
+ * If the returned array contains error terms, then the array will have
twice the normal length.
+ * See {@link #elements} for more discussion.
+ *
+ * <p>This method may return a direct reference to the internal array.
<strong>Do not modify.</strong></p>
+ */
+ private static double[] getExtendedElements(final Matrix matrix, final int
numRow, final int numCol) {
+ if (matrix instanceof MatrixSIS) {
+ return ((MatrixSIS) matrix).getExtendedElements();
+ }
+ final double[] elements = new double[numRow * numCol];
+ getElements(matrix, numRow, numCol, elements);
+ return elements;
+ }
+
+ /**
+ * Returns a direct reference to the internal array, which may or may not
contains error values
+ * for extended precision arithmetic.
+ */
+ @Override
+ final double[] getExtendedElements() {
+ return elements;
+ }
+
+ /**
* {@inheritDoc}
*/
@Override
@@ -251,12 +299,7 @@ class GeneralMatrix extends MatrixSIS {
* At this point, the 'double' values are those of an affine
transform.
* If this matrix uses extended precision, ensures that their
errors are zero.
*/
- for (i = base + numRow*numCol; i < elements.length; i++) {
- if (elements[i] != 0) {
- return false;
- }
- }
- return true;
+ return errorsAreZero(base + numRow*numCol);
}
}
return false;
@@ -287,8 +330,18 @@ class GeneralMatrix extends MatrixSIS {
* At this point, the 'double' values are those of an identity
transform.
* If this matrix uses extended precision, ensures that all errors are
zero.
*/
- for (int i=length; i<elements.length; i++) {
- if (elements[i] != 0) {
+ return errorsAreZero(length);
+ }
+
+ /**
+ * Returns {@code true} if this matrix has no error elements, or if all
error elements starting
+ * at the given index are zero.
+ *
+ * @param i Index of the first error elements to check (may be greater
than the array length).
+ */
+ private boolean errorsAreZero(int i) {
+ while (i < elements.length) {
+ if (elements[i++] != 0) {
return false;
}
}
@@ -305,7 +358,7 @@ class GeneralMatrix extends MatrixSIS {
public void transpose() {
final int numRow = this.numRow; // Protection against accidental
changes.
final int numCol = this.numCol;
- final int errors = (numRow * numCol) % elements.length; // Where error
values start, or 0 if none.
+ final int errors = indexOfErrors(numRow, numCol, elements); // Where
error values start, or 0 if none.
for (int j=0; j<numRow; j++) {
for (int i=0; i<j; i++) {
final int lo = j*numCol + i;
@@ -358,26 +411,47 @@ class GeneralMatrix extends MatrixSIS {
}
/**
- * {@inheritDoc}
+ * 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.
*/
- @Override
- public final MatrixSIS multiply(final Matrix matrix) {
+ final void setToProduct(final MatrixSIS A, final Matrix B) {
final int numRow = this.numRow; // Protection against accidental
changes.
final int numCol = this.numCol;
- final int nc = matrix.getNumCol();
- ensureNumRowMatch(numCol, matrix, nc);
- final MatrixSIS result = Matrices.createZero(numRow, nc);
- for (int j=0; j<numRow; j++) {
- final int srcOff = j * numCol;
- for (int i=0; i<nc; i++) {
- double sum = 0;
- for (int k=0; k<numCol; k++) {
- sum += elements[srcOff + k] * matrix.getElement(k, i);
+ final int nc = A.getNumCol();
+ assert B.getNumRow() == nc;
+ assert numRow == A.getNumRow() && numCol == B.getNumCol();
+ /*
+ * Get the matrix element values, together with the error terms if the
matrix
+ * use extended precision (double-double arithmetic).
+ */
+ final double[] eltA = A.getExtendedElements();
+ final double[] eltB = getExtendedElements(B, nc, numCol);
+ final int errors = indexOfErrors(numRow, numCol, elements); //
Where error values start, or 0 if none.
+ final int errA = indexOfErrors(numRow, nc, eltA);
+ final int errB = indexOfErrors(nc, numCol, eltB);
+ /*
+ * Compute the product, to be stored directly in 'this'.
+ */
+ final DoubleDouble dot = new DoubleDouble();
+ final DoubleDouble sum = new DoubleDouble();
+ for (int k=0,j=0; j<numRow; j++) {
+ for (int i=0; i<numCol; i++) {
+ sum.clear();
+ int iB = i; // Index of values in a single column of B.
+ int iA = j * nc; // Index of values in a single row of A.
+ final int nextRow = iA + nc;
+ while (iA < nextRow) {
+ dot.value = eltA[iA];
+ dot.error = (errA != 0) ? eltA[iA + errA] : 0;
+ dot.multiply(eltB[iB], (errB != 0) ? eltB[iB + errB] : 0);
+ sum.add(dot);
+ iB += numCol; // Move to next row of B.
+ iA++; // Move to next column of A.
}
- result.setElement(j, i, sum);
+ elements[k + errors] = sum.error;
+ elements[k++] = sum.value;
}
}
- return result;
}
/**
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java?rev=1528514&r1=1528513&r2=1528514&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
[UTF-8] Wed Oct 2 15:13:58 2013
@@ -269,19 +269,6 @@ public final class Matrix1 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS multiply(final Matrix matrix) {
- final int nc = matrix.getNumCol();
- ensureNumRowMatch(SIZE, matrix, nc);
- if (nc != SIZE) {
- return new NonSquareMatrix(this, matrix, 1);
- }
- return new Matrix1(m00 * matrix.getElement(0,0));
- }
-
- /**
* Returns {@code true} if the specified object is of type {@code Matrix1}
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/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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -286,23 +286,6 @@ public final class Matrix2 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS multiply(final Matrix matrix) {
- final int nc = matrix.getNumCol();
- ensureNumRowMatch(SIZE, matrix, nc);
- if (nc != SIZE) {
- return new NonSquareMatrix(this, matrix, 1);
- }
- final Matrix2 k = (matrix instanceof Matrix2) ? (Matrix2) matrix : new
Matrix2(matrix);
- return new Matrix2(m00 * k.m00 + m01 * k.m10,
- m00 * k.m01 + m01 * k.m11,
- m10 * k.m00 + m11 * k.m10,
- m10 * k.m01 + m11 * k.m11);
- }
-
- /**
* 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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -323,28 +323,6 @@ public final class Matrix3 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS multiply(final Matrix matrix) {
- final int nc = matrix.getNumCol();
- ensureNumRowMatch(SIZE, matrix, nc);
- if (nc != SIZE) {
- return new NonSquareMatrix(this, matrix, 1);
- }
- final Matrix3 k = (matrix instanceof Matrix3) ? (Matrix3) matrix : new
Matrix3(matrix);
- return new Matrix3(m00 * k.m00 + m01 * k.m10 + m02 * k.m20,
- m00 * k.m01 + m01 * k.m11 + m02 * k.m21,
- m00 * k.m02 + m01 * k.m12 + m02 * k.m22,
- m10 * k.m00 + m11 * k.m10 + m12 * k.m20,
- m10 * k.m01 + m11 * k.m11 + m12 * k.m21,
- m10 * k.m02 + m11 * k.m12 + m12 * k.m22,
- m20 * k.m00 + m21 * k.m10 + m22 * k.m20,
- m20 * k.m01 + m21 * k.m11 + m22 * k.m21,
- m20 * k.m02 + m21 * k.m12 + m22 * k.m22);
- }
-
- /**
* 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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -374,35 +374,6 @@ public final class Matrix4 extends Matri
}
/**
- * {@inheritDoc}
- */
- @Override
- public MatrixSIS multiply(final Matrix matrix) {
- final int nc = matrix.getNumCol();
- ensureNumRowMatch(SIZE, matrix, nc);
- if (nc != SIZE) {
- return new NonSquareMatrix(this, matrix, 1);
- }
- final Matrix4 k = (matrix instanceof Matrix4) ? (Matrix4) matrix : new
Matrix4(matrix);
- return new Matrix4(m00 * k.m00 + m01 * k.m10 + m02 * k.m20 + m03
* k.m30,
- m00 * k.m01 + m01 * k.m11 + m02 * k.m21 + m03
* k.m31,
- m00 * k.m02 + m01 * k.m12 + m02 * k.m22 + m03
* k.m32,
- m00 * k.m03 + m01 * k.m13 + m02 * k.m23 + m03
* k.m33,
- m10 * k.m00 + m11 * k.m10 + m12 * k.m20 + m13
* k.m30,
- m10 * k.m01 + m11 * k.m11 + m12 * k.m21 + m13
* k.m31,
- m10 * k.m02 + m11 * k.m12 + m12 * k.m22 + m13
* k.m32,
- m10 * k.m03 + m11 * k.m13 + m12 * k.m23 + m13
* k.m33,
- m20 * k.m00 + m21 * k.m10 + m22 * k.m20 + m23
* k.m30,
- m20 * k.m01 + m21 * k.m11 + m22 * k.m21 + m23
* k.m31,
- m20 * k.m02 + m21 * k.m12 + m22 * k.m22 + m23
* k.m32,
- m20 * k.m03 + m21 * k.m13 + m22 * k.m23 + m23
* k.m33,
- m30 * k.m00 + m31 * k.m10 + m32 * k.m20 + m33
* k.m30,
- m30 * k.m01 + m31 * k.m11 + m32 * k.m21 + m33
* k.m31,
- m30 * k.m02 + m31 * k.m12 + m32 * k.m22 + m33
* k.m32,
- m30 * k.m03 + m31 * k.m13 + m32 * k.m23 + m33
* k.m33);
- }
-
- /**
* 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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -126,6 +126,17 @@ public abstract class MatrixSIS implemen
}
/**
+ * Returns all matrix elements, possibly including the error terms for
extended-precision arithmetic.
+ * If the returned array contains error terms, then the array will have
twice the normal length.
+ * See {@link GeneralMatrix#elements} for more discussion.
+ *
+ * <p>This method may return a direct reference to the internal array.
<strong>Do not modify.</strong></p>
+ */
+ double[] getExtendedElements() {
+ return getElements();
+ }
+
+ /**
* Returns a copy of all matrix elements in a flat, row-major (column
indices vary fastest) array.
* The array length is <code>{@linkplain #getNumRow()} * {@linkplain
#getNumCol()}</code>.
*
@@ -204,7 +215,20 @@ public abstract class MatrixSIS implemen
* @throws MismatchedMatrixSizeException if the number of rows in the
given matrix is not equals to the
* number of columns in this matrix.
*/
- public abstract MatrixSIS multiply(Matrix matrix) throws
MismatchedMatrixSizeException;
+ public MatrixSIS multiply(final Matrix matrix) throws
MismatchedMatrixSizeException {
+ final int numRow = getNumRow();
+ final int numCol = getNumCol();
+ final int nc = matrix.getNumCol();
+ ensureNumRowMatch(numCol, matrix, nc);
+ final GeneralMatrix result;
+ if (numRow == nc) {
+ result = new GeneralMatrix(numRow, nc, false, 2);
+ } else {
+ result = new NonSquareMatrix(numRow, nc, false, 2);
+ }
+ result.setToProduct(this, matrix);
+ return result;
+ }
/**
* Returns the value of <var>U</var> which solves {@code this} Ã
<var>U</var> = {@code matrix}.
Modified:
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java?rev=1528514&r1=1528513&r2=1528514&view=diff
==============================================================================
---
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
[UTF-8] (original)
+++
sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
[UTF-8] Wed Oct 2 15:13:58 2013
@@ -82,35 +82,13 @@ final class NonSquareMatrix extends Gene
}
/**
- * Initializes this matrix to the product of the given matrices.
- * This constructor shall be invoked only when the result is known to be a
non-square matrix.
- */
- NonSquareMatrix(final Matrix A, final Matrix B, final int precision) {
- super(A.getNumRow(), B.getNumCol(), false, precision);
- final int numRow = this.numRow; // Protection against accidental
changes.
- final int numCol = this.numCol;
- final int common = A.getNumCol();
- ensureNumRowMatch(common, B, numCol);
- int offset = 0;
- for (int j=0; j<numRow; j++) {
- for (int i=0; i<numCol; i++) {
- double sum = 0;
- for (int k=0; k<common; k++) {
- sum += A.getElement(j, k) * B.getElement(k, i);
- }
- elements[offset++] = sum;
- }
- }
- }
-
- /**
* Sets the value of this matrix to its transpose.
*/
@Override
public void transpose() {
final short numRow = this.numRow; // Protection against accidental
changes before we are done.
final short numCol = this.numCol;
- final int errors = (numRow * numCol) % elements.length; // Where
error values start, or 0 if none.
+ final int errors = indexOfErrors(numRow, numCol, elements); // Where
error values start, or 0 if none.
final double[] copy = elements.clone();
int k = 0;
for (int j=0; j<numRow; j++) {
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=1528514&r1=1528513&r2=1528514&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 15:13:58 2013
@@ -206,6 +206,14 @@ public final class DoubleDouble extends
}
/**
+ * Resets the {@link #value} and {@link #error} terms to zero.
+ */
+ public void clear() {
+ value = 0;
+ error = 0;
+ }
+
+ /**
* Equivalent to a call to {@code setToQuickSum(value, error)} inlined.
* This is invoked after addition or multiplication operations.
*/
