diff --git a/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrix.java b/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrix.java new file mode 100644 index 000000000..3ba8928e5 --- /dev/null +++ b/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrix.java @@ -0,0 +1,2078 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.linear; + +import java.io.Serializable; + +import org.apache.commons.math.MathRuntimeException; + +/** + * Cache-friendly implementation of RealMatrix using recursive array layouts to store + * the matrix elements. + *

+ * As of 2009-02-13, this implementation does not work! The padding at left and bottom + * sides of the matrix should be cleared after some operations like scalerAdd + * and is not. Also there is a limitation in the multiplication that can only + * process matrices with sizes similar enough to have the same power of two + * number of tiles in all three matrices A, B and C such that C = A*B. These + * parts have not been fixed since the performance gain with respect to + * DenseRealMatrix are not very important, and the numerical stability is not + * good. This may well be due to a bad implementation. This code has been put + * in the experimental part for the record, putting it into production would + * require solving all these issues. + *

+ *

+ * This implementation is based on the 2002 paper: Recursive Array Layouts + * and Fast Matrix Multiplication by Siddhartha Chatterjee, Alvin R. Lebeck, + * Praveen K. Patnala and Mithuna Thottethodi. + *

+ *

+ * The matrix is split into several rectangular tiles. The tiles are laid out using + * a space-filling curve in a 2k×2k square. This + * implementation uses the Gray-Morton layout which starts as follows for a three-level + * recursion (i.e. an 8x8 matrix). The tiles size are adjusted in order to have the + * 2k×2k square. This may require padding at the right and + * bottom sides of the matrix (see above paper for a discussion of this padding feature). + *

+ * 
+ *                    |
+ *    00 01 | 06 07   |   24  25 | 30  31
+ *    03 02 | 05 04   |   27  26 | 29  28
+ *    ------+------   |   -------+-------
+ *    12 13 | 10 11   |   20  21 | 18  19
+ *    15 14 | 09 08   |   23  22 | 17  16
+ *                    |
+ * -------------------+--------------------
+ *                    |
+ *    48 49 | 54 55   |   40  41 | 46  47
+ *    51 50 | 53 52   |   43  42 | 45  44
+ *    ------+------   |   -------+-------
+ *    60 61 | 58 59   |   36  37 | 34  35
+ *    63 62 | 57 56   |   39  38 | 33  32
+ *                    |
+ *
+ * @version $Revision$ $Date$ + * @since 2.0 + */ +public class RecursiveLayoutRealMatrix extends AbstractRealMatrix implements Serializable { + + /** Serializable version identifier */ + private static final long serialVersionUID = 1607919006739190004L; + + /** Maximal allowed tile size in bytes. + *

In order to avoid cache miss during multiplication, + * a suggested value is cache_size/3.

+ */ + private static final int MAX_TILE_SIZE_BYTES = (64 * 1024) / 3; + //private static final int MAX_TILE_SIZE_BYTES = 32; + + /** Storage array for matrix elements. */ + private final double data[]; + + /** Number of rows of the matrix. */ + private final int rows; + + /** Number of columns of the matrix. */ + private final int columns; + + /** Number of terminal tiles along rows and columns (guaranteed to be a power of 2). */ + private final int tileNumber; + + /** Number of rows in each terminal tile. */ + private final int tileSizeRows; + + /** Number of columns in each terminal tile. */ + private final int tileSizeColumns; + + /** + * Create a new matrix with the supplied row and column dimensions. + * + * @param rows the number of rows in the new matrix + * @param columns the number of columns in the new matrix + * @throws IllegalArgumentException if row or column dimension is not + * positive + */ + public RecursiveLayoutRealMatrix(final int rows, final int columns) + throws IllegalArgumentException { + + super(rows, columns); + this.rows = rows; + this.columns = columns; + + // compute optimal layout + tileNumber = tilesNumber(rows, columns); + tileSizeRows = tileSize(rows, tileNumber); + tileSizeColumns = tileSize(columns, tileNumber); + + // create storage array + data = new double[tileNumber * tileNumber * tileSizeRows * tileSizeColumns]; + + } + + /** + * Create a new dense matrix copying entries from raw layout data. + *

The input array must be in raw layout.

+ *

Calling this constructor is equivalent to call: + * 

matrix = new RecursiveLayoutRealMatrix(rawData.length, rawData[0].length,
+     *                                             toRecursiveLayout(rawData), false);
+ *

+ * @param rawData data for new matrix, in raw layout + * + * @exception IllegalArgumentException if rawData shape is + * inconsistent with tile layout + * @see #DenseRealMatrix(int, int, double[][], boolean) + */ + public RecursiveLayoutRealMatrix(final double[][] rawData) + throws IllegalArgumentException { + this(rawData.length, rawData[0].length, toRecursiveLayout(rawData), false); + } + + /** + * Create a new dense matrix copying entries from recursive layout data. + *

The input array must already be in recursive layout.

+ * @param rows the number of rows in the new matrix + * @param columns the number of columns in the new matrix + * @param data data for new matrix, in recursive layout + * @param copyArray if true, the input array will be copied, otherwise + * it will be referenced + * + * @exception IllegalArgumentException if data size is + * inconsistent with matrix size + * @see #toRecursiveLayout(double[][]) + * @see #RecursiveLayoutRealMatrix(double[][]) + */ + public RecursiveLayoutRealMatrix(final int rows, final int columns, + final double[] data, final boolean copyArray) + throws IllegalArgumentException { + + super(rows, columns); + this.rows = rows; + this.columns = columns; + + // compute optimal layout + tileNumber = tilesNumber(rows, columns); + tileSizeRows = tileSize(rows, tileNumber); + tileSizeColumns = tileSize(columns, tileNumber); + + // create storage array + final int expectedLength = tileNumber * tileNumber * tileSizeRows * tileSizeColumns; + if (data.length != expectedLength) { + throw MathRuntimeException.createIllegalArgumentException("wrong array size (got {0}, expected {1})", + new Object[] { + data.length, + expectedLength + }); + } + + if (copyArray) { + // allocate storage array + this.data = data.clone(); + } else { + // reference existing array + this.data = data; + } + + } + + /** + * Convert a data array from raw layout to recursive layout. + *

+ * Raw layout is the straightforward layout where element at row i and + * column j is in array element rawData[i][j]. Recursive layout + * is the layout used in {@link RecursiveLayoutRealMatrix} instances, where the matrix + * is stored in a dimension 1 array using a space-filling curve to spread the matrix + * elements along the array. + *

+ * @param rawData data array in raw layout + * @return a new data array containing the same entries but in recursive layout + * @exception IllegalArgumentException if rawData is not rectangular + * (not all rows have the same length) + * @see #RecursiveLayoutRealMatrix(int, int, double[], boolean) + */ + public static double[] toRecursiveLayout(final double[][] rawData) + throws IllegalArgumentException { + + final int rows = rawData.length; + final int columns = rawData[0].length; + + // compute optimal layout + final int tileNumber = tilesNumber(rows, columns); + final int tileSizeRows = tileSize(rows, tileNumber); + final int tileSizeColumns = tileSize(columns, tileNumber); + + // safety checks + for (int i = 0; i < rawData.length; ++i) { + final int length = rawData[i].length; + if (length != columns) { + throw MathRuntimeException.createIllegalArgumentException( + "some rows have length {0} while others have length {1}", + new Object[] { columns, length }); + } + } + + // convert array row after row + final double[] data = new double[tileNumber * tileNumber * tileSizeRows * tileSizeColumns]; + for (int i = 0; i < rawData.length; ++i) { + final int iTile = i / tileSizeRows; + final double[] rawDataI = rawData[i]; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int tileStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns; + final int dataStart = tileStart + (i - iTile * tileSizeRows) * tileSizeColumns; + final int jStart = jTile * tileSizeColumns; + if (jStart < columns) { + final int jEnd = Math.min(jStart + tileSizeColumns, columns); + System.arraycopy(rawDataI, jStart, data, dataStart, jEnd - jStart); + } + } + } + + return data; + + } + + /** {@inheritDoc} */ + public RealMatrix createMatrix(final int rowDimension, final int columnDimension) + throws IllegalArgumentException { + return new RecursiveLayoutRealMatrix(rowDimension, columnDimension); + } + + /** {@inheritDoc} */ + public RealMatrix copy() { + return new RecursiveLayoutRealMatrix(rows, columns, data, true); + } + + /** {@inheritDoc} */ + public RealMatrix add(final RealMatrix m) + throws IllegalArgumentException { + try { + return add((RecursiveLayoutRealMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkAdditionCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // perform addition tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + + // perform addition on the current tile + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + out.data[k] = data[k] + m.getEntry(p, q); + } + } + + } + + return out; + + } + } + + /** + * Compute the sum of this and m. + * + * @param m matrix to be added + * @return this + m + * @throws IllegalArgumentException if m is not the same size as this + */ + public RecursiveLayoutRealMatrix add(final RecursiveLayoutRealMatrix m) + throws IllegalArgumentException { + + // safety check + checkAdditionCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // streamlined addition + for (int i = 0; i < data.length; ++i) { + out.data[i] = data[i] + m.data[i]; + } + + return out; + + } + + /** {@inheritDoc} */ + public RealMatrix subtract(final RealMatrix m) + throws IllegalArgumentException { + try { + return subtract((RecursiveLayoutRealMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkSubtractionCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // perform subtraction tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + + // perform addition on the current tile + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + out.data[k] = data[k] - m.getEntry(p, q); + } + } + + } + + return out; + + } + } + + /** + * Compute this minus m. + * + * @param m matrix to be subtracted + * @return this - m + * @throws IllegalArgumentException if m is not the same size as this + */ + public RecursiveLayoutRealMatrix subtract(final RecursiveLayoutRealMatrix m) + throws IllegalArgumentException { + + // safety check + checkSubtractionCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // streamlined subtraction + for (int i = 0; i < data.length; ++i) { + out.data[i] = data[i] - m.data[i]; + } + + return out; + + } + + /** {@inheritDoc} */ + public RealMatrix scalarAdd(final double d) + throws IllegalArgumentException { + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // streamlined addition + for (int i = 0; i < data.length; ++i) { + out.data[i] = data[i] + d; + } + + return out; + + } + + /** {@inheritDoc} */ + public RealMatrix scalarMultiply(final double d) + throws IllegalArgumentException { + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, columns); + + // streamlined multiplication + for (int i = 0; i < data.length; ++i) { + out.data[i] = data[i] * d; + } + + return out; + + } + + /** {@inheritDoc} */ + public RealMatrix multiply(final RealMatrix m) + throws IllegalArgumentException { + try { + return multiply((RecursiveLayoutRealMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkMultiplicationCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, m.getColumnDimension()); + + // perform multiplication tile-wise, to ensure good cache behavior + for (int index = 0; index < out.tileNumber * out.tileNumber; ++index) { + final int tileStart = index * out.tileSizeRows * out.tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int iStart = iTile * out.tileSizeRows; + final int iEnd = Math.min(iStart + out.tileSizeRows, out.rows); + final int jStart = jTile * out.tileSizeColumns; + final int jEnd = Math.min(jStart + out.tileSizeColumns, out.columns); + + // perform multiplication for current tile + for (int kTile = 0; kTile < tileNumber; ++kTile) { + final int kTileStart = tileIndex(iTile, kTile) * tileSizeRows * tileSizeColumns; + for (int i = iStart, lStart = kTileStart, oStart = tileStart; + i < iEnd; + ++i, lStart += tileSizeColumns, oStart += out.tileSizeColumns) { + final int lEnd = Math.min(lStart + tileSizeColumns, columns); + for (int j = jStart, o = oStart; j < jEnd; ++j, ++o) { + double sum = 0; + for (int l = lStart, k = kTile * tileSizeColumns; l < lEnd; ++l, ++k) { + sum += data[l] * m.getEntry(k, j); + } + out.data[o] += sum; + } + } + } + } + + return out; + + } + } + + /** + * Returns the result of postmultiplying this by m. + *

The Strassen matrix multiplication method is used here. This + * method computes C = A × B recursively by splitting all matrices + * in four quadrants and computing:

+ * 
+     * P1 = (A1,1 + A2,2) × (B1,1 + B2,2)
+     * P2 = (A2,1 + A2,2) × (B1,1)
+     * P3 = (A1,1) × (B1,2 - B2,2)
+     * P4 = (A2,2) × (B2,1 - B1,1)
+     * P5 = (A1,1 + A1,2) × B2,2
+     * P6 = (A2,1 - A1,1) × (B1,1 + B1,2)
+     * P7 = (A1,2 - A2,2) × (B2,1 + B2,2)
+     *
+     * C1,1 = P1 + P4 - P5 + P7
+     * C1,2 = P3 + P5
+     * C2,1 = P2 + P4
+     * C2,2 = P1 + P3 - P2 + P6
+     *
+ *

+ * This implementation is based on the 2002 paper: Recursive Array Layouts + * and Fast Matrix Multiplication by Siddhartha Chatterjee, Alvin R. Lebeck, + * Praveen K. Patnala and Mithuna Thottethodi. + *

+ * + * @param m matrix to postmultiply by + * @return this * m + * @throws IllegalArgumentException + * if columnDimension(this) != rowDimension(m) + */ + public RecursiveLayoutRealMatrix multiply(RecursiveLayoutRealMatrix m) + throws IllegalArgumentException { + + // safety check + checkMultiplicationCompatible(m); + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, m.columns); + if ((tileNumber != m.tileNumber) || (tileNumber != out.tileNumber)) { + // TODO get rid of this test + throw new RuntimeException("multiplication " + rows + "x" + columns + " * " + + m.rows + "x" + m.columns + " -> left matrix: " + tileNumber + + " tiles, right matrix: " + m.tileNumber + " tiles, result matrix " + + out.tileNumber + " tiles"); + } + strassenMultiply(data, 0, true, m.data, 0, true, out.data, 0, tileNumber, + tileSizeRows, m.tileSizeColumns, tileSizeColumns); + + return out; + + } + + /** + * Perform recursive multiplication using Strassen's algorithm. + * @param a left term of multiplication + * @param aStart start index in a + * @param aDirect direct/reversed orientation flag for a + * @param b right term of multiplication + * @param bStart start index in b + * @param bDirect direct/reversed orientation flag for b + * @param result result array (will have same orientation as b) + * @param resultStart start index in result + * @param nTiles number of elements to add + * @param bsRows number of rows in result tiles + * @param bsColumns number of columns in result tiles + * @param bsMultiplicands number of rows/columns in multiplicands + */ + private static void strassenMultiply(final double[] a, final int aStart, final boolean aDirect, + final double[] b, final int bStart, final boolean bDirect, + final double[] result, final int resultStart, final int nTiles, + final int bsRows, final int bsColumns, final int bsMultiplicands) { + if (nTiles == 1) { + // leaf recursion tile: perform traditional multiplication + final int bsColumns2 = 2 * bsColumns; + final int bsColumns3 = 3 * bsColumns; + final int bsColumns4 = 4 * bsColumns; + for (int i = 0; i < bsRows; ++i) { + for (int j = 0; j < bsColumns; ++j) { + double sum = 0; + int k = 0; + int aK = aStart + i * bsMultiplicands; + int bK = bStart + j; + while (k < bsMultiplicands - 3) { + sum += a[aK] * b[bK] + + a[aK + 1] * b[bK + bsColumns] + + a[aK + 2] * b[bK + bsColumns2] + + a[aK + 3] * b[bK + bsColumns3]; + k += 4; + aK += 4; + bK += bsColumns4; + } + while (k < bsMultiplicands) { + sum += a[aK] * b[bK]; + k += 1; + aK += 1; + bK += bsColumns; + } + result[resultStart + i * bsColumns + j] = sum; + } + } + } else { + // regular recursion node: use recursive Strassen implementation + final int n2 = nTiles / 2; + final int aQuadrantSize = bsRows * n2 * bsMultiplicands * n2; + final int bQuadrantSize = bsMultiplicands * n2 * bsColumns * n2; + final int cQuadrantSize = bsRows * n2 * bsColumns * n2; + final double[] sA = new double[aQuadrantSize]; + final double[] sB = new double[bQuadrantSize]; + final boolean nonLeafQuadrants = n2 > 1; + + // identify A quadrants start indices + final int a11Start, a12Start, a21Start, a22Start; + if (aDirect) { + a11Start = aStart; + a12Start = aStart + aQuadrantSize; + a21Start = aStart + 3 * aQuadrantSize; + a22Start = aStart + 2 * aQuadrantSize; + } else { + a11Start = aStart + 2 * aQuadrantSize; + a12Start = aStart + 3 * aQuadrantSize; + a21Start = aStart + aQuadrantSize; + a22Start = aStart; + } + + // identify B and C quadrants start indices + // (C is constructed with the same orientation as B) + final int b11Start, b12Start, b21Start, b22Start; + final int c11Start, c12Start, c21Start, c22Start; + if (bDirect) { + b11Start = bStart; + b12Start = bStart + bQuadrantSize; + b21Start = bStart + 3 * bQuadrantSize; + b22Start = bStart + 2 * bQuadrantSize; + c11Start = resultStart; + c12Start = resultStart + cQuadrantSize; + c21Start = resultStart + 3 * cQuadrantSize; + c22Start = resultStart + 2 * cQuadrantSize; + } else { + b11Start = bStart + 2 * bQuadrantSize; + b12Start = bStart + 3 * bQuadrantSize; + b21Start = bStart + bQuadrantSize; + b22Start = bStart; + c11Start = resultStart + 2 * cQuadrantSize; + c12Start = resultStart + 3 * cQuadrantSize; + c21Start = resultStart + cQuadrantSize; + c22Start = resultStart; + } + + // optimal order for cache efficiency: P3, P6, P2, P1, P5, P7, P4 + + // P3 = (A11)(B12 - B22) + // C12 = P3 + ... + tilesSubtract(b, b12Start, false, b, b22Start, false, sB, 0, + bQuadrantSize, nonLeafQuadrants); + strassenMultiply(a, a11Start, true, sB, 0, false, result, c12Start, + n2, bsRows, bsColumns, bsMultiplicands); + + // P6 = (A21 - A11)(B11 + B12) + // C22 = P3 + P6 + ... + final double[] p67 = new double[cQuadrantSize]; + tilesSubtract(a, a21Start, true, a, a11Start, true, sA, 0, + aQuadrantSize, nonLeafQuadrants); + tilesAdd(b, b11Start, true, b, b12Start, false, sB, 0, + bQuadrantSize, nonLeafQuadrants); + strassenMultiply(sA, 0, true, sB, 0, true, p67, 0, + n2, bsRows, bsColumns, bsMultiplicands); + tilesAdd(result, c12Start, false, p67, 0, true, result, c22Start, + cQuadrantSize, nonLeafQuadrants); + + // P2 = (A21 + A22)(B11) + // C21 = P2 + ... + // C22 = P3 + P6 - P2 + ... + tilesAdd(a, a21Start, true, a, a22Start, false, sA, 0, + aQuadrantSize, nonLeafQuadrants); + strassenMultiply(sA, 0, true, b, b11Start, true, result, c21Start, + n2, bsRows, bsColumns, bsMultiplicands); + tilesSelfSubtract(result, c22Start, false, result, c21Start, true, + cQuadrantSize, nonLeafQuadrants); + + // P1 = (A11 + A22)(B11 + B22) + // C11 = P1 + ... + // C22 = P3 + P6 - P2 + P1 + tilesAdd(a, a11Start, true, a, a22Start, false, sA, 0, + aQuadrantSize, nonLeafQuadrants); + tilesAdd(b, b11Start, true, b, b22Start, false, sB, 0, + bQuadrantSize, nonLeafQuadrants); + strassenMultiply(sA, 0, true, sB, 0, true, result, c11Start, + n2, bsRows, bsColumns, bsMultiplicands); + tilesSelfAdd(result, c22Start, false, result, c11Start, true, + cQuadrantSize, nonLeafQuadrants); + + // P5 = (A11 + A12)B22 + // beware: there is a sign error here in Chatterjee et al. paper + // in figure 1, table b they subtract A12 from A11 instead of adding it + // C12 = P3 + P5 + // C11 = P1 - P5 + ... + final double[] p45 = new double[cQuadrantSize]; + tilesAdd(a, a11Start, true, a, a12Start, false, sA, 0, + aQuadrantSize, nonLeafQuadrants); + strassenMultiply(sA, 0, true, b, b22Start, false, p45, 0, + n2, bsRows, bsColumns, bsMultiplicands); + tilesSelfAdd(result, c12Start, false, p45, 0, false, + cQuadrantSize, nonLeafQuadrants); + tilesSelfSubtract(result, c11Start, true, p45, 0, false, + cQuadrantSize, nonLeafQuadrants); + + // P7 = (A12 - A22)(B21 + B22) + // C11 = P1 - P5 + P7 + ... + tilesSubtract(a, a12Start, false, a, a22Start, false, sA, 0, + aQuadrantSize, nonLeafQuadrants); + tilesAdd(b, b21Start, true, b, b22Start, false, sB, 0, + bQuadrantSize, nonLeafQuadrants); + strassenMultiply(sA, 0, false, sB, 0, true, p67, 0, + n2, bsRows, bsColumns, bsMultiplicands); + tilesSelfAdd(result, c11Start, true, p67, 0, true, + cQuadrantSize, nonLeafQuadrants); + + // P4 = (A22)(B21 - B11) + // C11 = P1 - P5 + P7 + P4 + // C21 = P2 + P4 + tilesSubtract(b, b21Start, true, b, b11Start, true, sB, 0, + bQuadrantSize, nonLeafQuadrants); + strassenMultiply(a, a22Start, false, sB, 0, true, p45, 0, + n2, bsRows, bsColumns, bsMultiplicands); + tilesSelfAdd(result, c11Start, true, p45, 0, true, + cQuadrantSize, nonLeafQuadrants); + tilesSelfAdd(result, c21Start, true, p45, 0, true, + cQuadrantSize, nonLeafQuadrants); + + } + } + + /** + * Perform an addition on a few tiles in arrays. + * @param a left term of addition + * @param aStart start index in a + * @param aDirect direct/reversed orientation flag for a + * @param b right term of addition + * @param bStart start index in b + * @param bDirect direct/reversed orientation flag for b + * @param result result array (will have same orientation as a) + * @param resultStart start index in result + * @param n number of elements to add + * @param nonLeafQuadrants if true the quadrant can be further decomposed + */ + private static void tilesAdd(final double[] a, final int aStart, final boolean aDirect, + final double[] b, final int bStart, final boolean bDirect, + final double[] result, final int resultStart, + final int n, final boolean nonLeafQuadrants) { + if ((aDirect ^ bDirect) & nonLeafQuadrants) { + // a and b have different orientations + // perform addition in two half + final int n2 = n / 2; + addLoop(a, aStart, b, bStart + n2, result, resultStart, n2); + addLoop(a, aStart + n2, b, bStart, result, resultStart + n2, n2); + } else { + // a and b have same orientations + // perform addition in one loop + addLoop(a, aStart, b, bStart, result, resultStart, n); + } + } + + /** + * Perform an addition loop. + * @param a left term of addition + * @param aStart start index in a + * @param b right term of addition + * @param bStart start index in b + * @param result result array (will have same orientation as a) + * @param resultStart start index in result + * @param n number of elements to add + */ + private static void addLoop(final double[] a, final int aStart, + final double[] b, final int bStart, + final double[] result, final int resultStart, + final int n) { + int i = 0; + while (i < n - 3) { + final int r0 = resultStart + i; + final int a0 = aStart + i; + final int b0 = bStart + i; + result[r0] = a[a0] + b[b0]; + result[r0 + 1] = a[a0 + 1] + b[b0 + 1]; + result[r0 + 2] = a[a0 + 2] + b[b0 + 2]; + result[r0 + 3] = a[a0 + 3] + b[b0 + 3]; + i += 4; + } + while (i < n) { + result[resultStart + i] = a[aStart + i] + b[bStart + i]; + ++i; + } + } + + /** + * Perform a subtraction on a few tiles in arrays. + * @param a left term of subtraction + * @param aStart start index in a + * @param aDirect direct/reversed orientation flag for a + * @param b right term of subtraction + * @param bStart start index in b + * @param bDirect direct/reversed orientation flag for b + * @param result result array (will have same orientation as a) + * @param resultStart start index in result + * @param n number of elements to subtract + * @param nonLeafQuadrants if true the quadrant can be further decomposed + */ + private static void tilesSubtract(final double[] a, final int aStart, final boolean aDirect, + final double[] b, final int bStart, final boolean bDirect, + final double[] result, final int resultStart, + final int n, final boolean nonLeafQuadrants) { + if ((aDirect ^ bDirect) & nonLeafQuadrants) { + // a and b have different orientations + // perform subtraction in two half + final int n2 = n / 2; + subtractLoop(a, aStart, b, bStart + n2, result, resultStart, n2); + subtractLoop(a, aStart + n2, b, bStart, result, resultStart + n2, n2); + } else { + // a and b have same orientations + // perform subtraction in one loop + subtractLoop(a, aStart, b, bStart, result, resultStart, n); + } + } + + /** + * Perform a subtraction loop. + * @param a left term of subtraction + * @param aStart start index in a + * @param b right term of subtraction + * @param bStart start index in b + * @param result result array (will have same orientation as a) + * @param resultStart start index in result + * @param n number of elements to subtract + */ + private static void subtractLoop(final double[] a, final int aStart, + final double[] b, final int bStart, + final double[] result, final int resultStart, + final int n) { + int i = 0; + while (i < n - 3) { + final int r0 = resultStart + i; + final int a0 = aStart + i; + final int b0 = bStart + i; + result[r0] = a[a0] - b[b0]; + result[r0 + 1] = a[a0 + 1] - b[b0 + 1]; + result[r0 + 2] = a[a0 + 2] - b[b0 + 2]; + result[r0 + 3] = a[a0 + 3] - b[b0 + 3]; + i += 4; + } + while (i < n) { + result[resultStart + i] = a[aStart + i] - b[bStart + i]; + ++i; + } + } + + /** + * Perform a self-addition on a few tiles in arrays. + * @param a left term of addition (will be overwritten with result) + * @param aStart start index in a + * @param aDirect direct/reversed orientation flag for a + * @param b right term of addition + * @param bStart start index in b + * @param bDirect direct/reversed orientation flag for b + * @param n number of elements to add + * @param nonLeafQuadrants if true the quadrant can be further decomposed + */ + private static void tilesSelfAdd(final double[] a, final int aStart, final boolean aDirect, + final double[] b, final int bStart, final boolean bDirect, + final int n, final boolean nonLeafQuadrants) { + if ((aDirect ^ bDirect) & nonLeafQuadrants) { + // a and b have different orientations + // perform addition in two half + final int n2 = n / 2; + selfAddLoop(a, aStart, b, bStart + n2, n2); + selfAddLoop(a, aStart + n2, b, bStart, n2); + } else { + // a and b have same orientations + // perform addition in one loop + selfAddLoop(a, aStart, b, bStart, n); + } + } + + /** + * Perform a self-addition loop. + * @param a left term of addition (will be overwritten with result) + * @param aStart start index in a + * @param b right term of addition + * @param bStart start index in b + * @param n number of elements to add + */ + private static void selfAddLoop(final double[] a, final int aStart, + final double[] b, final int bStart, + final int n) { + int i = 0; + while (i < n - 3) { + final int a0 = aStart + i; + final int b0 = bStart + i; + a[a0] += b[b0]; + a[a0 + 1] += b[b0 + 1]; + a[a0 + 2] += b[b0 + 2]; + a[a0 + 3] += b[b0 + 3]; + i += 4; + } + while (i < n) { + a[aStart + i] += b[bStart + i]; + ++i; + } + } + + /** + * Perform a self-subtraction on a few tiles in arrays. + * @param a left term of subtraction (will be overwritten with result) + * @param aStart start index in a + * @param aDirect direct/reversed orientation flag for a + * @param b right term of subtraction + * @param bStart start index in b + * @param bDirect direct/reversed orientation flag for b + * @param n number of elements to subtract + * @param nonLeafQuadrants if true the quadrant can be further decomposed + */ + private static void tilesSelfSubtract(final double[] a, final int aStart, final boolean aDirect, + final double[] b, final int bStart, final boolean bDirect, + final int n, final boolean nonLeafQuadrants) { + if ((aDirect ^ bDirect) & nonLeafQuadrants) { + // a and b have different orientations + // perform subtraction in two half + final int n2 = n / 2; + selfSubtractLoop(a, aStart, b, bStart + n2, n2); + selfSubtractLoop(a, aStart + n2, b, bStart, n2); + } else { + // a and b have same orientations + // perform subtraction in one loop + selfSubtractLoop(a, aStart, b, bStart, n); + } + } + + /** + * Perform a self-subtraction loop. + * @param a left term of subtraction (will be overwritten with result) + * @param aStart start index in a + * @param b right term of subtraction + * @param bStart start index in b + * @param n number of elements to subtract + */ + private static void selfSubtractLoop(final double[] a, final int aStart, + final double[] b, final int bStart, + final int n) { + int i = 0; + while (i < n - 3) { + final int a0 = aStart + i; + final int b0 = bStart + i; + a[a0] -= b[b0]; + a[a0 + 1] -= b[b0 + 1]; + a[a0 + 2] -= b[b0 + 2]; + a[a0 + 3] -= b[b0 + 3]; + i += 4; + } + while (i < n) { + a[aStart + i] -= b[bStart + i]; + ++i; + } + } + + /** {@inheritDoc} */ + public double[][] getData() { + + final double[][] out = new double[rows][columns]; + + // perform extraction tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + + // perform extraction on the current tile + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int qStart = jTile * tileSizeColumns; + if (pStart < rows && qStart < columns) { + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + int tileRowStart = tileStart; + for (int p = pStart; p < pEnd; ++p) { + System.arraycopy(data, tileRowStart, out[p], qStart, qEnd - qStart); + tileRowStart += tileSizeColumns; + } + } + + } + + return out; + + } + + /** {@inheritDoc} */ + public double getFrobeniusNorm() { + double sum2 = 0; + for (final double entry : data) { + sum2 += entry * entry; + } + return Math.sqrt(sum2); + } + + /** {@inheritDoc} */ + public RealMatrix getSubMatrix(final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException { + + // safety checks + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + + // create the output matrix + final RecursiveLayoutRealMatrix out = + new RecursiveLayoutRealMatrix(endRow - startRow + 1, endColumn - startColumn + 1); + + // perform extraction tile-wise, to ensure good cache behavior + for (int iTile = 0; iTile < out.tileNumber; ++iTile) { + final int iStart = startRow + iTile * out.tileSizeRows; + final int iEnd = Math.min(startRow + Math.min((iTile + 1) * out.tileSizeRows, out.rows), + endRow + 1); + for (int jTile = 0; jTile < out.tileNumber; ++jTile) { + final int jStart = startColumn + jTile * out.tileSizeColumns; + final int jEnd = Math.min(startColumn + Math.min((jTile + 1) * out.tileSizeColumns, out.columns), + endColumn + 1); + + // the current output tile may expand on more than one instance tile + for (int pTile = iStart / tileSizeRows; pTile * tileSizeRows < iEnd; ++pTile) { + final int p0 = pTile * tileSizeRows; + final int pStart = Math.max(p0, iStart); + final int pEnd = Math.min(Math.min(p0 + tileSizeRows, endRow + 1), iEnd); + for (int qTile = jStart / tileSizeColumns; qTile * tileSizeColumns < jEnd; ++qTile) { + final int q0 = qTile * tileSizeColumns; + final int qStart = Math.max(q0, jStart); + final int qEnd = Math.min(Math.min(q0 + tileSizeColumns, endColumn + 1), jEnd); + + // copy the overlapping part of instance and output tiles + int outIndex = tileIndex(iTile, jTile) * out.tileSizeRows * out.tileSizeColumns + + (pStart - iStart) * out.tileSizeColumns + (qStart - jStart); + int index = tileIndex(pTile, qTile) * tileSizeRows * tileSizeColumns + + (pStart - p0) * tileSizeColumns + (qStart - q0); + for (int p = pStart; p < pEnd; ++p) { + System.arraycopy(data, index, out.data, outIndex, qEnd - qStart); + outIndex += out.tileSizeColumns; + index += tileSizeColumns; + } + + + } + } + + } + } + + return out; + + } + + /** {@inheritDoc} */ + public void setSubMatrix(final double[][] subMatrix, final int row, final int column) + throws MatrixIndexException { + + // safety checks + final int refLength = subMatrix[0].length; + if (refLength < 1) { + throw MathRuntimeException.createIllegalArgumentException("matrix must have at least one column", + null); + } + final int endRow = row + subMatrix.length - 1; + final int endColumn = column + refLength - 1; + checkSubMatrixIndex(row, endRow, column, endColumn); + for (final double[] subRow : subMatrix) { + if (subRow.length != refLength) { + throw MathRuntimeException.createIllegalArgumentException("some rows have length {0} while others have length {1}", + new Object[] { + refLength, subRow.length + }); + } + } + + // compute tiles bounds + final int tileStartRow = row / tileSizeRows; + final int tileEndRow = (endRow + tileSizeRows) / tileSizeRows; + final int tileStartColumn = column / tileSizeColumns; + final int tileEndColumn = (endColumn + tileSizeColumns) / tileSizeColumns; + + // perform copy tile-wise, to ensure good cache behavior + for (int iTile = tileStartRow; iTile < tileEndRow; ++iTile) { + final int firstRow = iTile * tileSizeRows; + final int iStart = Math.max(row, firstRow); + final int iEnd = Math.min(endRow + 1, firstRow + tileSizeRows); + + for (int jTile = tileStartColumn; jTile < tileEndColumn; ++jTile) { + final int firstColumn = jTile * tileSizeColumns; + final int jStart = Math.max(column, firstColumn); + final int jEnd = Math.min(endColumn + 1, firstColumn + tileSizeColumns); + final int jLength = jEnd - jStart; + final int tileStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns; + + // handle one tile, row by row + for (int i = iStart; i < iEnd; ++i) { + System.arraycopy(subMatrix[i - row], jStart - column, + data, tileStart + (i - firstRow) * tileSizeColumns + (jStart - firstColumn), + jLength); + } + + } + } + } + + /** {@inheritDoc} */ + public RealMatrix getRowMatrix(final int row) + throws MatrixIndexException { + + checkRowIndex(row); + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(1, columns); + + // a row matrix has always only one large tile, + // because a single row cannot be split into 2^k tiles + // perform copy tile-wise, to ensure good cache behavior + final int iTile = row / tileSizeRows; + final int rowOffset = row - iTile * tileSizeRows; + int outIndex = 0; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + rowOffset * tileSizeColumns; + final int length = Math.min(outIndex + tileSizeColumns, columns) - outIndex; + System.arraycopy(data, kStart, out.data, outIndex, length); + outIndex += length; + } + + return out; + + } + + /** {@inheritDoc} */ + public void setRowMatrix(final int row, final RealMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + try { + setRowMatrix(row, (RecursiveLayoutRealMatrix) matrix); + } catch (ClassCastException cce) { + super.setRowMatrix(row, matrix); + } + } + + /** + * Sets the entries in row number row + * as a row matrix. Row indices start at 0. + * + * @param row the row to be set + * @param matrix row matrix (must have one row and the same number of columns + * as the instance) + * @throws MatrixIndexException if the specified row index is invalid + * @throws InvalidMatrixException if the matrix dimensions do not match one + * instance row + */ + public void setRowMatrix(final int row, final RecursiveLayoutRealMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + + checkRowIndex(row); + final int nCols = getColumnDimension(); + if ((matrix.getRowDimension() != 1) || + (matrix.getColumnDimension() != nCols)) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + matrix.getRowDimension(), + matrix.getColumnDimension(), + 1, nCols + }); + } + + // a row matrix has always only one large tile, + // because a single row cannot be split into 2^k tiles + // perform copy tile-wise, to ensure good cache behavior + final int iTile = row / tileSizeRows; + final int rowOffset = row - iTile * tileSizeRows; + int outIndex = 0; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + rowOffset * tileSizeColumns; + final int length = Math.min(outIndex + tileSizeColumns, columns) - outIndex; + System.arraycopy(matrix.data, outIndex, data, kStart, length); + outIndex += length; + } + + } + + /** {@inheritDoc} */ + public RealMatrix getColumnMatrix(final int column) + throws MatrixIndexException { + + checkColumnIndex(column); + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(rows, 1); + + // a column matrix has always only one large tile, + // because a single column cannot be split into 2^k tiles + // perform copy tile-wise, to ensure good cache behavior + final int jTile = column / tileSizeColumns; + final int columnOffset = column - jTile * tileSizeColumns; + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + columnOffset; + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + out.data[p] = data[k]; + } + } + + return out; + + } + + /** {@inheritDoc} */ + public void setColumnMatrix(final int column, final RealMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + try { + setColumnMatrix(column, (RecursiveLayoutRealMatrix) matrix); + } catch (ClassCastException cce) { + super.setColumnMatrix(column, matrix); + } + } + + /** + * Sets the entries in column number column + * as a column matrix. Column indices start at 0. + * + * @param column the column to be set + * @param matrix column matrix (must have one column and the same number of rows + * as the instance) + * @throws MatrixIndexException if the specified column index is invalid + * @throws InvalidMatrixException if the matrix dimensions do not match one + * instance column + */ + void setColumnMatrix(final int column, final RecursiveLayoutRealMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + + checkColumnIndex(column); + final int nRows = getRowDimension(); + if ((matrix.getRowDimension() != nRows) || + (matrix.getColumnDimension() != 1)) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + matrix.getRowDimension(), + matrix.getColumnDimension(), + nRows, 1 + }); + } + + // a column matrix has always only one large tile, + // because a single column cannot be split into 2^k tiles + // perform copy tile-wise, to ensure good cache behavior + final int jTile = column / tileSizeColumns; + final int columnOffset = column - jTile * tileSizeColumns; + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + columnOffset; + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + data[k] = matrix.data[p]; + } + } + + } + + /** {@inheritDoc} */ + public void setRowVector(final int row, final RealVector vector) + throws MatrixIndexException, InvalidMatrixException { + try { + setRow(row, ((RealVectorImpl) vector).getDataRef()); + } catch (ClassCastException cce) { + checkRowIndex(row); + if (vector.getDimension() != columns) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + 1, vector.getDimension(), + 1, columns + }); + } + + // perform copy tile-wise, to ensure good cache behavior + final int iTile = row / tileSizeRows; + final int rowOffset = row - iTile * tileSizeRows; + int outIndex = 0; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + rowOffset * tileSizeColumns; + final int length = Math.min(outIndex + tileSizeColumns, columns) - outIndex; + for (int l = 0; l < length; ++l) { + data[kStart + l] = vector.getEntry(outIndex + l); + } + outIndex += length; + } + } + } + + /** {@inheritDoc} */ + public void setColumnVector(final int column, final RealVector vector) + throws MatrixIndexException, InvalidMatrixException { + try { + setColumn(column, ((RealVectorImpl) vector).getDataRef()); + } catch (ClassCastException cce) { + checkColumnIndex(column); + if (vector.getDimension() != rows) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + vector.getDimension(), 1, + rows, 1 + }); + } + + // perform copy tile-wise, to ensure good cache behavior + final int jTile = column / tileSizeColumns; + final int columnOffset = column - jTile * tileSizeColumns; + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + columnOffset; + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + data[k] = vector.getEntry(p); + } + } + } + } + + /** {@inheritDoc} */ + public double[] getRow(final int row) + throws MatrixIndexException { + + checkRowIndex(row); + final double[] out = new double[columns]; + + // perform copy tile-wise, to ensure good cache behavior + final int iTile = row / tileSizeRows; + final int rowOffset = row - iTile * tileSizeRows; + int outIndex = 0; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + rowOffset * tileSizeColumns; + final int length = Math.min(outIndex + tileSizeColumns, columns) - outIndex; + System.arraycopy(data, kStart, out, outIndex, length); + outIndex += length; + } + + return out; + + } + + /** {@inheritDoc} */ + public void setRow(final int row, final double[] array) + throws MatrixIndexException, InvalidMatrixException { + + checkRowIndex(row); + if (array.length != columns) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + 1, array.length, + 1, columns + }); + } + + // perform copy tile-wise, to ensure good cache behavior + final int iTile = row / tileSizeRows; + final int rowOffset = row - iTile * tileSizeRows; + int outIndex = 0; + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + rowOffset * tileSizeColumns; + final int length = Math.min(outIndex + tileSizeColumns, columns) - outIndex; + System.arraycopy(array, outIndex, data, kStart, length); + outIndex += length; + } + + } + + /** {@inheritDoc} */ + public double[] getColumn(final int column) + throws MatrixIndexException { + + checkColumnIndex(column); + final double[] out = new double[rows]; + + // perform copy tile-wise, to ensure good cache behavior + final int jTile = column / tileSizeColumns; + final int columnOffset = column - jTile * tileSizeColumns; + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + columnOffset; + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + out[p] = data[k]; + } + } + + return out; + + } + + /** {@inheritDoc} */ + public void setColumn(final int column, final double[] array) + throws MatrixIndexException, InvalidMatrixException { + + checkColumnIndex(column); + if (array.length != rows) { + throw new InvalidMatrixException("dimensions mismatch: got {0}x{1} but expected {2}x{3}", + new Object[] { + array.length, 1, + rows, 1 + }); + } + + // perform copy tile-wise, to ensure good cache behavior + final int jTile = column / tileSizeColumns; + final int columnOffset = column - jTile * tileSizeColumns; + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int kStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns + + columnOffset; + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + data[k] = array[p]; + } + } + + } + + /** {@inheritDoc} */ + public double getEntry(final int row, final int column) + throws MatrixIndexException { + if ((row < 0) || (row >= rows) || (column < 0) || (column >= columns)) { + throw new MatrixIndexException("no entry at indices ({0}, {1}) in a {2}x{3} matrix", + new Object[] { + row, column, + getRowDimension(), getColumnDimension() + }); + } + return data[index(row, column)]; + } + + /** {@inheritDoc} */ + public void setEntry(final int row, final int column, final double value) + throws MatrixIndexException { + if ((row < 0) || (row >= rows) || (column < 0) || (column >= columns)) { + throw new MatrixIndexException("no entry at indices ({0}, {1}) in a {2}x{3} matrix", + new Object[] { + row, column, + getRowDimension(), getColumnDimension() + }); + } + data[index(row, column)] = value; + } + + /** {@inheritDoc} */ + public void addToEntry(final int row, final int column, final double increment) + throws MatrixIndexException { + if ((row < 0) || (row >= rows) || (column < 0) || (column >= columns)) { + throw new MatrixIndexException("no entry at indices ({0}, {1}) in a {2}x{3} matrix", + new Object[] { + row, column, + getRowDimension(), getColumnDimension() + }); + } + data[index(row, column)] += increment; + } + + /** {@inheritDoc} */ + public void multiplyEntry(final int row, final int column, final double factor) + throws MatrixIndexException { + if ((row < 0) || (row >= rows) || (column < 0) || (column >= columns)) { + throw new MatrixIndexException("no entry at indices ({0}, {1}) in a {2}x{3} matrix", + new Object[] { + row, column, + getRowDimension(), getColumnDimension() + }); + } + data[index(row, column)] *= factor; + } + + /** {@inheritDoc} */ + public RealMatrix transpose() { + + final RecursiveLayoutRealMatrix out = new RecursiveLayoutRealMatrix(columns, rows); + + // perform transpose tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int outJTile = (int) (indices >> 32); // iTile in the instance + final int outITile = (int) (indices & 0xffffffff); // jTile in the instance + final int outIndex = tileIndex(outITile, outJTile); + final int outTileStart = outIndex * tileSizeRows * tileSizeColumns; + + // transpose current tile + final int outPStart = outITile * tileSizeColumns; + final int outPEnd = Math.min(outPStart + tileSizeColumns, columns); + final int outQStart = outJTile * tileSizeRows; + final int outQEnd = Math.min(outQStart + tileSizeRows, rows); + for (int outP = outPStart; outP < outPEnd; ++outP) { + final int dP = outP - outPStart; + int k = outTileStart + dP * tileSizeRows; + int l = tileStart + dP; + for (int outQ = outQStart; outQ < outQEnd; ++outQ) { + out.data[k++] = data[l]; + l+= tileSizeColumns; + } + } + + } + + return out; + + } + + /** {@inheritDoc} */ + public int getRowDimension() { + return rows; + } + + /** {@inheritDoc} */ + public int getColumnDimension() { + return columns; + } + + /** {@inheritDoc} */ + public double[] operate(final double[] v) + throws IllegalArgumentException { + + if (v.length != columns) { + throw MathRuntimeException.createIllegalArgumentException("vector length mismatch:" + + " got {0} but expected {1}", + new Object[] { + v.length, columns + }); + } + final double[] out = new double[rows]; + + // perform multiplication tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int p = pStart, k = tileStart; p < pEnd; ++p) { + double sum = 0; + int q = qStart; + while (q < qEnd - 3) { + sum += data[k] * v[q] + + data[k + 1] * v[q + 1] + + data[k + 2] * v[q + 2] + + data[k + 3] * v[q + 3]; + k += 4; + q += 4; + } + while (q < qEnd) { + sum += data[k++] * v[q++]; + } + out[p] += sum; + } + } + + return out; + + } + + /** {@inheritDoc} */ + public double[] preMultiply(final double[] v) + throws IllegalArgumentException { + + if (v.length != rows) { + throw MathRuntimeException.createIllegalArgumentException("vector length mismatch:" + + " got {0} but expected {1}", + new Object[] { + v.length, rows + }); + } + final double[] out = new double[columns]; + + final int offset1 = tileSizeColumns; + final int offset2 = offset1 + offset1; + final int offset3 = offset2 + offset1; + final int offset4 = offset3 + offset1; + + // perform multiplication tile-wise, to ensure good cache behavior + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int q = qStart; q < qEnd; ++q) { + int k = tileStart + q - qStart; + double sum = 0; + int p = pStart; + while (p < pEnd - 3) { + sum += data[k] * v[p] + + data[k + offset1] * v[p + 1] + + data[k + offset2] * v[p + 2] + + data[k + offset3] * v[p + 3]; + k += offset4; + p += 4; + } + while (p < pEnd) { + sum += data[k] * v[p++]; + k += offset1; + } + out[q] += sum; + } + } + + return out; + + } + + /** {@inheritDoc} */ + public double walkInRowOrder(final RealMatrixChangingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + for (int p = pStart; p < pEnd; ++p) { + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInRowOrder(final RealMatrixPreservingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + for (int p = pStart; p < pEnd; ++p) { + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInRowOrder(final RealMatrixChangingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(rows, columns, startRow, endRow, startColumn, endColumn); + for (int iTile = startRow / tileSizeRows; iTile < 1 + endRow / tileSizeRows; ++iTile) { + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + for (int p = pStart; p < pEnd; ++p) { + for (int jTile = startColumn / tileSizeColumns; jTile < 1 + endColumn / tileSizeColumns; ++jTile) { + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (p - p0) * tileSizeColumns + (qStart - q0); + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInRowOrder(final RealMatrixPreservingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(rows, columns, startRow, endRow, startColumn, endColumn); + for (int iTile = startRow / tileSizeRows; iTile < 1 + endRow / tileSizeRows; ++iTile) { + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + for (int p = pStart; p < pEnd; ++p) { + for (int jTile = startColumn / tileSizeColumns; jTile < 1 + endColumn / tileSizeColumns; ++jTile) { + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (p - p0) * tileSizeColumns + (qStart - q0); + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInColumnOrder(final RealMatrixChangingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int q = qStart; q < qEnd; ++q) { + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (q - qStart); + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInColumnOrder(final RealMatrixPreservingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int q = qStart; q < qEnd; ++q) { + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (q - qStart); + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInColumnOrder(final RealMatrixChangingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(getRowDimension(), getColumnDimension(), + startRow, endRow, startColumn, endColumn); + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + for (int q = qStart; q < qEnd; ++q) { + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (pStart - p0) * tileSizeColumns + (q - q0); + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInColumnOrder(final RealMatrixPreservingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(getRowDimension(), getColumnDimension(), + startRow, endRow, startColumn, endColumn); + for (int jTile = 0; jTile < tileNumber; ++jTile) { + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + for (int q = qStart; q < qEnd; ++q) { + for (int iTile = 0; iTile < tileNumber; ++iTile) { + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + final int tileStart = tileIndex(iTile, jTile) * + tileSizeRows * tileSizeColumns; + final int kStart = tileStart + (pStart - p0) * tileSizeColumns + (q - q0); + for (int p = pStart, k = kStart; p < pEnd; ++p, k += tileSizeColumns) { + visitor.visit(p, q, data[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int pStart = iTile * tileSizeRows; + final int pEnd = Math.min(pStart + tileSizeRows, rows); + final int qStart = jTile * tileSizeColumns; + final int qEnd = Math.min(qStart + tileSizeColumns, columns); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - pStart) * tileSizeColumns; + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + visitor.visit(p, q, data[k]); + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(rows, columns, startRow, endRow, startColumn, endColumn); + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - p0) * tileSizeColumns + (qStart - q0); + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + data[k] = visitor.visit(p, q, data[k]); + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor, + final int startRow, final int endRow, + final int startColumn, final int endColumn) + throws MatrixIndexException, MatrixVisitorException { + checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); + visitor.start(rows, columns, startRow, endRow, startColumn, endColumn); + for (int index = 0; index < tileNumber * tileNumber; ++index) { + final int tileStart = index * tileSizeRows * tileSizeColumns; + final long indices = tilesIndices(index); + final int iTile = (int) (indices >> 32); + final int jTile = (int) (indices & 0xffffffff); + final int p0 = iTile * tileSizeRows; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iTile + 1) * tileSizeRows, 1 + endRow); + final int q0 = jTile * tileSizeColumns; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jTile + 1) * tileSizeColumns, 1 + endColumn); + for (int p = pStart; p < pEnd; ++p) { + final int kStart = tileStart + (p - p0) * tileSizeColumns + (qStart - q0); + for (int q = qStart, k = kStart; q < qEnd; ++q, ++k) { + visitor.visit(p, q, data[k]); + } + } + } + return visitor.end(); + } + + /** + * Get the index of an element. + * @param row row index of the element + * @param column column index of the element + * @return index of the element + */ + private int index(final int row, final int columns) { + final int iTile = row / tileSizeRows; + final int jTile = columns / tileSizeColumns; + final int tileStart = tileIndex(iTile, jTile) * tileSizeRows * tileSizeColumns; + final int indexInTile = (row % tileSizeRows) * tileSizeColumns + + (columns % tileSizeColumns); + return tileStart + indexInTile; + } + + /** + * Get the index of a tile. + * @param iTile row index of the tile + * @param jTile column index of the tile + * @return index of the tile + */ + private static int tileIndex(int iTile, int jTile) { + + // compute n = 2^k such that a nxn square contains the indices + int n = Integer.highestOneBit(Math.max(iTile, jTile)) << 1; + + // start recursion by noting the index is somewhere in the nxn + // square whose lowest index is 0 and which has direct orientation + int lowIndex = 0; + boolean direct = true; + + // the tail-recursion on the square size is replaced by an iteration here + while (n > 1) { + + // reduce square to 4 quadrants + n >>= 1; + final int n2 = n * n; + + // check in which quadrant the element is, + // updating the lowest index of the quadrant and its orientation + if (iTile < n) { + if (jTile < n) { + // the element is in the top-left quadrant + if (!direct) { + lowIndex += 2 * n2; + direct = true; + } + } else { + // the element is in the top-right quadrant + jTile -= n; + if (direct) { + lowIndex += n2; + direct = false; + } else { + lowIndex += 3 * n2; + } + } + } else { + iTile -= n; + if (jTile < n) { + // the element is in the bottom-left quadrant + if (direct) { + lowIndex += 3 * n2; + } else { + lowIndex += n2; + direct = true; + } + } else { + // the element is in the bottom-right quadrant + jTile -= n; + if (direct) { + lowIndex += 2 * n2; + direct = false; + } + } + } + } + + // the lowest index of the remaining 1x1 quadrant is the requested index + return lowIndex; + + } + + /** + * Get the row and column tile indices of a tile. + * @param index index of the tile in the layout + * @return row and column indices packed in one long (row tile index + * in 32 high order bits, column tile index in low order bits) + */ + private static long tilesIndices(int index) { + + // compute n = 2^k such that a nxn square contains the index + int n = Integer.highestOneBit((int) Math.sqrt(index)) << 1; + + // start recursion by noting the index is somewhere in the nxn + // square whose lowest index is 0 and which has direct orientation + int iLow = 0; + int jLow = 0; + boolean direct = true; + + // the tail-recursion on the square size is replaced by an iteration here + while (n > 1) { + + // reduce square to 4 quadrants + n >>= 1; + final int n2 = n * n; + + // check in which quadrant the element is, + // updating the low indices of the quadrant and its orientation + switch (index / n2) { + case 0 : + if (!direct) { + iLow += n; + jLow += n; + } + break; + case 1 : + if (direct) { + jLow += n; + } else { + iLow += n; + } + index -= n2; + direct = !direct; + break; + case 2 : + if (direct) { + iLow += n; + jLow += n; + } + index -= 2 * n2; + direct = !direct; + break; + default : + if (direct) { + iLow += n; + } else { + jLow += n; + } + index -= 3 * n2; + } + + } + + // the lowest indices of the remaining 1x1 quadrant are the requested indices + return (((long) iLow) << 32) | (long) jLow; + + } + + /** + * Compute the power of two number of tiles for a matrix. + * @param rows number of rows + * @param columns number of columns + * @return power of two number of tiles + */ + private static int tilesNumber(final int rows, final int columns) { + + // find the minimal number of tiles, given that one double variable is 8 bytes + final int nbElements = rows * columns; + final int maxElementsPerTile = MAX_TILE_SIZE_BYTES / 8; + final int minTiles = nbElements / maxElementsPerTile; + + // the number of tiles must be a 2^k x 2^k square + int twoK = 1; + for (int nTiles = minTiles; nTiles != 0; nTiles >>= 2) { + twoK <<= 1; + } + + // make sure the tiles have at least one row and one column each + // (this may lead to tile sizes greater than MAX_BLOCK_SIZE_BYTES, + // in degenerate cases like a 3000x1 matrix) + while (twoK > Math.min(rows, columns)) { + twoK >>= 1; + } + + return twoK; + + } + + /** + * Compute optimal tile size for a row or column count. + * @param count row or column count + * @param twoK optimal tile number (must be a power of 2) + * @return optimal tile size + */ + private static int tileSize(final int count, final int twoK) { + return (count + twoK - 1) / twoK; + } + +} diff --git a/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrixTest.java b/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrixTest.java new file mode 100644 index 000000000..2e8719097 --- /dev/null +++ b/src/experimental/org/apache/commons/math/linear/RecursiveLayoutRealMatrixTest.java @@ -0,0 +1,1242 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.linear; + +import java.util.Arrays; +import java.util.Random; + +import junit.framework.Test; +import junit.framework.TestCase; +import junit.framework.TestSuite; + +/** + * Test cases for the {@link RecursiveLayoutRealMatrix} class. + * + * @version $Revision$ $Date$ + */ + +public final class RecursiveLayoutRealMatrixTest extends TestCase { + + // 3 x 3 identity matrix + protected double[][] id = { {1d,0d,0d}, {0d,1d,0d}, {0d,0d,1d} }; + + // Test data for group operations + protected double[][] testData = { {1d,2d,3d}, {2d,5d,3d}, {1d,0d,8d} }; + protected double[][] testDataLU = {{2d, 5d, 3d}, {.5d, -2.5d, 6.5d}, {0.5d, 0.2d, .2d}}; + protected double[][] testDataPlus2 = { {3d,4d,5d}, {4d,7d,5d}, {3d,2d,10d} }; + protected double[][] testDataMinus = { {-1d,-2d,-3d}, {-2d,-5d,-3d}, + {-1d,0d,-8d} }; + protected double[] testDataRow1 = {1d,2d,3d}; + protected double[] testDataCol3 = {3d,3d,8d}; + protected double[][] testDataInv = + { {-40d,16d,9d}, {13d,-5d,-3d}, {5d,-2d,-1d} }; + protected double[] preMultTest = {8,12,33}; + protected double[][] testData2 ={ {1d,2d,3d}, {2d,5d,3d}}; + protected double[][] testData2T = { {1d,2d}, {2d,5d}, {3d,3d}}; + protected double[][] testDataPlusInv = + { {-39d,18d,12d}, {15d,0d,0d}, {6d,-2d,7d} }; + + // lu decomposition tests + protected double[][] luData = { {2d,3d,3d}, {0d,5d,7d}, {6d,9d,8d} }; + protected double[][] luDataLUDecomposition = { {6d,9d,8d}, {0d,5d,7d}, + {0.33333333333333,0d,0.33333333333333} }; + + // singular matrices + protected double[][] singular = { {2d,3d}, {2d,3d} }; + protected double[][] bigSingular = {{1d,2d,3d,4d}, {2d,5d,3d,4d}, + {7d,3d,256d,1930d}, {3d,7d,6d,8d}}; // 4th row = 1st + 2nd + protected double[][] detData = { {1d,2d,3d}, {4d,5d,6d}, {7d,8d,10d} }; + protected double[][] detData2 = { {1d, 3d}, {2d, 4d}}; + + // vectors + protected double[] testVector = {1,2,3}; + protected double[] testVector2 = {1,2,3,4}; + + // submatrix accessor tests + protected double[][] subTestData = {{1, 2, 3, 4}, {1.5, 2.5, 3.5, 4.5}, + {2, 4, 6, 8}, {4, 5, 6, 7}}; + // array selections + protected double[][] subRows02Cols13 = { {2, 4}, {4, 8}}; + protected double[][] subRows03Cols12 = { {2, 3}, {5, 6}}; + protected double[][] subRows03Cols123 = { {2, 3, 4} , {5, 6, 7}}; + // effective permutations + protected double[][] subRows20Cols123 = { {4, 6, 8} , {2, 3, 4}}; + protected double[][] subRows31Cols31 = {{7, 5}, {4.5, 2.5}}; + // contiguous ranges + protected double[][] subRows01Cols23 = {{3,4} , {3.5, 4.5}}; + protected double[][] subRows23Cols00 = {{2} , {4}}; + protected double[][] subRows00Cols33 = {{4}}; + // row matrices + protected double[][] subRow0 = {{1,2,3,4}}; + protected double[][] subRow3 = {{4,5,6,7}}; + // column matrices + protected double[][] subColumn1 = {{2}, {2.5}, {4}, {5}}; + protected double[][] subColumn3 = {{4}, {4.5}, {8}, {7}}; + + // tolerances + protected double entryTolerance = 10E-16; + protected double normTolerance = 10E-14; + + public RecursiveLayoutRealMatrixTest(String name) { + super(name); + } + + public void setUp() { + + } + + public static Test suite() { + TestSuite suite = new TestSuite(RecursiveLayoutRealMatrixTest.class); + suite.setName("RecursiveLayoutRealMatrix Tests"); + return suite; + } + + /** test dimensions */ + public void testDimensions() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testData2); + assertEquals("testData row dimension",3,m.getRowDimension()); + assertEquals("testData column dimension",3,m.getColumnDimension()); + assertTrue("testData is square",m.isSquare()); + assertEquals("testData2 row dimension",m2.getRowDimension(),2); + assertEquals("testData2 column dimension",m2.getColumnDimension(),3); + assertTrue("testData2 is not square",!m2.isSquare()); + } + + /** test copy functions */ + public void testCopyFunctions() { + Random r = new Random(66636328996002l); + RecursiveLayoutRealMatrix m1 = createRandomMatrix(r, 47, 83); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(m1.getData()); + assertEquals(m1, m2); + RecursiveLayoutRealMatrix m3 = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m4 = new RecursiveLayoutRealMatrix(m3.getData()); + assertEquals(m3, m4); + } + + /** test add */ + public void testAdd() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix mInv = new RecursiveLayoutRealMatrix(testDataInv); + RealMatrix mPlusMInv = m.add(mInv); + double[][] sumEntries = mPlusMInv.getData(); + for (int row = 0; row < m.getRowDimension(); row++) { + for (int col = 0; col < m.getColumnDimension(); col++) { + assertEquals("sum entry entry", + testDataPlusInv[row][col],sumEntries[row][col], + entryTolerance); + } + } + } + + /** test add failure */ + public void testAddFail() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testData2); + try { + m.add(m2); + fail("IllegalArgumentException expected"); + } catch (IllegalArgumentException ex) { + ; + } + } + + /** test norm */ + public void testNorm() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testData2); + assertEquals("testData norm",14d,m.getNorm(),entryTolerance); + assertEquals("testData2 norm",7d,m2.getNorm(),entryTolerance); + } + + /** test Frobenius norm */ + public void testFrobeniusNorm() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testData2); + assertEquals("testData Frobenius norm", Math.sqrt(117.0), m.getFrobeniusNorm(), entryTolerance); + assertEquals("testData2 Frobenius norm", Math.sqrt(52.0), m2.getFrobeniusNorm(), entryTolerance); + } + + /** test m-n = m + -n */ + public void testPlusMinus() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testDataInv); + assertClose(m.subtract(m2), m2.scalarMultiply(-1d).add(m), entryTolerance); + try { + m.subtract(new RecursiveLayoutRealMatrix(testData2)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + ; + } + } + + /** test multiply */ + public void testMultiply() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix mInv = new RecursiveLayoutRealMatrix(testDataInv); + RecursiveLayoutRealMatrix identity = new RecursiveLayoutRealMatrix(id); + RecursiveLayoutRealMatrix m2 = new RecursiveLayoutRealMatrix(testData2); + assertClose(m.multiply(mInv), identity, entryTolerance); + assertClose(mInv.multiply(m), identity, entryTolerance); + assertClose(m.multiply(identity), m, entryTolerance); + assertClose(identity.multiply(mInv), mInv, entryTolerance); + assertClose(m2.multiply(identity), m2, entryTolerance); + try { + m.multiply(new RecursiveLayoutRealMatrix(bigSingular)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + // expected + } + } + + public void testSeveralBlocks() { + + RealMatrix m = new RecursiveLayoutRealMatrix(35, 71); + for (int i = 0; i < m.getRowDimension(); ++i) { + for (int j = 0; j < m.getColumnDimension(); ++j) { + m.setEntry(i, j, i + j / 1024.0); + } + } + + RealMatrix mT = m.transpose(); + assertEquals(m.getRowDimension(), mT.getColumnDimension()); + assertEquals(m.getColumnDimension(), mT.getRowDimension()); + for (int i = 0; i < mT.getRowDimension(); ++i) { + for (int j = 0; j < mT.getColumnDimension(); ++j) { + assertEquals(m.getEntry(j, i), mT.getEntry(i, j), 0); + } + } + + RealMatrix mPm = m.add(m); + for (int i = 0; i < mPm.getRowDimension(); ++i) { + for (int j = 0; j < mPm.getColumnDimension(); ++j) { + assertEquals(2 * m.getEntry(i, j), mPm.getEntry(i, j), 0); + } + } + + RealMatrix mPmMm = mPm.subtract(m); + for (int i = 0; i < mPmMm.getRowDimension(); ++i) { + for (int j = 0; j < mPmMm.getColumnDimension(); ++j) { + assertEquals(m.getEntry(i, j), mPmMm.getEntry(i, j), 0); + } + } + + RealMatrix mTm = mT.multiply(m); + for (int i = 0; i < mTm.getRowDimension(); ++i) { + for (int j = 0; j < mTm.getColumnDimension(); ++j) { + double sum = 0; + for (int k = 0; k < mT.getColumnDimension(); ++k) { + sum += (k + i / 1024.0) * (k + j / 1024.0); + } + assertEquals(sum, mTm.getEntry(i, j), 0); + } + } + + RealMatrix mmT = m.multiply(mT); + for (int i = 0; i < mmT.getRowDimension(); ++i) { + for (int j = 0; j < mmT.getColumnDimension(); ++j) { + double sum = 0; + for (int k = 0; k < m.getColumnDimension(); ++k) { + sum += (i + k / 1024.0) * (j + k / 1024.0); + } + assertEquals(sum, mmT.getEntry(i, j), 0); + } + } + + RealMatrix sub1 = m.getSubMatrix(2, 9, 5, 20); + for (int i = 0; i < sub1.getRowDimension(); ++i) { + for (int j = 0; j < sub1.getColumnDimension(); ++j) { + assertEquals((i + 2) + (j + 5) / 1024.0, sub1.getEntry(i, j), 0); + } + } + + RealMatrix sub2 = m.getSubMatrix(10, 12, 3, 70); + for (int i = 0; i < sub2.getRowDimension(); ++i) { + for (int j = 0; j < sub2.getColumnDimension(); ++j) { + assertEquals((i + 10) + (j + 3) / 1024.0, sub2.getEntry(i, j), 0); + } + } + + RealMatrix sub3 = m.getSubMatrix(30, 34, 0, 5); + for (int i = 0; i < sub3.getRowDimension(); ++i) { + for (int j = 0; j < sub3.getColumnDimension(); ++j) { + assertEquals((i + 30) + (j + 0) / 1024.0, sub3.getEntry(i, j), 0); + } + } + + RealMatrix sub4 = m.getSubMatrix(30, 32, 62, 65); + for (int i = 0; i < sub4.getRowDimension(); ++i) { + for (int j = 0; j < sub4.getColumnDimension(); ++j) { + assertEquals((i + 30) + (j + 62) / 1024.0, sub4.getEntry(i, j), 0); + } + } + + } + + //Additional Test for RecursiveLayoutRealMatrixTest.testMultiply + + private double[][] d3 = new double[][] {{1,2,3,4},{5,6,7,8}}; + private double[][] d4 = new double[][] {{1},{2},{3},{4}}; + private double[][] d5 = new double[][] {{30},{70}}; + + public void testMultiply2() { + RealMatrix m3 = new RecursiveLayoutRealMatrix(d3); + RealMatrix m4 = new RecursiveLayoutRealMatrix(d4); + RealMatrix m5 = new RecursiveLayoutRealMatrix(d5); + assertClose(m3.multiply(m4), m5, entryTolerance); + } + + /** test trace */ + public void testTrace() { + RealMatrix m = new RecursiveLayoutRealMatrix(id); + assertEquals("identity trace",3d,m.getTrace(),entryTolerance); + m = new RecursiveLayoutRealMatrix(testData2); + try { + m.getTrace(); + fail("Expecting NonSquareMatrixException"); + } catch (NonSquareMatrixException ex) { + ; + } + } + + /** test scalarAdd */ + public void testScalarAdd() { + RealMatrix m = new RecursiveLayoutRealMatrix(testData); + assertClose(new RecursiveLayoutRealMatrix(testDataPlus2), m.scalarAdd(2d), entryTolerance); + } + + /** test operate */ + public void testOperate() { + RealMatrix m = new RecursiveLayoutRealMatrix(id); + assertClose(testVector, m.operate(testVector), entryTolerance); + assertClose(testVector, m.operate(new RealVectorImpl(testVector)).getData(), entryTolerance); + m = new RecursiveLayoutRealMatrix(bigSingular); + try { + m.operate(testVector); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + ; + } + } + + public void testMultiplyMedium() { + RealMatrix m1 = new RecursiveLayoutRealMatrix( + new double[][] { + { 80, 45, 13, 77, -82 }, + { -90, 33, 98, 80, 74 }, + { 24, -37, 36, -8, -69 }, + { -74, 2, 32, -67, -65 }, + { -29, -81, 44, 54, -65 }, + { 17, 58, -36, -98, 25 }, + { 48, -64, -95, -75, 34 } + }); + RealMatrix m2 = new RecursiveLayoutRealMatrix( + new double[][] { + { 81, 58, 70, 18, 5, -57 }, + { -54, 33, 87, 68, -22, 73 }, + { -78, -5, 34, -7, -3, -31 }, + { -16, -82, -68, 7, 10, -47 }, + { 51, 4, 92, 15, 32, -51 } + }); + RealMatrix m1m2 = m1.multiply(m2); + RealMatrix reference = new RecursiveLayoutRealMatrix( + new double[][]{ + { -2378, -582, -2823, 3718, -2483, -1115 }, + { -14222, -10885, 1271, 1608, 1698, -3033 }, + { -2257, 371, -6119, -3427, -1462, -1290 }, + { -10841, 848, -5342, -2864, -3260, 9836 }, + { -5586, -9263, -17233, -6935, -35, -4847 }, + { 3896, 11216, 13976, 4191, -1263, 7712 }, + { 17688, 7433, 2790, -2838, 2271, -2672 } + }); + assertEquals(0, m1m2.subtract(reference).getNorm(), 0.0); + } + + public void testOperateLarge() { + int testBlockSize = 64; + int p = (7 * testBlockSize) / 2; + int q = (5 * testBlockSize) / 2; + int r = 3 * testBlockSize; + Random random = new Random(111007463902334l); + RealMatrix m1 = createRandomMatrix(random, p, q); + RealMatrix m2 = createRandomMatrix(random, q, r); + RealMatrix m1m2 = m1.multiply(m2); + for (int i = 0; i < r; ++i) { + checkArrays(m1m2.getColumn(i), m1.operate(m2.getColumn(i))); + } + } + + public void testOperatePremultiplyLarge() { + int testBlockSize = 64; + int p = (7 * testBlockSize) / 2; + int q = (5 * testBlockSize) / 2; + int r = 3 * testBlockSize; + Random random = new Random(111007463902334l); + RealMatrix m1 = createRandomMatrix(random, p, q); + RealMatrix m2 = createRandomMatrix(random, q, r); + RealMatrix m1m2 = m1.multiply(m2); + for (int i = 0; i < p; ++i) { + checkArrays(m1m2.getRow(i), m2.preMultiply(m1.getRow(i))); + } + } + + /** test issue MATH-209 */ + public void testMath209() { + RealMatrix a = new RecursiveLayoutRealMatrix(new double[][] { + { 1, 2 }, { 3, 4 }, { 5, 6 } + }); + double[] b = a.operate(new double[] { 1, 1 }); + assertEquals(a.getRowDimension(), b.length); + assertEquals( 3.0, b[0], 1.0e-12); + assertEquals( 7.0, b[1], 1.0e-12); + assertEquals(11.0, b[2], 1.0e-12); + } + + /** test transpose */ + public void testTranspose() { + RealMatrix m = new RecursiveLayoutRealMatrix(testData); + RealMatrix mIT = new LUDecompositionImpl(m).getSolver().getInverse().transpose(); + RealMatrix mTI = new LUDecompositionImpl(m.transpose()).getSolver().getInverse(); + assertClose(mIT, mTI, normTolerance); + m = new RecursiveLayoutRealMatrix(testData2); + RealMatrix mt = new RecursiveLayoutRealMatrix(testData2T); + assertClose(mt, m.transpose(), normTolerance); + } + + /** test preMultiply by vector */ + public void testPremultiplyVector() { + RealMatrix m = new RecursiveLayoutRealMatrix(testData); + assertClose(m.preMultiply(testVector), preMultTest, normTolerance); + assertClose(m.preMultiply(new RealVectorImpl(testVector).getData()), + preMultTest, normTolerance); + m = new RecursiveLayoutRealMatrix(bigSingular); + try { + m.preMultiply(testVector); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException ex) { + ; + } + } + + public void testPremultiply() { + RealMatrix m3 = new RecursiveLayoutRealMatrix(d3); + RealMatrix m4 = new RecursiveLayoutRealMatrix(d4); + RealMatrix m5 = new RecursiveLayoutRealMatrix(d5); + assertClose(m4.preMultiply(m3), m5, entryTolerance); + + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix mInv = new RecursiveLayoutRealMatrix(testDataInv); + RecursiveLayoutRealMatrix identity = new RecursiveLayoutRealMatrix(id); + assertClose(m.preMultiply(mInv), identity, entryTolerance); + assertClose(mInv.preMultiply(m), identity, entryTolerance); + assertClose(m.preMultiply(identity), m, entryTolerance); + assertClose(identity.preMultiply(mInv), mInv, entryTolerance); + try { + m.preMultiply(new RecursiveLayoutRealMatrix(bigSingular)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + ; + } + } + + public void testGetVectors() { + RealMatrix m = new RecursiveLayoutRealMatrix(testData); + assertClose(m.getRow(0), testDataRow1, entryTolerance); + assertClose(m.getColumn(2), testDataCol3, entryTolerance); + try { + m.getRow(10); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + ; + } + try { + m.getColumn(-1); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + ; + } + } + + public void testGetEntry() { + RealMatrix m = new RecursiveLayoutRealMatrix(testData); + assertEquals("get entry",m.getEntry(0,1),2d,entryTolerance); + try { + m.getEntry(10, 4); + fail ("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + /** test examples in user guide */ + public void testExamples() { + // Create a real matrix with two rows and three columns + double[][] matrixData = { {1d,2d,3d}, {2d,5d,3d}}; + RealMatrix m = new RecursiveLayoutRealMatrix(matrixData); + // One more with three rows, two columns + double[][] matrixData2 = { {1d,2d}, {2d,5d}, {1d, 7d}}; + RealMatrix n = new RecursiveLayoutRealMatrix(matrixData2); + // Now multiply m by n + RealMatrix p = m.multiply(n); + assertEquals(2, p.getRowDimension()); + assertEquals(2, p.getColumnDimension()); + // Invert p + RealMatrix pInverse = new LUDecompositionImpl(p).getSolver().getInverse(); + assertEquals(2, pInverse.getRowDimension()); + assertEquals(2, pInverse.getColumnDimension()); + + // Solve example + double[][] coefficientsData = {{2, 3, -2}, {-1, 7, 6}, {4, -3, -5}}; + RealMatrix coefficients = new RecursiveLayoutRealMatrix(coefficientsData); + double[] constants = {1, -2, 1}; + double[] solution = new LUDecompositionImpl(coefficients).getSolver().solve(constants); + assertEquals(2 * solution[0] + 3 * solution[1] -2 * solution[2], constants[0], 1E-12); + assertEquals(-1 * solution[0] + 7 * solution[1] + 6 * solution[2], constants[1], 1E-12); + assertEquals(4 * solution[0] - 3 * solution[1] -5 * solution[2], constants[2], 1E-12); + + } + + // test submatrix accessors + public void testGetSubMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + checkGetSubMatrix(m, subRows23Cols00, 2 , 3 , 0, 0, false); + checkGetSubMatrix(m, subRows00Cols33, 0 , 0 , 3, 3, false); + checkGetSubMatrix(m, subRows01Cols23, 0 , 1 , 2, 3, false); + checkGetSubMatrix(m, subRows02Cols13, new int[] { 0, 2 }, new int[] { 1, 3 }, false); + checkGetSubMatrix(m, subRows03Cols12, new int[] { 0, 3 }, new int[] { 1, 2 }, false); + checkGetSubMatrix(m, subRows03Cols123, new int[] { 0, 3 }, new int[] { 1, 2, 3 }, false); + checkGetSubMatrix(m, subRows20Cols123, new int[] { 2, 0 }, new int[] { 1, 2, 3 }, false); + checkGetSubMatrix(m, subRows31Cols31, new int[] { 3, 1 }, new int[] { 3, 1 }, false); + checkGetSubMatrix(m, subRows31Cols31, new int[] { 3, 1 }, new int[] { 3, 1 }, false); + checkGetSubMatrix(m, null, 1, 0, 2, 4, true); + checkGetSubMatrix(m, null, -1, 1, 2, 2, true); + checkGetSubMatrix(m, null, 1, 0, 2, 2, true); + checkGetSubMatrix(m, null, 1, 0, 2, 4, true); + checkGetSubMatrix(m, null, new int[] {}, new int[] { 0 }, true); + checkGetSubMatrix(m, null, new int[] { 0 }, new int[] { 4 }, true); + } + + private void checkGetSubMatrix(RealMatrix m, double[][] reference, + int startRow, int endRow, int startColumn, int endColumn, + boolean mustFail) { + try { + RealMatrix sub = m.getSubMatrix(startRow, endRow, startColumn, endColumn); + assertEquals(new RecursiveLayoutRealMatrix(reference), sub); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + private void checkGetSubMatrix(RealMatrix m, double[][] reference, + int[] selectedRows, int[] selectedColumns, + boolean mustFail) { + try { + RealMatrix sub = m.getSubMatrix(selectedRows, selectedColumns); + assertEquals(new RecursiveLayoutRealMatrix(reference), sub); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + public void testGetSetMatrixLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + RealMatrix sub = new RecursiveLayoutRealMatrix(n - 4, n - 4).scalarAdd(1); + + m.setSubMatrix(sub.getData(), 2, 2); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if ((i < 2) || (i > n - 3) || (j < 2) || (j > n - 3)) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + assertEquals(sub, m.getSubMatrix(2, n - 3, 2, n - 3)); + + } + + public void testCopySubMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + checkCopy(m, subRows23Cols00, 2 , 3 , 0, 0, false); + checkCopy(m, subRows00Cols33, 0 , 0 , 3, 3, false); + checkCopy(m, subRows01Cols23, 0 , 1 , 2, 3, false); + checkCopy(m, subRows02Cols13, new int[] { 0, 2 }, new int[] { 1, 3 }, false); + checkCopy(m, subRows03Cols12, new int[] { 0, 3 }, new int[] { 1, 2 }, false); + checkCopy(m, subRows03Cols123, new int[] { 0, 3 }, new int[] { 1, 2, 3 }, false); + checkCopy(m, subRows20Cols123, new int[] { 2, 0 }, new int[] { 1, 2, 3 }, false); + checkCopy(m, subRows31Cols31, new int[] { 3, 1 }, new int[] { 3, 1 }, false); + checkCopy(m, subRows31Cols31, new int[] { 3, 1 }, new int[] { 3, 1 }, false); + + checkCopy(m, null, 1, 0, 2, 4, true); + checkCopy(m, null, -1, 1, 2, 2, true); + checkCopy(m, null, 1, 0, 2, 2, true); + checkCopy(m, null, 1, 0, 2, 4, true); + checkCopy(m, null, new int[] {}, new int[] { 0 }, true); + checkCopy(m, null, new int[] { 0 }, new int[] { 4 }, true); + } + + private void checkCopy(RealMatrix m, double[][] reference, + int startRow, int endRow, int startColumn, int endColumn, + boolean mustFail) { + try { + double[][] sub = (reference == null) ? + new double[1][1] : + new double[reference.length][reference[0].length]; + m.copySubMatrix(startRow, endRow, startColumn, endColumn, sub); + assertEquals(new RecursiveLayoutRealMatrix(reference), new RecursiveLayoutRealMatrix(sub)); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + private void checkCopy(RealMatrix m, double[][] reference, + int[] selectedRows, int[] selectedColumns, + boolean mustFail) { + try { + double[][] sub = (reference == null) ? + new double[1][1] : + new double[reference.length][reference[0].length]; + m.copySubMatrix(selectedRows, selectedColumns, sub); + assertEquals(new RecursiveLayoutRealMatrix(reference), new RecursiveLayoutRealMatrix(sub)); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + public void testGetRowMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealMatrix mRow0 = new RecursiveLayoutRealMatrix(subRow0); + RealMatrix mRow3 = new RecursiveLayoutRealMatrix(subRow3); + assertEquals("Row0", mRow0, m.getRowMatrix(0)); + assertEquals("Row3", mRow3, m.getRowMatrix(3)); + try { + m.getRowMatrix(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getRowMatrix(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetRowMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealMatrix mRow3 = new RecursiveLayoutRealMatrix(subRow3); + assertNotSame(mRow3, m.getRowMatrix(0)); + m.setRowMatrix(0, mRow3); + assertEquals(mRow3, m.getRowMatrix(0)); + try { + m.setRowMatrix(-1, mRow3); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setRowMatrix(0, m); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetRowMatrixLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + RealMatrix sub = new RecursiveLayoutRealMatrix(1, n).scalarAdd(1); + + m.setRowMatrix(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + assertEquals(sub, m.getRowMatrix(2)); + + } + + public void testGetColumnMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealMatrix mColumn1 = new RecursiveLayoutRealMatrix(subColumn1); + RealMatrix mColumn3 = new RecursiveLayoutRealMatrix(subColumn3); + assertEquals(mColumn1, m.getColumnMatrix(1)); + assertEquals(mColumn3, m.getColumnMatrix(3)); + try { + m.getColumnMatrix(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getColumnMatrix(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetColumnMatrix() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealMatrix mColumn3 = new RecursiveLayoutRealMatrix(subColumn3); + assertNotSame(mColumn3, m.getColumnMatrix(1)); + m.setColumnMatrix(1, mColumn3); + assertEquals(mColumn3, m.getColumnMatrix(1)); + try { + m.setColumnMatrix(-1, mColumn3); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setColumnMatrix(0, m); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetColumnMatrixLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + RealMatrix sub = new RecursiveLayoutRealMatrix(n, 1).scalarAdd(1); + + m.setColumnMatrix(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + assertEquals(sub, m.getColumnMatrix(2)); + + } + + public void testGetRowVector() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealVector mRow0 = new RealVectorImpl(subRow0[0]); + RealVector mRow3 = new RealVectorImpl(subRow3[0]); + assertEquals(mRow0, m.getRowVector(0)); + assertEquals(mRow3, m.getRowVector(3)); + try { + m.getRowVector(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getRowVector(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetRowVector() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealVector mRow3 = new RealVectorImpl(subRow3[0]); + assertNotSame(mRow3, m.getRowMatrix(0)); + m.setRowVector(0, mRow3); + assertEquals(mRow3, m.getRowVector(0)); + try { + m.setRowVector(-1, mRow3); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setRowVector(0, new RealVectorImpl(5)); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetRowVectorLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + RealVector sub = new RealVectorImpl(n, 1.0); + + m.setRowVector(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + assertEquals(sub, m.getRowVector(2)); + + } + + public void testGetColumnVector() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealVector mColumn1 = columnToVector(subColumn1); + RealVector mColumn3 = columnToVector(subColumn3); + assertEquals(mColumn1, m.getColumnVector(1)); + assertEquals(mColumn3, m.getColumnVector(3)); + try { + m.getColumnVector(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getColumnVector(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetColumnVector() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + RealVector mColumn3 = columnToVector(subColumn3); + assertNotSame(mColumn3, m.getColumnVector(1)); + m.setColumnVector(1, mColumn3); + assertEquals(mColumn3, m.getColumnVector(1)); + try { + m.setColumnVector(-1, mColumn3); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setColumnVector(0, new RealVectorImpl(5)); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetColumnVectorLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + RealVector sub = new RealVectorImpl(n, 1.0); + + m.setColumnVector(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + assertEquals(sub, m.getColumnVector(2)); + + } + + private RealVector columnToVector(double[][] column) { + double[] data = new double[column.length]; + for (int i = 0; i < data.length; ++i) { + data[i] = column[i][0]; + } + return new RealVectorImpl(data, false); + } + + public void testGetRow() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + checkArrays(subRow0[0], m.getRow(0)); + checkArrays(subRow3[0], m.getRow(3)); + try { + m.getRow(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getRow(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetRow() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + assertTrue(subRow3[0][0] != m.getRow(0)[0]); + m.setRow(0, subRow3[0]); + checkArrays(subRow3[0], m.getRow(0)); + try { + m.setRow(-1, subRow3[0]); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setRow(0, new double[5]); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetRowLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + double[] sub = new double[n]; + Arrays.fill(sub, 1.0); + + m.setRow(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + checkArrays(sub, m.getRow(2)); + + } + + public void testGetColumn() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + double[] mColumn1 = columnToArray(subColumn1); + double[] mColumn3 = columnToArray(subColumn3); + checkArrays(mColumn1, m.getColumn(1)); + checkArrays(mColumn3, m.getColumn(3)); + try { + m.getColumn(-1); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.getColumn(4); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + } + + public void testSetColumn() { + RealMatrix m = new RecursiveLayoutRealMatrix(subTestData); + double[] mColumn3 = columnToArray(subColumn3); + assertTrue(mColumn3[0] != m.getColumn(1)[0]); + m.setColumn(1, mColumn3); + checkArrays(mColumn3, m.getColumn(1)); + try { + m.setColumn(-1, mColumn3); + fail("Expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // expected + } + try { + m.setColumn(0, new double[5]); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetColumnLarge() { + int n = 3 * 64; + RealMatrix m = new RecursiveLayoutRealMatrix(n, n); + double[] sub = new double[n]; + Arrays.fill(sub, 1.0); + + m.setColumn(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(0.0, m.getEntry(i, j), 0.0); + } else { + assertEquals(1.0, m.getEntry(i, j), 0.0); + } + } + } + checkArrays(sub, m.getColumn(2)); + + } + + private double[] columnToArray(double[][] column) { + double[] data = new double[column.length]; + for (int i = 0; i < data.length; ++i) { + data[i] = column[i][0]; + } + return data; + } + + private void checkArrays(double[] expected, double[] actual) { + assertEquals(expected.length, actual.length); + for (int i = 0; i < expected.length; ++i) { + assertEquals(expected[i], actual[i], 1.0e-9 * Math.abs(expected[i])); + } + } + + public void testEqualsAndHashCode() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + RecursiveLayoutRealMatrix m1 = (RecursiveLayoutRealMatrix) m.copy(); + RecursiveLayoutRealMatrix mt = (RecursiveLayoutRealMatrix) m.transpose(); + assertTrue(m.hashCode() != mt.hashCode()); + assertEquals(m.hashCode(), m1.hashCode()); + assertEquals(m, m); + assertEquals(m, m1); + assertFalse(m.equals(null)); + assertFalse(m.equals(mt)); + assertFalse(m.equals(new RecursiveLayoutRealMatrix(bigSingular))); + } + + public void testToString() { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + assertEquals("RecursiveLayoutRealMatrix{{1.0,2.0,3.0},{2.0,5.0,3.0},{1.0,0.0,8.0}}", + m.toString()); + } + + public void testSetSubMatrix() throws Exception { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(testData); + m.setSubMatrix(detData2,1,1); + RealMatrix expected = new RecursiveLayoutRealMatrix + (new double[][] {{1.0,2.0,3.0},{2.0,1.0,3.0},{1.0,2.0,4.0}}); + assertEquals(expected, m); + + m.setSubMatrix(detData2,0,0); + expected = new RecursiveLayoutRealMatrix + (new double[][] {{1.0,3.0,3.0},{2.0,4.0,3.0},{1.0,2.0,4.0}}); + assertEquals(expected, m); + + m.setSubMatrix(testDataPlus2,0,0); + expected = new RecursiveLayoutRealMatrix + (new double[][] {{3.0,4.0,5.0},{4.0,7.0,5.0},{3.0,2.0,10.0}}); + assertEquals(expected, m); + + // javadoc example + RecursiveLayoutRealMatrix matrix = new RecursiveLayoutRealMatrix + (new double[][] {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 0, 1 , 2}}); + matrix.setSubMatrix(new double[][] {{3, 4}, {5, 6}}, 1, 1); + expected = new RecursiveLayoutRealMatrix + (new double[][] {{1, 2, 3, 4}, {5, 3, 4, 8}, {9, 5 ,6, 2}}); + assertEquals(expected, matrix); + + // dimension overflow + try { + m.setSubMatrix(testData,1,1); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException e) { + // expected + } + // dimension underflow + try { + m.setSubMatrix(testData,-1,1); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException e) { + // expected + } + try { + m.setSubMatrix(testData,1,-1); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException e) { + // expected + } + + // null + try { + m.setSubMatrix(null,1,1); + fail("expecting NullPointerException"); + } catch (NullPointerException e) { + // expected + } + + // ragged + try { + m.setSubMatrix(new double[][] {{1}, {2, 3}}, 0, 0); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException e) { + // expected + } + + // empty + try { + m.setSubMatrix(new double[][] {{}}, 0, 0); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException e) { + // expected + } + + } + + public void testWalk() { + int rows = 150; + int columns = 75; + + RealMatrix m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInRowOrder(new SetVisitor()); + GetVisitor getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInRowOrder(new SetVisitor(), 1, rows - 2, 1, columns - 2); + getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor, 1, rows - 2, 1, columns - 2); + assertEquals((rows - 2) * (columns - 2), getVisitor.getCount()); + for (int i = 0; i < rows; ++i) { + assertEquals(0.0, m.getEntry(i, 0), 0); + assertEquals(0.0, m.getEntry(i, columns - 1), 0); + } + for (int j = 0; j < columns; ++j) { + assertEquals(0.0, m.getEntry(0, j), 0); + assertEquals(0.0, m.getEntry(rows - 1, j), 0); + } + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInColumnOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInColumnOrder(new SetVisitor(), 1, rows - 2, 1, columns - 2); + getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor, 1, rows - 2, 1, columns - 2); + assertEquals((rows - 2) * (columns - 2), getVisitor.getCount()); + for (int i = 0; i < rows; ++i) { + assertEquals(0.0, m.getEntry(i, 0), 0); + assertEquals(0.0, m.getEntry(i, columns - 1), 0); + } + for (int j = 0; j < columns; ++j) { + assertEquals(0.0, m.getEntry(0, j), 0); + assertEquals(0.0, m.getEntry(rows - 1, j), 0); + } + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInOptimizedOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInRowOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInOptimizedOrder(new SetVisitor(), 1, rows - 2, 1, columns - 2); + getVisitor = new GetVisitor(); + m.walkInRowOrder(getVisitor, 1, rows - 2, 1, columns - 2); + assertEquals((rows - 2) * (columns - 2), getVisitor.getCount()); + for (int i = 0; i < rows; ++i) { + assertEquals(0.0, m.getEntry(i, 0), 0); + assertEquals(0.0, m.getEntry(i, columns - 1), 0); + } + for (int j = 0; j < columns; ++j) { + assertEquals(0.0, m.getEntry(0, j), 0); + assertEquals(0.0, m.getEntry(rows - 1, j), 0); + } + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInOptimizedOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInColumnOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new RecursiveLayoutRealMatrix(rows, columns); + m.walkInOptimizedOrder(new SetVisitor(), 1, rows - 2, 1, columns - 2); + getVisitor = new GetVisitor(); + m.walkInColumnOrder(getVisitor, 1, rows - 2, 1, columns - 2); + assertEquals((rows - 2) * (columns - 2), getVisitor.getCount()); + for (int i = 0; i < rows; ++i) { + assertEquals(0.0, m.getEntry(i, 0), 0); + assertEquals(0.0, m.getEntry(i, columns - 1), 0); + } + for (int j = 0; j < columns; ++j) { + assertEquals(0.0, m.getEntry(0, j), 0); + assertEquals(0.0, m.getEntry(rows - 1, j), 0); + } + + } + + private static class SetVisitor extends DefaultRealMatrixChangingVisitor { + private static final long serialVersionUID = 1773444180892369386L; + public double visit(int i, int j, double value) { + return i + j / 1024.0; + } + } + + private static class GetVisitor extends DefaultRealMatrixPreservingVisitor { + private static final long serialVersionUID = -7745543227178932689L; + private int count = 0; + public void visit(int i, int j, double value) { + ++count; + assertEquals(i + j / 1024.0, value, 0.0); + } + public int getCount() { + return count; + } + }; + + //--------------- -----------------Protected methods + + /** verifies that two matrices are close (1-norm) */ + protected void assertClose(RealMatrix m, RealMatrix n, double tolerance) { + assertTrue(m.subtract(n).getNorm() < tolerance); + } + + /** verifies that two vectors are close (sup norm) */ + protected void assertClose(double[] m, double[] n, double tolerance) { + if (m.length != n.length) { + fail("vectors not same length"); + } + for (int i = 0; i < m.length; i++) { + assertEquals(m[i], n[i], tolerance); + } + } + + private RecursiveLayoutRealMatrix createRandomMatrix(Random r, int rows, int columns) { + RecursiveLayoutRealMatrix m = new RecursiveLayoutRealMatrix(rows, columns); + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < columns; ++j) { + m.setEntry(i, j, 200 * r.nextDouble() - 100); + } + } + return m; + } + +} +