sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1528868 - in /sis/branches/JDK7/core: sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/ sis-utility/src/main/java/org/apache/sis/internal/util/
Date Thu, 03 Oct 2013 13:56:09 GMT
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.
      *



Mime
View raw message