Author: desruisseaux
Date: Thu Oct 3 13:56:09 2013
New Revision: 1528868
URL: http://svn.apache.org/r1528868
Log:
Complete the use 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/Solver.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=1528868&r1=1528867&r2=1528868&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] Thu Oct 3 13:56:09 2013
@@ -258,14 +258,20 @@ class GeneralMatrix extends MatrixSIS {
* Returns all elements of the given matrix followed by the error terms
for extended-precision arithmetic.
* 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>
+ * <p>This method may return a direct reference to the internal array.
<strong>Do not modify.</strong>,
+ * unless the {@code copy} argument is {@code true}.</p>
+ *
+ * @param copy If {@code true}, then the returned array is guaranteed to
be a copy, never the internal array.
*/
- static double[] getExtendedElements(final Matrix matrix, final int numRow,
final int numCol) {
+ static double[] getExtendedElements(final Matrix matrix, final int numRow,
final int numCol, final boolean copy) {
double[] elements;
final int length = numRow * numCol * 2;
if (matrix instanceof GeneralMatrix) {
elements = ((GeneralMatrix) matrix).elements;
if (elements.length == length) {
+ if (copy) {
+ elements = elements.clone();
+ }
return elements; // Internal array already uses extended
precision.
} else {
elements = Arrays.copyOf(elements, length);
@@ -427,8 +433,8 @@ class GeneralMatrix extends MatrixSIS {
* Get the matrix element values, together with the error terms if the
matrix
* use extended precision (double-double arithmetic).
*/
- final double[] eltA = getExtendedElements(A, numRow, nc);
- final double[] eltB = getExtendedElements(B, nc, numCol);
+ final double[] eltA = getExtendedElements(A, numRow, nc, false);
+ final double[] eltB = getExtendedElements(B, nc, numCol, false);
final int errorOffset = numRow * numCol; // Where error terms start.
final int errA = numRow * nc;
final int errB = nc * numCol;
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=1528868&r1=1528867&r2=1528868&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] Thu Oct 3 13:56:09 2013
@@ -102,7 +102,7 @@ final class Solver implements Matrix {
*/
static MatrixSIS inverse(final MatrixSIS X) throws
NoninvertibleMatrixException {
final int size = X.getNumRow();
- return solve(X, IDENTITY, size, size);
+ return solve(X, IDENTITY, null, size, size);
}
/**
@@ -123,6 +123,7 @@ final class Solver implements Matrix {
* }
*
* @param X The matrix to invert.
+ * @param Y The desired result of {@code X} × <var>U</var>.
* @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.
@@ -130,11 +131,30 @@ final class Solver implements Matrix {
static MatrixSIS solve(final MatrixSIS X, final Matrix Y, final int size,
final int innerSize)
throws NoninvertibleMatrixException
{
- /*
- * Use a "left-looking", dot-product, Crout/Doolittle algorithm.
- */
+ double[] eltY = null;
+ if (Y instanceof GeneralMatrix) {
+ eltY = ((GeneralMatrix) Y).elements;
+ if (eltY.length == size * innerSize) {
+ eltY = null; // Matrix does not contains error terms.
+ }
+ }
+ return solve(X, Y, eltY, size, innerSize);
+ }
+
+ /**
+ * Implementation of {@code solve} and {@code inverse} methods.
+ * Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+ *
+ * @param eltY Elements and error terms of the {@code Y} matrix, or {@code
null} if not available.
+ */
+ private static MatrixSIS solve(final MatrixSIS X, final Matrix Y, final
double[] eltY,
+ final int size, final int innerSize) throws
NoninvertibleMatrixException
+ {
+ assert (X.getNumRow() == size && X.getNumCol() == size) : size;
+ assert (Y.getNumRow() == size && Y.getNumCol() == innerSize) || (Y
instanceof Solver);
+
final int errorLU = size * size;
- final double[] LU = GeneralMatrix.getExtendedElements(X, size, size);
+ final double[] LU = GeneralMatrix.getExtendedElements(X, size, size,
true);
assert errorLU == GeneralMatrix.indexOfErrors(size, size, LU);
final int[] pivot = new int[size];
for (int j=0; j<size; j++) {
@@ -153,7 +173,14 @@ final class Solver implements Matrix {
column[j + size] = LU[k + errorLU]; // Error
}
/*
- * Apply previous transformations.
+ * Apply previous transformations. This part is equivalent to the
following code,
+ * but using double-double arithmetic instead than the primitive
'double' type:
+ *
+ * double sum = 0;
+ * for (int k=0; k<kmax; k++) {
+ * sum += LU[rowOffset + k] * column[k];
+ * }
+ * LU[rowOffset + i] = (column[j] -= sum);
*/
for (int j=0; j<size; j++) {
final int rowOffset = j*size;
@@ -170,7 +197,8 @@ final class Solver implements Matrix {
sum.storeTo(LU, rowOffset + i, errorLU);
}
/*
- * Find pivot and exchange if necessary.
+ * Find pivot and exchange if necessary. There is no
floating-point arithmetic here
+ * (ignoring the comparison for magnitude order), only work on
index values.
*/
int p = i;
for (int j=i; ++j < size;) {
@@ -187,15 +215,23 @@ final class Solver implements Matrix {
ArraysExt.swap(pivot, p, i);
}
/*
- * Compute multipliers.
+ * Compute multipliers. This part is equivalent to the following
code, but
+ * using double-double arithmetic instead than the primitive
'double' type:
+ *
+ * final double sum = LU[i*size + i];
+ * if (sum != 0.0) {
+ * for (int j=i; ++j < size;) {
+ * LU[j*size + i] /= sum;
+ * }
+ * }
*/
sum.setFrom(LU, i*size + i, errorLU);
if (!sum.isZero()) {
for (int j=i; ++j < size;) {
final int t = j*size + i;
- tmp.setFrom(LU, t, errorLU);
- tmp.divide(sum);
- tmp.storeTo(LU, t, errorLU);
+ tmp.setFrom(sum);
+ tmp.inverseDivide(LU, t, errorLU);
+ tmp.storeTo (LU, t, errorLU);
}
}
}
@@ -210,44 +246,77 @@ final class Solver implements Matrix {
}
}
/*
- * Copy right hand side with pivoting.
- * We will write the result of this method directly in the elements
array.
+ * Copy right hand side with pivoting. Write the result directly in
the elements array
+ * of the result matrix. This block does not perform floating-point
arithmetic operations.
*/
final GeneralMatrix result =
GeneralMatrix.createExtendedPrecision(size, innerSize);
final double[] elements = result.elements;
+ final int errorOffset = size * innerSize;
for (int k=0,j=0; j<size; j++) {
final int p = pivot[j];
for (int i=0; i<innerSize; i++) {
- elements[k++] = Y.getElement(p, i);
+ if (eltY != null) {
+ final int t = p*innerSize + i;
+ elements[k] = eltY[t];
+ elements[k + errorOffset] = eltY[t + errorOffset];
+ } else {
+ elements[k] = Y.getElement(p, i);
+ }
+ k++;
}
}
/*
- * Solve L*Y = B(pivot, :)
+ * Solve L*Y = B(pivot, :). The inner block is equivalent to the
following line,
+ * but using double-double arithmetic instead of 'double' primitive
type:
+ *
+ * elements[loRowOffset + i] -= (elements[rowOffset + i] *
LU[luRowOffset + k]);
*/
for (int k=0; k<size; k++) {
final int rowOffset = k*innerSize; // Offset of row
computed by current iteration.
for (int j=k; ++j < size;) {
- final int loRowOffset = j*innerSize; // Offset of a row
after (locate lower) the current row.
+ final int loRowOffset = j*innerSize; // Offset of some row
after the current row.
final int luRowOffset = j*size; // Offset of the
corresponding row in the LU matrix.
for (int i=0; i<innerSize; i++) {
- elements[loRowOffset + i] -= (elements[rowOffset + i] *
LU[luRowOffset + k]);
+ sum.setFrom (elements, loRowOffset + i, errorOffset);
+ tmp.setFrom (elements, rowOffset + i, errorOffset);
+ tmp.multiply(LU, luRowOffset + k, errorLU);
+ sum.subtract(tmp);
+ sum.storeTo (elements, loRowOffset + i, errorOffset);
}
}
}
/*
- * Solve U*X = Y
+ * Solve U*X = Y. The content of the loop is equivalent to the
following line,
+ * but using double-double arithmetic instead of 'double' primitive
type:
+ *
+ * double sum = LU[k*size + k];
+ * for (int i=0; i<innerSize; i++) {
+ * elements[rowOffset + i] /= sum;
+ * }
+ * for (int j=0; j<k; j++) {
+ * sum = LU[j*size + k];
+ * for (int i=0; i<innerSize; i++) {
+ * elements[upRowOffset + i] -= (elements[rowOffset + i] *
sum);
+ * }
+ * }
*/
for (int k=size; --k >= 0;) {
final int rowOffset = k*innerSize; // Offset of row
computed by current iteration.
- final double d = LU[k*size + k]; // A diagonal element on the
current row.
+ sum.setFrom(LU, k*size + k, errorLU); // A diagonal element
on the current row.
for (int i=0; i<innerSize; i++) { // Apply to all
columns in the current row.
- elements[rowOffset + i] /= d;
+ tmp.setFrom(sum);
+ tmp.inverseDivide(elements, rowOffset + i, errorOffset);
+ tmp.storeTo (elements, rowOffset + i, errorOffset);
}
for (int j=0; j<k; j++) {
final int upRowOffset = j*innerSize; // Offset of a row
before (locate upper) the current row.
- final double c = LU[j*size + k]; // Same column than the
diagonal element, but in the upper row.
+ sum.setFrom(LU, j*size + k, errorLU); // Same column than
the diagonal element, but in the upper row.
for (int i=0; i<innerSize; i++) { // Apply to all
columns in the upper row.
- elements[upRowOffset + i] -= (elements[rowOffset + i] * c);
+ tmp.setFrom(elements, rowOffset + i, errorOffset);
+ tmp.multiply(sum);
+ tmp.subtract(elements, upRowOffset + i, errorOffset);
+ tmp.negate();
+ tmp.storeTo(elements, upRowOffset + i, errorOffset);
}
}
}
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=1528868&r1=1528867&r2=1528868&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] Thu Oct 3 13:56:09 2013
@@ -241,6 +241,16 @@ public final class DoubleDouble extends
}
/**
+ * Sets this {@code DoubleDouble} to the same value than the given
instance.
+ *
+ * @param other The instance to copy.
+ */
+ public void setFrom(final DoubleDouble other) {
+ value = other.value;
+ error = other.error;
+ }
+
+ /**
* Sets the {@link #value} and {@link #error} terms to values read from
the given array.
* This is a convenience method for a frequently used operation,
implemented as below:
*
@@ -576,6 +586,36 @@ public final class DoubleDouble extends
}
/**
+ * Divides this {@code DoubleDouble} by an other double-double value
stored in the given array.
+ * This is a convenience method for a frequently used operation,
implemented as below:
+ *
+ * {@preformat java
+ * divide(array[index], array[index + errorOffset]);
+ * }
+ *
+ * @param array The array from which to get the value and error.
+ * @param index Index of the value in the given array.
+ * @param errorOffset Offset to add to {@code index} in order to get the
index of the error in the given array.
+ */
+ public void divide(final double[] array, final int index, final int
errorOffset) {
+ divide(array[index], array[index + errorOffset]);
+ }
+
+ /**
+ * Divides the given double-double value by this {@code DoubleDouble}.
+ * This is a convenience method for:
+ *
+ * {@preformat java
+ * inverseDivide(other.value, other.error);
+ * }
+ *
+ * @param other The other value to add to this value.
+ */
+ public void inverseDivide(final DoubleDouble other) {
+ inverseDivide(other.value, other.error);
+ }
+
+ /**
* Divides the given double-double value by this {@code DoubleDouble}.
* The result is stored in this instance.
*
@@ -592,6 +632,11 @@ public final class DoubleDouble extends
* @param numeratorError The error of the other value to divide by this
{@code DoubleDouble}.
*/
public void inverseDivide(final double numeratorValue, final double
numeratorError) {
+ if (DISABLED) {
+ value = numeratorValue / value;
+ error = 0;
+ return;
+ }
final double denominatorValue = value;
/*
* The 'b * (a.value / b.value)' part in the method javadoc.
@@ -614,6 +659,22 @@ public final class DoubleDouble extends
}
/**
+ * Divides the given double-double value by this {@code DoubleDouble}.
+ * This is a convenience method for a frequently used operation,
implemented as below:
+ *
+ * {@preformat java
+ * inverseDivide(array[index], array[index + errorOffset]);
+ * }
+ *
+ * @param array The array from which to get the value and error.
+ * @param index Index of the value in the given array.
+ * @param errorOffset Offset to add to {@code index} in order to get the
index of the error in the given array.
+ */
+ public void inverseDivide(final double[] array, final int index, final int
errorOffset) {
+ inverseDivide(array[index], array[index + errorOffset]);
+ }
+
+ /**
* Returns a string representation of this number for debugging purpose.
* The returned string does not need to contains all digits that this
{@code DoubleDouble} can handle.
*
