sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1528514 - 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 Wed, 02 Oct 2013 15:13:58 GMT
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.
      */



Mime
View raw message