From b70bdcde87a60421221d3b2bb3d48e7316e84b74 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Mon, 20 Apr 2009 20:19:49 +0000 Subject: [PATCH] added DenseFieldMatrix git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@766852 13f79535-47bb-0310-9956-ffa450edef68 --- .../commons/math/linear/DenseFieldMatrix.java | 1613 +++++++++++++++++ .../commons/math/linear/MatrixUtils.java | 186 +- .../math/linear/DenseFieldMatrixTest.java | 1294 +++++++++++++ 3 files changed, 3075 insertions(+), 18 deletions(-) create mode 100644 src/java/org/apache/commons/math/linear/DenseFieldMatrix.java create mode 100644 src/test/org/apache/commons/math/linear/DenseFieldMatrixTest.java diff --git a/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java b/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java new file mode 100644 index 000000000..de5e309d3 --- /dev/null +++ b/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java @@ -0,0 +1,1613 @@ +/* + * 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 org.apache.commons.math.Field; +import org.apache.commons.math.FieldElement; +import org.apache.commons.math.MathRuntimeException; + +/** + * Cache-friendly implementation of FieldMatrix using a flat arrays to store + * square blocks of the matrix. + *

+ * This implementation is specially designed to be cache-friendly. Square blocks are + * stored as small arrays and allow efficient traversal of data both in row major direction + * and columns major direction, one block at a time. This greatly increases performances + * for algorithms that use crossed directions loops like multiplication or transposition. + *

+ *

+ * The size of square blocks is a static parameter. It may be tuned according to the cache + * size of the target computer processor. As a rule of thumbs, it should be the largest + * value that allows three blocks to be simultaneously cached (this is necessary for example + * for matrix multiplication). The default value is to use 36x36 blocks. + *

+ *

+ * The regular blocks represent {@link #BLOCK_SIZE} x {@link #BLOCK_SIZE} squares. Blocks + * at right hand side and bottom side which may be smaller to fit matrix dimensions. The square + * blocks are flattened in row major order in single dimension arrays which are therefore + * {@link #BLOCK_SIZE}2 elements long for regular blocks. The blocks are themselves + * organized in row major order. + *

+ *

+ * As an example, for a block size of 36x36, a 100x60 matrix would be stored in 6 blocks. + * Block 0 would be a Field[1296] array holding the upper left 36x36 square, block 1 would be + * a Field[1296] array holding the upper center 36x36 square, block 2 would be a Field[1008] + * array holding the upper right 36x28 rectangle, block 3 would be a Field[864] array holding + * the lower left 24x36 rectangle, block 4 would be a Field[864] array holding the lower center + * 24x36 rectangle and block 5 would be a Field[672] array holding the lower right 24x28 + * rectangle. + *

+ *

+ * The layout complexity overhead versus simple mapping of matrices to java + * arrays is negligible for small matrices (about 1%). The gain from cache efficiency leads + * to up to 3-fold improvements for matrices of moderate to large size. + *

+ * @version $Revision$ $Date$ + * @since 2.0 + */ +public class DenseFieldMatrix> extends AbstractFieldMatrix { + + /** Serializable version identifier */ + private static final long serialVersionUID = -4602336630143123183L; + + /** Block size. */ + public static final int BLOCK_SIZE = 36; + + /** Blocks of matrix entries. */ + private final T blocks[][]; + + /** Number of rows of the matrix. */ + private final int rows; + + /** Number of columns of the matrix. */ + private final int columns; + + /** Number of block rows of the matrix. */ + private final int blockRows; + + /** Number of block columns of the matrix. */ + private final int blockColumns; + + /** + * Create a new matrix with the supplied row and column dimensions. + * + * @param field field to which the elements belong + * @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 DenseFieldMatrix(final Field field, final int rows, final int columns) + throws IllegalArgumentException { + + super(field, rows, columns); + this.rows = rows; + this.columns = columns; + + // number of blocks + blockRows = (rows + BLOCK_SIZE - 1) / BLOCK_SIZE; + blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; + + // allocate storage blocks, taking care of smaller ones at right and bottom + blocks = createBlocksLayout(field, rows, columns); + + } + + /** + * Create a new dense matrix copying entries from raw layout data. + *

The input array must already be in raw layout.

+ *

Calling this constructor is equivalent to call: + * 

matrix = new DenseFieldMatrix(getField(), rawData.length, rawData[0].length,
+     *                                   toBlocksLayout(rawData), false);
+ *

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

The input array must already be in blocks layout.

+ * @param rows the number of rows in the new matrix + * @param columns the number of columns in the new matrix + * @param blockData data for new matrix + * @param copyArray if true, the input array will be copied, otherwise + * it will be referenced + * + * @exception IllegalArgumentException if blockData shape is + * inconsistent with block layout + * @see #createBlocksLayout(int, int) + * @see #toBlocksLayout(T[][]) + * @see #DenseFieldMatrix(T[][]) + */ + public DenseFieldMatrix(final int rows, final int columns, + final T[][] blockData, final boolean copyArray) + throws IllegalArgumentException { + + super(extractField(blockData), rows, columns); + this.rows = rows; + this.columns = columns; + + // number of blocks + blockRows = (rows + BLOCK_SIZE - 1) / BLOCK_SIZE; + blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; + + if (copyArray) { + // allocate storage blocks, taking care of smaller ones at right and bottom + blocks = buildArray(getField(), blockRows * blockColumns, -1); + } else { + // reference existing array + blocks = blockData; + } + + int index = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++index) { + if (blockData[index].length != iHeight * blockWidth(jBlock)) { + throw MathRuntimeException.createIllegalArgumentException( + "wrong array shape (block length = {0}, expected {1})", + blockData[index].length, iHeight * blockWidth(jBlock)); + } + if (copyArray) { + blocks[index] = blockData[index].clone(); + } + } + } + + } + + /** + * Convert a data array from raw layout to blocks layout. + *

+ * Raw layout is the straightforward layout where element at row i and + * column j is in array element rawData[i][j]. Blocks layout + * is the layout used in {@link DenseFieldMatrix} instances, where the matrix + * is split in square blocks (except at right and bottom side where blocks may + * be rectangular to fit matrix size) and each block is stored in a flattened + * one-dimensional array. + *

+ *

+ * This method creates an array in blocks layout from an input array in raw layout. + * It can be used to provide the array argument of the {@link + * DenseFieldMatrix#DenseFieldMatrix(int, int, T[][], boolean)} constructor. + *

+ * @param rawData data array in raw layout + * @return a new data array containing the same entries but in blocks layout + * @exception IllegalArgumentException if rawData is not rectangular + * (not all rows have the same length) + * @see #createBlocksLayout(int, int) + * @see #DenseFieldMatrix(int, int, T[][], boolean) + */ + public static > T[][] toBlocksLayout(final T[][] rawData) + throws IllegalArgumentException { + + final int rows = rawData.length; + final int columns = rawData[0].length; + final int blockRows = (rows + BLOCK_SIZE - 1) / BLOCK_SIZE; + final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; + + // 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}", + columns, length); + } + } + + // convert array + final Field field = extractField(rawData); + final T[][] blocks = buildArray(field, blockRows * blockColumns, -1); + for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + final int iHeight = pEnd - pStart; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final int jWidth = qEnd - qStart; + + // allocate new block + final T[] block = buildArray(field, iHeight * jWidth); + blocks[blockIndex] = block; + + // copy data + for (int p = pStart, index = 0; p < pEnd; ++p, index += jWidth) { + System.arraycopy(rawData[p], qStart, block, index, jWidth); + } + + } + } + + return blocks; + + } + + /** + * Create a data array in blocks layout. + *

+ * This method can be used to create the array argument of the {@link + * DenseFieldMatrix#DenseFieldMatrix(int, int, T[][], boolean)} constructor. + *

+ * @param rows the number of rows in the new matrix + * @param columns the number of columns in the new matrix + * @return a new data array in blocks layout + * @see #toBlocksLayout(T[][]) + * @see #DenseFieldMatrix(int, int, T[][], boolean) + */ + public static > T[][] createBlocksLayout(final Field field, + final int rows, final int columns) { + + final int blockRows = (rows + BLOCK_SIZE - 1) / BLOCK_SIZE; + final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; + + final T[][] blocks = buildArray(field, blockRows * blockColumns, -1); + for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + final int iHeight = pEnd - pStart; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final int jWidth = qEnd - qStart; + blocks[blockIndex] = buildArray(field, iHeight * jWidth); + } + } + + return blocks; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix createMatrix(final int rowDimension, final int columnDimension) + throws IllegalArgumentException { + return new DenseFieldMatrix(getField(), rowDimension, columnDimension); + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix copy() { + + // create an empty matrix + DenseFieldMatrix copied = new DenseFieldMatrix(getField(), rows, columns); + + // copy the blocks + for (int i = 0; i < blocks.length; ++i) { + System.arraycopy(blocks[i], 0, copied.blocks[i], 0, blocks[i].length); + } + + return copied; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix add(final FieldMatrix m) + throws IllegalArgumentException { + try { + return add((DenseFieldMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkAdditionCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform addition block-wise, to ensure good cache behavior + int blockIndex = 0; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { + + // perform addition on the current block + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + for (int p = pStart, k = 0; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q, ++k) { + outBlock[k] = tBlock[k].add(m.getEntry(p, q)); + } + } + + // go to next block + ++blockIndex; + + } + } + + 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 DenseFieldMatrix add(final DenseFieldMatrix m) + throws IllegalArgumentException { + + // safety check + checkAdditionCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform addition block-wise, to ensure good cache behavior + for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) { + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + final T[] mBlock = m.blocks[blockIndex]; + for (int k = 0; k < outBlock.length; ++k) { + outBlock[k] = tBlock[k].add(mBlock[k]); + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix subtract(final FieldMatrix m) + throws IllegalArgumentException { + try { + return subtract((DenseFieldMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkSubtractionCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform subtraction block-wise, to ensure good cache behavior + int blockIndex = 0; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { + + // perform subtraction on the current block + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + for (int p = pStart, k = 0; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q, ++k) { + outBlock[k] = tBlock[k].subtract(m.getEntry(p, q)); + } + } + + // go to next block + ++blockIndex; + + } + } + + 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 DenseFieldMatrix subtract(final DenseFieldMatrix m) + throws IllegalArgumentException { + + // safety check + checkSubtractionCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform subtraction block-wise, to ensure good cache behavior + for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) { + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + final T[] mBlock = m.blocks[blockIndex]; + for (int k = 0; k < outBlock.length; ++k) { + outBlock[k] = tBlock[k].subtract(mBlock[k]); + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix scalarAdd(final T d) + throws IllegalArgumentException { + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform subtraction block-wise, to ensure good cache behavior + for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) { + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + for (int k = 0; k < outBlock.length; ++k) { + outBlock[k] = tBlock[k].add(d); + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix scalarMultiply(final T d) + throws IllegalArgumentException { + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, columns); + + // perform subtraction block-wise, to ensure good cache behavior + for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) { + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[blockIndex]; + for (int k = 0; k < outBlock.length; ++k) { + outBlock[k] = tBlock[k].multiply(d); + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix multiply(final FieldMatrix m) + throws IllegalArgumentException { + try { + return multiply((DenseFieldMatrix) m); + } catch (ClassCastException cce) { + + // safety check + checkMultiplicationCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, m.getColumnDimension()); + final T zero = getField().getZero(); + + // perform multiplication block-wise, to ensure good cache behavior + int blockIndex = 0; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { + + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { + + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, m.getColumnDimension()); + + // select current block + final T[] outBlock = out.blocks[blockIndex]; + + // perform multiplication on current block + for (int kBlock = 0; kBlock < blockColumns; ++kBlock) { + final int kWidth = blockWidth(kBlock); + final T[] tBlock = blocks[iBlock * blockColumns + kBlock]; + final int rStart = kBlock * BLOCK_SIZE; + for (int p = pStart, k = 0; p < pEnd; ++p) { + final int lStart = (p - pStart) * kWidth; + final int lEnd = lStart + kWidth; + for (int q = qStart; q < qEnd; ++q) { + T sum = zero; + for (int l = lStart, r = rStart; l < lEnd; ++l, ++r) { + sum = sum.add(tBlock[l].multiply(m.getEntry(r, q))); + } + outBlock[k] = outBlock[k].add(sum); + ++k; + } + } + } + + // go to next block + ++blockIndex; + + } + } + + return out; + + } + } + + /** + * Returns the result of postmultiplying this by m. + * + * @param m matrix to postmultiply by + * @return this * m + * @throws IllegalArgumentException + * if columnDimension(this) != rowDimension(m) + */ + public DenseFieldMatrix multiply(DenseFieldMatrix m) throws IllegalArgumentException { + + // safety check + checkMultiplicationCompatible(m); + + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, m.columns); + final T zero = getField().getZero(); + + // perform multiplication block-wise, to ensure good cache behavior + int blockIndex = 0; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { + + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { + final int jWidth = out.blockWidth(jBlock); + final int jWidth2 = jWidth + jWidth; + final int jWidth3 = jWidth2 + jWidth; + final int jWidth4 = jWidth3 + jWidth; + + // select current block + final T[] outBlock = out.blocks[blockIndex]; + + // perform multiplication on current block + for (int kBlock = 0; kBlock < blockColumns; ++kBlock) { + final int kWidth = blockWidth(kBlock); + final T[] tBlock = blocks[iBlock * blockColumns + kBlock]; + final T[] mBlock = m.blocks[kBlock * m.blockColumns + jBlock]; + for (int p = pStart, k = 0; p < pEnd; ++p) { + final int lStart = (p - pStart) * kWidth; + final int lEnd = lStart + kWidth; + for (int nStart = 0; nStart < jWidth; ++nStart) { + T sum = zero; + int l = lStart; + int n = nStart; + while (l < lEnd - 3) { + sum = sum. + add(tBlock[l].multiply(mBlock[n])). + add(tBlock[l + 1].multiply(mBlock[n + jWidth])). + add(tBlock[l + 2].multiply(mBlock[n + jWidth2])). + add(tBlock[l + 3].multiply(mBlock[n + jWidth3])); + l += 4; + n += jWidth4; + } + while (l < lEnd) { + sum = sum.add(tBlock[l++].multiply(mBlock[n])); + n += jWidth; + } + outBlock[k] = outBlock[k].add(sum); + ++k; + } + } + } + + // go to next block + ++blockIndex; + + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public T[][] getData() { + + final T[][] data = buildArray(getField(), getRowDimension(), getColumnDimension()); + final int lastColumns = columns - (blockColumns - 1) * BLOCK_SIZE; + + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + int regularPos = 0; + int lastPos = 0; + for (int p = pStart; p < pEnd; ++p) { + final T[] dataP = data[p]; + int blockIndex = iBlock * blockColumns; + int dataPos = 0; + for (int jBlock = 0; jBlock < blockColumns - 1; ++jBlock) { + System.arraycopy(blocks[blockIndex++], regularPos, dataP, dataPos, BLOCK_SIZE); + dataPos += BLOCK_SIZE; + } + System.arraycopy(blocks[blockIndex], lastPos, dataP, dataPos, lastColumns); + regularPos += BLOCK_SIZE; + lastPos += lastColumns; + } + } + + return data; + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix 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 DenseFieldMatrix out = + new DenseFieldMatrix(getField(), endRow - startRow + 1, endColumn - startColumn + 1); + + // compute blocks shifts + final int blockStartRow = startRow / BLOCK_SIZE; + final int rowsShift = startRow % BLOCK_SIZE; + final int blockStartColumn = startColumn / BLOCK_SIZE; + final int columnsShift = startColumn % BLOCK_SIZE; + + // perform extraction block-wise, to ensure good cache behavior + for (int iBlock = 0, pBlock = blockStartRow; iBlock < out.blockRows; ++iBlock, ++pBlock) { + final int iHeight = out.blockHeight(iBlock); + for (int jBlock = 0, qBlock = blockStartColumn; jBlock < out.blockColumns; ++jBlock, ++qBlock) { + final int jWidth = out.blockWidth(jBlock); + + // handle one block of the output matrix + final int outIndex = iBlock * out.blockColumns + jBlock; + final T[] outBlock = out.blocks[outIndex]; + final int index = pBlock * blockColumns + qBlock; + final int width = blockWidth(qBlock); + + final int heightExcess = iHeight + rowsShift - BLOCK_SIZE; + final int widthExcess = jWidth + columnsShift - BLOCK_SIZE; + if (heightExcess > 0) { + // the submatrix block spans on two blocks rows from the original matrix + if (widthExcess > 0) { + // the submatrix block spans on two blocks columns from the original matrix + final int width2 = blockWidth(qBlock + 1); + copyBlockPart(blocks[index], width, + rowsShift, BLOCK_SIZE, + columnsShift, BLOCK_SIZE, + outBlock, jWidth, 0, 0); + copyBlockPart(blocks[index + 1], width2, + rowsShift, BLOCK_SIZE, + 0, widthExcess, + outBlock, jWidth, 0, jWidth - widthExcess); + copyBlockPart(blocks[index + blockColumns], width, + 0, heightExcess, + columnsShift, BLOCK_SIZE, + outBlock, jWidth, iHeight - heightExcess, 0); + copyBlockPart(blocks[index + blockColumns + 1], width2, + 0, heightExcess, + 0, widthExcess, + outBlock, jWidth, iHeight - heightExcess, jWidth - widthExcess); + } else { + // the submatrix block spans on one block column from the original matrix + copyBlockPart(blocks[index], width, + rowsShift, BLOCK_SIZE, + columnsShift, jWidth + columnsShift, + outBlock, jWidth, 0, 0); + copyBlockPart(blocks[index + blockColumns], width, + 0, heightExcess, + columnsShift, jWidth + columnsShift, + outBlock, jWidth, iHeight - heightExcess, 0); + } + } else { + // the submatrix block spans on one block row from the original matrix + if (widthExcess > 0) { + // the submatrix block spans on two blocks columns from the original matrix + final int width2 = blockWidth(qBlock + 1); + copyBlockPart(blocks[index], width, + rowsShift, iHeight + rowsShift, + columnsShift, BLOCK_SIZE, + outBlock, jWidth, 0, 0); + copyBlockPart(blocks[index + 1], width2, + rowsShift, iHeight + rowsShift, + 0, widthExcess, + outBlock, jWidth, 0, jWidth - widthExcess); + } else { + // the submatrix block spans on one block column from the original matrix + copyBlockPart(blocks[index], width, + rowsShift, iHeight + rowsShift, + columnsShift, jWidth + columnsShift, + outBlock, jWidth, 0, 0); + } + } + + } + } + + return out; + + } + + /** + * Copy a part of a block into another one + *

This method can be called only when the specified part fits in both + * blocks, no verification is done here.

+ * @param srcBlock source block + * @param srcWidth source block width ({@link #BLOCK_SIZE} or smaller) + * @param srcStartRow start row in the source block + * @param srcEndRow end row (exclusive) in the source block + * @param srcStartColumn start column in the source block + * @param srcEndColumn end column (exclusive) in the source block + * @param dstBlock destination block + * @param dstWidth destination block width ({@link #BLOCK_SIZE} or smaller) + * @param dstStartRow start row in the destination block + * @param dstStartColumn start column in the destination block + */ + private void copyBlockPart(final T[] srcBlock, final int srcWidth, + final int srcStartRow, final int srcEndRow, + final int srcStartColumn, final int srcEndColumn, + final T[] dstBlock, final int dstWidth, + final int dstStartRow, final int dstStartColumn) { + final int length = srcEndColumn - srcStartColumn; + int srcPos = srcStartRow * srcWidth + srcStartColumn; + int dstPos = dstStartRow * dstWidth + dstStartColumn; + for (int srcRow = srcStartRow; srcRow < srcEndRow; ++srcRow) { + System.arraycopy(srcBlock, srcPos, dstBlock, dstPos, length); + srcPos += srcWidth; + dstPos += dstWidth; + } + } + + /** {@inheritDoc} */ + @Override + public void setSubMatrix(final T[][] 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"); + } + final int endRow = row + subMatrix.length - 1; + final int endColumn = column + refLength - 1; + checkSubMatrixIndex(row, endRow, column, endColumn); + for (final T[] subRow : subMatrix) { + if (subRow.length != refLength) { + throw MathRuntimeException.createIllegalArgumentException( + "some rows have length {0} while others have length {1}", + refLength, subRow.length); + } + } + + // compute blocks bounds + final int blockStartRow = row / BLOCK_SIZE; + final int blockEndRow = (endRow + BLOCK_SIZE) / BLOCK_SIZE; + final int blockStartColumn = column / BLOCK_SIZE; + final int blockEndColumn = (endColumn + BLOCK_SIZE) / BLOCK_SIZE; + + // perform copy block-wise, to ensure good cache behavior + for (int iBlock = blockStartRow; iBlock < blockEndRow; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final int firstRow = iBlock * BLOCK_SIZE; + final int iStart = Math.max(row, firstRow); + final int iEnd = Math.min(endRow + 1, firstRow + iHeight); + + for (int jBlock = blockStartColumn; jBlock < blockEndColumn; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int firstColumn = jBlock * BLOCK_SIZE; + final int jStart = Math.max(column, firstColumn); + final int jEnd = Math.min(endColumn + 1, firstColumn + jWidth); + final int jLength = jEnd - jStart; + + // handle one block, row by row + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = iStart; i < iEnd; ++i) { + System.arraycopy(subMatrix[i - row], jStart - column, + block, (i - firstRow) * jWidth + (jStart - firstColumn), + jLength); + } + + } + } + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix getRowMatrix(final int row) + throws MatrixIndexException { + + checkRowIndex(row); + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), 1, columns); + + // perform copy block-wise, to ensure good cache behavior + final int iBlock = row / BLOCK_SIZE; + final int iRow = row - iBlock * BLOCK_SIZE; + int outBlockIndex = 0; + int outIndex = 0; + T[] outBlock = out.blocks[outBlockIndex]; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + final int available = outBlock.length - outIndex; + if (jWidth > available) { + System.arraycopy(block, iRow * jWidth, outBlock, outIndex, available); + outBlock = out.blocks[++outBlockIndex]; + System.arraycopy(block, iRow * jWidth, outBlock, 0, jWidth - available); + outIndex = jWidth - available; + } else { + System.arraycopy(block, iRow * jWidth, outBlock, outIndex, jWidth); + outIndex += jWidth; + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public void setRowMatrix(final int row, final FieldMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + try { + setRowMatrix(row, (DenseFieldMatrix) 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 DenseFieldMatrix 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}", + matrix.getRowDimension(), matrix.getColumnDimension(), + 1, nCols); + } + + // perform copy block-wise, to ensure good cache behavior + final int iBlock = row / BLOCK_SIZE; + final int iRow = row - iBlock * BLOCK_SIZE; + int mBlockIndex = 0; + int mIndex = 0; + T[] mBlock = matrix.blocks[mBlockIndex]; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + final int available = mBlock.length - mIndex; + if (jWidth > available) { + System.arraycopy(mBlock, mIndex, block, iRow * jWidth, available); + mBlock = matrix.blocks[++mBlockIndex]; + System.arraycopy(mBlock, 0, block, iRow * jWidth, jWidth - available); + mIndex = jWidth - available; + } else { + System.arraycopy(mBlock, mIndex, block, iRow * jWidth, jWidth); + mIndex += jWidth; + } + } + + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix getColumnMatrix(final int column) + throws MatrixIndexException { + + checkColumnIndex(column); + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), rows, 1); + + // perform copy block-wise, to ensure good cache behavior + final int jBlock = column / BLOCK_SIZE; + final int jColumn = column - jBlock * BLOCK_SIZE; + final int jWidth = blockWidth(jBlock); + int outBlockIndex = 0; + int outIndex = 0; + T[] outBlock = out.blocks[outBlockIndex]; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = 0; i < iHeight; ++i) { + if (outIndex >= outBlock.length) { + outBlock = out.blocks[++outBlockIndex]; + outIndex = 0; + } + outBlock[outIndex++] = block[i * jWidth + jColumn]; + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public void setColumnMatrix(final int column, final FieldMatrix matrix) + throws MatrixIndexException, InvalidMatrixException { + try { + setColumnMatrix(column, (DenseFieldMatrix) 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 DenseFieldMatrix 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}", + matrix.getRowDimension(), matrix.getColumnDimension(), + nRows, 1); + } + + // perform copy block-wise, to ensure good cache behavior + final int jBlock = column / BLOCK_SIZE; + final int jColumn = column - jBlock * BLOCK_SIZE; + final int jWidth = blockWidth(jBlock); + int mBlockIndex = 0; + int mIndex = 0; + T[] mBlock = matrix.blocks[mBlockIndex]; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = 0; i < iHeight; ++i) { + if (mIndex >= mBlock.length) { + mBlock = matrix.blocks[++mBlockIndex]; + mIndex = 0; + } + block[i * jWidth + jColumn] = mBlock[mIndex++]; + } + } + + } + + /** {@inheritDoc} */ + @Override + public FieldVector getRowVector(final int row) + throws MatrixIndexException { + + checkRowIndex(row); + final T[] outData = buildArray(getField(), columns); + + // perform copy block-wise, to ensure good cache behavior + final int iBlock = row / BLOCK_SIZE; + final int iRow = row - iBlock * BLOCK_SIZE; + int outIndex = 0; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + System.arraycopy(block, iRow * jWidth, outData, outIndex, jWidth); + outIndex += jWidth; + } + + return new FieldVectorImpl(outData, false); + + } + + /** {@inheritDoc} */ + @Override + public void setRowVector(final int row, final FieldVector vector) + throws MatrixIndexException, InvalidMatrixException { + try { + setRow(row, ((FieldVectorImpl) vector).getDataRef()); + } catch (ClassCastException cce) { + super.setRowVector(row, vector); + } + } + + /** {@inheritDoc} */ + @Override + public FieldVector getColumnVector(final int column) + throws MatrixIndexException { + + checkColumnIndex(column); + final T[] outData = buildArray(getField(), rows); + + // perform copy block-wise, to ensure good cache behavior + final int jBlock = column / BLOCK_SIZE; + final int jColumn = column - jBlock * BLOCK_SIZE; + final int jWidth = blockWidth(jBlock); + int outIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = 0; i < iHeight; ++i) { + outData[outIndex++] = block[i * jWidth + jColumn]; + } + } + + return new FieldVectorImpl(outData, false); + + } + + /** {@inheritDoc} */ + @Override + public void setColumnVector(final int column, final FieldVector vector) + throws MatrixIndexException, InvalidMatrixException { + try { + setColumn(column, ((FieldVectorImpl) vector).getDataRef()); + } catch (ClassCastException cce) { + super.setColumnVector(column, vector); + } + } + + /** {@inheritDoc} */ + @Override + public T[] getRow(final int row) + throws MatrixIndexException { + + checkRowIndex(row); + final T[] out = buildArray(getField(), columns); + + // perform copy block-wise, to ensure good cache behavior + final int iBlock = row / BLOCK_SIZE; + final int iRow = row - iBlock * BLOCK_SIZE; + int outIndex = 0; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + System.arraycopy(block, iRow * jWidth, out, outIndex, jWidth); + outIndex += jWidth; + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public void setRow(final int row, final T[] array) + throws MatrixIndexException, InvalidMatrixException { + + checkRowIndex(row); + final int nCols = getColumnDimension(); + if (array.length != nCols) { + throw new InvalidMatrixException( + "dimensions mismatch: got {0}x{1} but expected {2}x{3}", + 1, array.length, 1, nCols); + } + + // perform copy block-wise, to ensure good cache behavior + final int iBlock = row / BLOCK_SIZE; + final int iRow = row - iBlock * BLOCK_SIZE; + int outIndex = 0; + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + System.arraycopy(array, outIndex, block, iRow * jWidth, jWidth); + outIndex += jWidth; + } + + } + + /** {@inheritDoc} */ + @Override + public T[] getColumn(final int column) + throws MatrixIndexException { + + checkColumnIndex(column); + final T[] out = buildArray(getField(), rows); + + // perform copy block-wise, to ensure good cache behavior + final int jBlock = column / BLOCK_SIZE; + final int jColumn = column - jBlock * BLOCK_SIZE; + final int jWidth = blockWidth(jBlock); + int outIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = 0; i < iHeight; ++i) { + out[outIndex++] = block[i * jWidth + jColumn]; + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public void setColumn(final int column, final T[] array) + throws MatrixIndexException, InvalidMatrixException { + + checkColumnIndex(column); + final int nRows = getRowDimension(); + if (array.length != nRows) { + throw new InvalidMatrixException( + "dimensions mismatch: got {0}x{1} but expected {2}x{3}", + array.length, 1, nRows, 1); + } + + // perform copy block-wise, to ensure good cache behavior + final int jBlock = column / BLOCK_SIZE; + final int jColumn = column - jBlock * BLOCK_SIZE; + final int jWidth = blockWidth(jBlock); + int outIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int iHeight = blockHeight(iBlock); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int i = 0; i < iHeight; ++i) { + block[i * jWidth + jColumn] = array[outIndex++]; + } + } + + } + + /** {@inheritDoc} */ + @Override + public T getEntry(final int row, final int column) + throws MatrixIndexException { + try { + final int iBlock = row / BLOCK_SIZE; + final int jBlock = column / BLOCK_SIZE; + final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) + + (column - jBlock * BLOCK_SIZE); + return blocks[iBlock * blockColumns + jBlock][k]; + } catch (ArrayIndexOutOfBoundsException e) { + throw new MatrixIndexException( + "no entry at indices ({0}, {1}) in a {2}x{3} matrix", + row, column, getRowDimension(), getColumnDimension()); + } + } + + /** {@inheritDoc} */ + @Override + public void setEntry(final int row, final int column, final T value) + throws MatrixIndexException { + try { + final int iBlock = row / BLOCK_SIZE; + final int jBlock = column / BLOCK_SIZE; + final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) + + (column - jBlock * BLOCK_SIZE); + blocks[iBlock * blockColumns + jBlock][k] = value; + } catch (ArrayIndexOutOfBoundsException e) { + throw new MatrixIndexException( + "no entry at indices ({0}, {1}) in a {2}x{3} matrix", + row, column, getRowDimension(), getColumnDimension()); + } + } + + /** {@inheritDoc} */ + @Override + public void addToEntry(final int row, final int column, final T increment) + throws MatrixIndexException { + try { + final int iBlock = row / BLOCK_SIZE; + final int jBlock = column / BLOCK_SIZE; + final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) + + (column - jBlock * BLOCK_SIZE); + final T[] blockIJ = blocks[iBlock * blockColumns + jBlock]; + blockIJ[k] = blockIJ[k].add(increment); + } catch (ArrayIndexOutOfBoundsException e) { + throw new MatrixIndexException( + "no entry at indices ({0}, {1}) in a {2}x{3} matrix", + row, column, getRowDimension(), getColumnDimension()); + } + } + + /** {@inheritDoc} */ + @Override + public void multiplyEntry(final int row, final int column, final T factor) + throws MatrixIndexException { + try { + final int iBlock = row / BLOCK_SIZE; + final int jBlock = column / BLOCK_SIZE; + final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) + + (column - jBlock * BLOCK_SIZE); + final T[] blockIJ = blocks[iBlock * blockColumns + jBlock]; + blockIJ[k] = blockIJ[k].multiply(factor); + } catch (ArrayIndexOutOfBoundsException e) { + throw new MatrixIndexException( + "no entry at indices ({0}, {1}) in a {2}x{3} matrix", + row, column, getRowDimension(), getColumnDimension()); + } + } + + /** {@inheritDoc} */ + @Override + public FieldMatrix transpose() { + + final int nRows = getRowDimension(); + final int nCols = getColumnDimension(); + final DenseFieldMatrix out = new DenseFieldMatrix(getField(), nCols, nRows); + + // perform transpose block-wise, to ensure good cache behavior + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockColumns; ++iBlock) { + for (int jBlock = 0; jBlock < blockRows; ++jBlock) { + + // transpose current block + final T[] outBlock = out.blocks[blockIndex]; + final T[] tBlock = blocks[jBlock * blockColumns + iBlock]; + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, columns); + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, rows); + for (int p = pStart, k = 0; p < pEnd; ++p) { + final int lInc = pEnd - pStart; + for (int q = qStart, l = p - pStart; q < qEnd; ++q, l+= lInc) { + outBlock[k++] = tBlock[l]; + } + } + + // go to next block + ++blockIndex; + + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public int getRowDimension() { + return rows; + } + + /** {@inheritDoc} */ + @Override + public int getColumnDimension() { + return columns; + } + + /** {@inheritDoc} */ + @Override + public T[] operate(final T[] v) + throws IllegalArgumentException { + + if (v.length != columns) { + throw MathRuntimeException.createIllegalArgumentException( + "vector length mismatch: got {0} but expected {1}", + v.length, columns); + } + final T[] out = buildArray(getField(), rows); + final T zero = getField().getZero(); + + // perform multiplication block-wise, to ensure good cache behavior + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final T[] block = blocks[iBlock * blockColumns + jBlock]; + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + for (int p = pStart, k = 0; p < pEnd; ++p) { + T sum = zero; + int q = qStart; + while (q < qEnd - 3) { + sum = sum. + add(block[k].multiply(v[q])). + add(block[k + 1].multiply(v[q + 1])). + add(block[k + 2].multiply(v[q + 2])). + add(block[k + 3].multiply(v[q + 3])); + k += 4; + q += 4; + } + while (q < qEnd) { + sum = sum.add(block[k++].multiply(v[q++])); + } + out[p] = out[p].add(sum); + } + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public T[] preMultiply(final T[] v) + throws IllegalArgumentException { + + if (v.length != rows) { + throw MathRuntimeException.createIllegalArgumentException( + "vector length mismatch: got {0} but expected {1}", + v.length, rows); + } + final T[] out = buildArray(getField(), columns); + final T zero = getField().getZero(); + + // perform multiplication block-wise, to ensure good cache behavior + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int jWidth2 = jWidth + jWidth; + final int jWidth3 = jWidth2 + jWidth; + final int jWidth4 = jWidth3 + jWidth; + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final T[] block = blocks[iBlock * blockColumns + jBlock]; + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int q = qStart; q < qEnd; ++q) { + int k = q - qStart; + T sum = zero; + int p = pStart; + while (p < pEnd - 3) { + sum = sum. + add(block[k].multiply(v[p])). + add(block[k + jWidth].multiply(v[p + 1])). + add(block[k + jWidth2].multiply(v[p + 2])). + add(block[k + jWidth3].multiply(v[p + 3])); + k += jWidth4; + p += 4; + } + while (p < pEnd) { + sum = sum.add(block[k].multiply(v[p++])); + k += jWidth; + } + out[q] = out[q].add(sum); + } + } + } + + return out; + + } + + /** {@inheritDoc} */ + @Override + public T walkInRowOrder(final FieldMatrixChangingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int p = pStart; p < pEnd; ++p) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + block[k] = visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInRowOrder(final FieldMatrixPreservingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int p = pStart; p < pEnd; ++p) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInRowOrder(final FieldMatrixChangingVisitor 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 iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) { + final int p0 = iBlock * BLOCK_SIZE; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow); + for (int p = pStart; p < pEnd; ++p) { + for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int q0 = jBlock * BLOCK_SIZE; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + block[k] = visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInRowOrder(final FieldMatrixPreservingVisitor 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 iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) { + final int p0 = iBlock * BLOCK_SIZE; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow); + for (int p = pStart; p < pEnd; ++p) { + for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int q0 = jBlock * BLOCK_SIZE; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInOptimizedOrder(final FieldMatrixChangingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final T[] block = blocks[blockIndex]; + for (int p = pStart, k = 0; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q, ++k) { + block[k] = visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor visitor) + throws MatrixVisitorException { + visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); + for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + final int pStart = iBlock * BLOCK_SIZE; + final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); + for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + final int qStart = jBlock * BLOCK_SIZE; + final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); + final T[] block = blocks[blockIndex]; + for (int p = pStart, k = 0; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q, ++k) { + visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInOptimizedOrder(final FieldMatrixChangingVisitor 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 iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) { + final int p0 = iBlock * BLOCK_SIZE; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow); + for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int q0 = jBlock * BLOCK_SIZE; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + block[k] = visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** {@inheritDoc} */ + @Override + public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor 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 iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) { + final int p0 = iBlock * BLOCK_SIZE; + final int pStart = Math.max(startRow, p0); + final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow); + for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) { + final int jWidth = blockWidth(jBlock); + final int q0 = jBlock * BLOCK_SIZE; + final int qStart = Math.max(startColumn, q0); + final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); + final T[] block = blocks[iBlock * blockColumns + jBlock]; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + visitor.visit(p, q, block[k]); + } + } + } + } + return visitor.end(); + } + + /** + * Get the height of a block. + * @param blockRow row index (in block sense) of the block + * @return height (number of rows) of the block + */ + private int blockHeight(final int blockRow) { + return (blockRow == blockRows - 1) ? rows - blockRow * BLOCK_SIZE : BLOCK_SIZE; + } + + /** + * Get the width of a block. + * @param blockColumn column index (in block sense) of the block + * @return width (number of columns) of the block + */ + private int blockWidth(final int blockColumn) { + return (blockColumn == blockColumns - 1) ? columns - blockColumn * BLOCK_SIZE : BLOCK_SIZE; + } + +} diff --git a/src/java/org/apache/commons/math/linear/MatrixUtils.java b/src/java/org/apache/commons/math/linear/MatrixUtils.java index 0e468650e..6299ac0e1 100644 --- a/src/java/org/apache/commons/math/linear/MatrixUtils.java +++ b/src/java/org/apache/commons/math/linear/MatrixUtils.java @@ -17,9 +17,13 @@ package org.apache.commons.math.linear; +import java.lang.reflect.Array; import java.math.BigDecimal; import java.util.Arrays; +import org.apache.commons.math.Field; +import org.apache.commons.math.FieldElement; + /** * A collection of static methods that operate on or return matrices. * @@ -46,6 +50,22 @@ public class MatrixUtils { return new DenseRealMatrix(rows, columns); } + /** + * Returns a {@link FieldMatrix} with specified dimensions. + *

The matrix elements are all set to field.getZero().

+ * @param field field to which the matrix elements belong + * @param rows number of rows of the matrix + * @param columns number of columns of the matrix + * @return FieldMatrix with specified dimensions + * @see #createFieldMatrix(FieldElement[][]) + * @since 2.0 + */ + public static > FieldMatrix createFieldMatrix(final Field field, + final int rows, + final int columns) { + return new DenseFieldMatrix(field, rows, columns); + } + /** * Returns a {@link RealMatrix} whose entries are the the values in the * the input array. The input array is copied, not referenced. @@ -61,6 +81,24 @@ public class MatrixUtils { return new DenseRealMatrix(data); } + /** + * Returns a {@link FieldMatrix} whose entries are the the values in the + * the input array. + *

+ * The input array is copied, not referenced. + *

+ * @param data input array + * @return RealMatrix containing the values of the array + * @throws IllegalArgumentException if data is not rectangular + * (not all rows have the same length) or empty + * @throws NullPointerException if data is null + * @see #createFieldMatrix(Field, int, int) + * @since 2.0 + */ + public static > FieldMatrix createFieldMatrix(T[][] data) { + return new DenseFieldMatrix(data); + } + /** * Returns dimension x dimension identity matrix. * @@ -76,6 +114,48 @@ public class MatrixUtils { } return m; } + + /** + * Returns dimension x dimension identity matrix. + * + * @param dimension dimension of identity matrix to generate + * @return identity matrix + * @throws IllegalArgumentException if dimension is not positive + * @since 2.0 + */ + @SuppressWarnings("unchecked") + public static > FieldMatrix + createFieldIdentityMatrix(final Field field, final int dimension) { + final T zero = field.getZero(); + final T one = field.getOne(); + final T[][] d = (T[][]) Array.newInstance(zero.getClass(), dimension, dimension); + for (int row = 0; row < dimension; row++) { + final T[] dRow = d[row]; + Arrays.fill(dRow, zero); + dRow[row] = one; + } + return new FieldMatrixImpl(d, false); + } + + /** + * Returns dimension x dimension identity matrix. + * + * @param dimension dimension of identity matrix to generate + * @return identity matrix + * @throws IllegalArgumentException if dimension is not positive + * @since 1.1 + * @deprecated since 2.0, replaced by {@link #createFieldIdentityMatrix(Field, int)} + */ + @Deprecated + public static BigMatrix createBigIdentityMatrix(int dimension) { + final BigDecimal[][] d = new BigDecimal[dimension][dimension]; + for (int row = 0; row < dimension; row++) { + final BigDecimal[] dRow = d[row]; + Arrays.fill(dRow, BigMatrixImpl.ZERO); + dRow[row] = BigMatrixImpl.ONE; + } + return new BigMatrixImpl(d, false); + } /** * Returns a diagonal matrix with specified elements. @@ -93,6 +173,24 @@ public class MatrixUtils { return m; } + /** + * Returns a diagonal matrix with specified elements. + * + * @param diagonal diagonal elements of the matrix (the array elements + * will be copied) + * @return diagonal matrix + * @since 2.0 + */ + public static > FieldMatrix + createFieldDiagonalMatrix(final T[] diagonal) { + final FieldMatrix m = + createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length); + for (int i = 0; i < diagonal.length; ++i) { + m.setEntry(i, i, diagonal[i]); + } + return m; + } + /** * Returns a {@link BigMatrix} whose entries are the the values in the * the input array. The input array is copied, not referenced. @@ -102,7 +200,9 @@ public class MatrixUtils { * @throws IllegalArgumentException if data is not rectangular * (not all rows have the same length) or empty * @throws NullPointerException if data is null + * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])} */ + @Deprecated public static BigMatrix createBigMatrix(double[][] data) { return new BigMatrixImpl(data); } @@ -116,7 +216,9 @@ public class MatrixUtils { * @throws IllegalArgumentException if data is not rectangular * (not all rows have the same length) or empty * @throws NullPointerException if data is null + * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])} */ + @Deprecated public static BigMatrix createBigMatrix(BigDecimal[][] data) { return new BigMatrixImpl(data); } @@ -136,7 +238,9 @@ public class MatrixUtils { * (not all rows have the same length) or empty * @throws NullPointerException if data is null * @see #createRealMatrix(double[][]) + * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])} */ + @Deprecated public static BigMatrix createBigMatrix(BigDecimal[][] data, boolean copyArray) { return new BigMatrixImpl(data, copyArray); } @@ -150,7 +254,9 @@ public class MatrixUtils { * @throws IllegalArgumentException if data is not rectangular * (not all rows have the same length) or empty * @throws NullPointerException if data is null + * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])} */ + @Deprecated public static BigMatrix createBigMatrix(String[][] data) { return new BigMatrixImpl(data); } @@ -167,6 +273,18 @@ public class MatrixUtils { return new RealVectorImpl(data, true); } + /** + * Creates a {@link FieldVector} using the data from the input array. + * + * @param data the input data + * @return a data.length FieldVector + * @throws IllegalArgumentException if data is empty + * @throws NullPointerException if datais null + */ + public static > FieldVector createFieldVector(final T[] data) { + return new FieldVectorImpl(data, true); + } + /** * Creates a row {@link RealMatrix} using the data from the input * array. @@ -185,6 +303,25 @@ public class MatrixUtils { return m; } + /** + * Creates a row {@link FieldMatrix} using the data from the input + * array. + * + * @param rowData the input row data + * @return a 1 x rowData.length FieldMatrix + * @throws IllegalArgumentException if rowData is empty + * @throws NullPointerException if rowDatais null + */ + public static > FieldMatrix + createRowFieldMatrix(final T[] rowData) { + final int nCols = rowData.length; + final FieldMatrix m = createFieldMatrix(rowData[0].getField(), 1, nCols); + for (int i = 0; i < nCols; ++i) { + m.setEntry(0, i, rowData[i]); + } + return m; + } + /** * Creates a row {@link BigMatrix} using the data from the input * array. @@ -193,7 +330,9 @@ public class MatrixUtils { * @return a 1 x rowData.length BigMatrix * @throws IllegalArgumentException if rowData is empty * @throws NullPointerException if rowDatais null + * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createRowBigMatrix(double[] rowData) { final int nCols = rowData.length; final BigDecimal[][] data = new BigDecimal[1][nCols]; @@ -211,7 +350,9 @@ public class MatrixUtils { * @return a 1 x rowData.length BigMatrix * @throws IllegalArgumentException if rowData is empty * @throws NullPointerException if rowDatais null + * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createRowBigMatrix(BigDecimal[] rowData) { final int nCols = rowData.length; final BigDecimal[][] data = new BigDecimal[1][nCols]; @@ -227,7 +368,9 @@ public class MatrixUtils { * @return a 1 x rowData.length BigMatrix * @throws IllegalArgumentException if rowData is empty * @throws NullPointerException if rowDatais null + * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createRowBigMatrix(String[] rowData) { final int nCols = rowData.length; final BigDecimal[][] data = new BigDecimal[1][nCols]; @@ -255,6 +398,25 @@ public class MatrixUtils { return m; } + /** + * Creates a column {@link FieldMatrix} using the data from the input + * array. + * + * @param columnData the input column data + * @return a columnData x 1 FieldMatrix + * @throws IllegalArgumentException if columnData is empty + * @throws NullPointerException if columnDatais null + */ + public static > FieldMatrix + createColumnFieldMatrix(final T[] columnData) { + final int nRows = columnData.length; + final FieldMatrix m = createFieldMatrix(columnData[0].getField(), nRows, 1); + for (int i = 0; i < nRows; ++i) { + m.setEntry(i, 0, columnData[i]); + } + return m; + } + /** * Creates a column {@link BigMatrix} using the data from the input * array. @@ -263,7 +425,9 @@ public class MatrixUtils { * @return a columnData x 1 BigMatrix * @throws IllegalArgumentException if columnData is empty * @throws NullPointerException if columnDatais null + * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createColumnBigMatrix(double[] columnData) { final int nRows = columnData.length; final BigDecimal[][] data = new BigDecimal[nRows][1]; @@ -281,7 +445,9 @@ public class MatrixUtils { * @return a columnData x 1 BigMatrix * @throws IllegalArgumentException if columnData is empty * @throws NullPointerException if columnDatais null + * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createColumnBigMatrix(BigDecimal[] columnData) { final int nRows = columnData.length; final BigDecimal[][] data = new BigDecimal[nRows][1]; @@ -299,7 +465,9 @@ public class MatrixUtils { * @return a columnData x 1 BigMatrix * @throws IllegalArgumentException if columnData is empty * @throws NullPointerException if columnDatais null + * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])} */ + @Deprecated public static BigMatrix createColumnBigMatrix(String[] columnData) { int nRows = columnData.length; final BigDecimal[][] data = new BigDecimal[nRows][1]; @@ -309,23 +477,5 @@ public class MatrixUtils { return new BigMatrixImpl(data, false); } - /** - * Returns dimension x dimension identity matrix. - * - * @param dimension dimension of identity matrix to generate - * @return identity matrix - * @throws IllegalArgumentException if dimension is not positive - * @since 1.1 - */ - public static BigMatrix createBigIdentityMatrix(int dimension) { - final BigDecimal[][] d = new BigDecimal[dimension][dimension]; - for (int row = 0; row < dimension; row++) { - final BigDecimal[] dRow = d[row]; - Arrays.fill(dRow, BigMatrixImpl.ZERO); - dRow[row] = BigMatrixImpl.ONE; - } - return new BigMatrixImpl(d, false); - } - } diff --git a/src/test/org/apache/commons/math/linear/DenseFieldMatrixTest.java b/src/test/org/apache/commons/math/linear/DenseFieldMatrixTest.java new file mode 100644 index 000000000..6ffe68066 --- /dev/null +++ b/src/test/org/apache/commons/math/linear/DenseFieldMatrixTest.java @@ -0,0 +1,1294 @@ +/* + * 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; + +import org.apache.commons.math.TestUtils; +import org.apache.commons.math.fraction.Fraction; +import org.apache.commons.math.fraction.FractionField; +import org.apache.commons.math.linear.decomposition.FieldLUDecompositionImpl; +import org.apache.commons.math.linear.decomposition.NonSquareMatrixException; + +/** + * Test cases for the {@link DenseFieldMatrix} class. + * + * @version $Revision$ $Date$ + */ + +public final class DenseFieldMatrixTest extends TestCase { + + // 3 x 3 identity matrix + protected Fraction[][] id = { + {new Fraction(1),new Fraction(0),new Fraction(0)}, + {new Fraction(0),new Fraction(1),new Fraction(0)}, + {new Fraction(0),new Fraction(0),new Fraction(1)} + }; + + // Test data for group operations + protected Fraction[][] testData = { + {new Fraction(1),new Fraction(2),new Fraction(3)}, + {new Fraction(2),new Fraction(5),new Fraction(3)}, + {new Fraction(1),new Fraction(0),new Fraction(8)} + }; + protected Fraction[][] testDataLU = { + {new Fraction(2), new Fraction(5), new Fraction(3)}, + {new Fraction(1, 2), new Fraction(-5, 2), new Fraction(13, 2)}, + {new Fraction(1, 2), new Fraction(1, 5), new Fraction(1, 5)} + }; + protected Fraction[][] testDataPlus2 = { + {new Fraction(3),new Fraction(4),new Fraction(5)}, + {new Fraction(4),new Fraction(7),new Fraction(5)}, + {new Fraction(3),new Fraction(2),new Fraction(10)} + }; + protected Fraction[][] testDataMinus = { + {new Fraction(-1),new Fraction(-2),new Fraction(-3)}, + {new Fraction(-2),new Fraction(-5),new Fraction(-3)}, + {new Fraction(-1),new Fraction(0),new Fraction(-8)} + }; + protected Fraction[] testDataRow1 = {new Fraction(1),new Fraction(2),new Fraction(3)}; + protected Fraction[] testDataCol3 = {new Fraction(3),new Fraction(3),new Fraction(8)}; + protected Fraction[][] testDataInv = { + {new Fraction(-40),new Fraction(16),new Fraction(9)}, + {new Fraction(13),new Fraction(-5),new Fraction(-3)}, + {new Fraction(5),new Fraction(-2),new Fraction(-1)} + }; + protected Fraction[] preMultTest = {new Fraction(8), new Fraction(12), new Fraction(33)}; + protected Fraction[][] testData2 = { + {new Fraction(1),new Fraction(2),new Fraction(3)}, + {new Fraction(2),new Fraction(5),new Fraction(3)} + }; + protected Fraction[][] testData2T = { + {new Fraction(1),new Fraction(2)}, + {new Fraction(2),new Fraction(5)}, + {new Fraction(3),new Fraction(3)} + }; + protected Fraction[][] testDataPlusInv = { + {new Fraction(-39),new Fraction(18),new Fraction(12)}, + {new Fraction(15),new Fraction(0),new Fraction(0)}, + {new Fraction(6),new Fraction(-2),new Fraction(7)} + }; + + // lu decomposition tests + protected Fraction[][] luData = { + {new Fraction(2),new Fraction(3),new Fraction(3)}, + {new Fraction(0),new Fraction(5),new Fraction(7)}, + {new Fraction(6),new Fraction(9),new Fraction(8)} + }; + protected Fraction[][] luDataLUDecomposition = { + {new Fraction(6),new Fraction(9),new Fraction(8)}, + {new Fraction(0),new Fraction(5),new Fraction(7)}, + {new Fraction(1, 3),new Fraction(0),new Fraction(1, 3)} + }; + + // singular matrices + protected Fraction[][] singular = { {new Fraction(2),new Fraction(3)}, {new Fraction(2),new Fraction(3)} }; + protected Fraction[][] bigSingular = { + {new Fraction(1),new Fraction(2),new Fraction(3),new Fraction(4)}, + {new Fraction(2),new Fraction(5),new Fraction(3),new Fraction(4)}, + {new Fraction(7),new Fraction(3),new Fraction(256),new Fraction(1930)}, + {new Fraction(3),new Fraction(7),new Fraction(6),new Fraction(8)} + }; // 4th row = 1st + 2nd + protected Fraction[][] detData = { + {new Fraction(1),new Fraction(2),new Fraction(3)}, + {new Fraction(4),new Fraction(5),new Fraction(6)}, + {new Fraction(7),new Fraction(8),new Fraction(10)} + }; + protected Fraction[][] detData2 = { {new Fraction(1), new Fraction(3)}, {new Fraction(2), new Fraction(4)}}; + + // vectors + protected Fraction[] testVector = {new Fraction(1),new Fraction(2),new Fraction(3)}; + protected Fraction[] testVector2 = {new Fraction(1),new Fraction(2),new Fraction(3),new Fraction(4)}; + + // submatrix accessor tests + protected Fraction[][] subTestData = { + {new Fraction(1), new Fraction(2), new Fraction(3), new Fraction(4)}, + {new Fraction(3, 2), new Fraction(5, 2), new Fraction(7, 2), new Fraction(9, 2)}, + {new Fraction(2), new Fraction(4), new Fraction(6), new Fraction(8)}, + {new Fraction(4), new Fraction(5), new Fraction(6), new Fraction(7)} + }; + // array selections + protected Fraction[][] subRows02Cols13 = { {new Fraction(2), new Fraction(4)}, {new Fraction(4), new Fraction(8)}}; + protected Fraction[][] subRows03Cols12 = { {new Fraction(2), new Fraction(3)}, {new Fraction(5), new Fraction(6)}}; + protected Fraction[][] subRows03Cols123 = { + {new Fraction(2), new Fraction(3), new Fraction(4)}, + {new Fraction(5), new Fraction(6), new Fraction(7)} + }; + // effective permutations + protected Fraction[][] subRows20Cols123 = { + {new Fraction(4), new Fraction(6), new Fraction(8)}, + {new Fraction(2), new Fraction(3), new Fraction(4)} + }; + protected Fraction[][] subRows31Cols31 = {{new Fraction(7), new Fraction(5)}, {new Fraction(9, 2), new Fraction(5, 2)}}; + // contiguous ranges + protected Fraction[][] subRows01Cols23 = {{new Fraction(3),new Fraction(4)} , {new Fraction(7, 2), new Fraction(9, 2)}}; + protected Fraction[][] subRows23Cols00 = {{new Fraction(2)} , {new Fraction(4)}}; + protected Fraction[][] subRows00Cols33 = {{new Fraction(4)}}; + // row matrices + protected Fraction[][] subRow0 = {{new Fraction(1),new Fraction(2),new Fraction(3),new Fraction(4)}}; + protected Fraction[][] subRow3 = {{new Fraction(4),new Fraction(5),new Fraction(6),new Fraction(7)}}; + // column matrices + protected Fraction[][] subColumn1 = {{new Fraction(2)}, {new Fraction(5, 2)}, {new Fraction(4)}, {new Fraction(5)}}; + protected Fraction[][] subColumn3 = {{new Fraction(4)}, {new Fraction(9, 2)}, {new Fraction(8)}, {new Fraction(7)}}; + + // tolerances + protected double entryTolerance = 10E-16; + protected double normTolerance = 10E-14; + + public DenseFieldMatrixTest(String name) { + super(name); + } + + public static Test suite() { + TestSuite suite = new TestSuite(DenseFieldMatrixTest.class); + suite.setName("DenseFieldMatrix Tests"); + return suite; + } + + /** test dimensions */ + public void testDimensions() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix m2 = new DenseFieldMatrix(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); + DenseFieldMatrix m1 = createRandomMatrix(r, 47, 83); + DenseFieldMatrix m2 = new DenseFieldMatrix(m1.getData()); + assertEquals(m1, m2); + DenseFieldMatrix m3 = new DenseFieldMatrix(testData); + DenseFieldMatrix m4 = new DenseFieldMatrix(m3.getData()); + assertEquals(m3, m4); + } + + /** test add */ + public void testAdd() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix mInv = new DenseFieldMatrix(testDataInv); + FieldMatrix mPlusMInv = m.add(mInv); + Fraction[][] sumEntries = mPlusMInv.getData(); + for (int row = 0; row < m.getRowDimension(); row++) { + for (int col = 0; col < m.getColumnDimension(); col++) { + assertEquals(testDataPlusInv[row][col],sumEntries[row][col]); + } + } + } + + /** test add failure */ + public void testAddFail() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix m2 = new DenseFieldMatrix(testData2); + try { + m.add(m2); + fail("IllegalArgumentException expected"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + + /** test m-n = m + -n */ + public void testPlusMinus() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix m2 = new DenseFieldMatrix(testDataInv); + TestUtils.assertEquals(m.subtract(m2), m2.scalarMultiply(new Fraction(-1)).add(m)); + try { + m.subtract(new DenseFieldMatrix(testData2)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + + /** test multiply */ + public void testMultiply() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix mInv = new DenseFieldMatrix(testDataInv); + DenseFieldMatrix identity = new DenseFieldMatrix(id); + DenseFieldMatrix m2 = new DenseFieldMatrix(testData2); + TestUtils.assertEquals(m.multiply(mInv), identity); + TestUtils.assertEquals(mInv.multiply(m), identity); + TestUtils.assertEquals(m.multiply(identity), m); + TestUtils.assertEquals(identity.multiply(mInv), mInv); + TestUtils.assertEquals(m2.multiply(identity), m2); + try { + m.multiply(new DenseFieldMatrix(bigSingular)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + // expected + } + } + + public void testSeveralBlocks() { + + FieldMatrix m = + new DenseFieldMatrix(FractionField.getInstance(), 37, 41); + for (int i = 0; i < m.getRowDimension(); ++i) { + for (int j = 0; j < m.getColumnDimension(); ++j) { + m.setEntry(i, j, new Fraction(i * 11 + j, 11)); + } + } + + FieldMatrix 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)); + } + } + + FieldMatrix mPm = m.add(m); + for (int i = 0; i < mPm.getRowDimension(); ++i) { + for (int j = 0; j < mPm.getColumnDimension(); ++j) { + assertEquals(m.getEntry(i, j).multiply(new Fraction(2)), mPm.getEntry(i, j)); + } + } + + FieldMatrix 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)); + } + } + + FieldMatrix mTm = mT.multiply(m); + for (int i = 0; i < mTm.getRowDimension(); ++i) { + for (int j = 0; j < mTm.getColumnDimension(); ++j) { + Fraction sum = Fraction.ZERO; + for (int k = 0; k < mT.getColumnDimension(); ++k) { + sum = sum.add(new Fraction(k * 11 + i, 11).multiply(new Fraction(k * 11 + j, 11))); + } + assertEquals(sum, mTm.getEntry(i, j)); + } + } + + FieldMatrix mmT = m.multiply(mT); + for (int i = 0; i < mmT.getRowDimension(); ++i) { + for (int j = 0; j < mmT.getColumnDimension(); ++j) { + Fraction sum = Fraction.ZERO; + for (int k = 0; k < m.getColumnDimension(); ++k) { + sum = sum.add(new Fraction(i * 11 + k, 11).multiply(new Fraction(j * 11 + k, 11))); + } + assertEquals(sum, mmT.getEntry(i, j)); + } + } + + FieldMatrix sub1 = m.getSubMatrix(2, 9, 5, 20); + for (int i = 0; i < sub1.getRowDimension(); ++i) { + for (int j = 0; j < sub1.getColumnDimension(); ++j) { + assertEquals(new Fraction((i + 2) * 11 + (j + 5), 11), sub1.getEntry(i, j)); + } + } + + FieldMatrix sub2 = m.getSubMatrix(10, 12, 3, 40); + for (int i = 0; i < sub2.getRowDimension(); ++i) { + for (int j = 0; j < sub2.getColumnDimension(); ++j) { + assertEquals(new Fraction((i + 10) * 11 + (j + 3), 11), sub2.getEntry(i, j)); + } + } + + FieldMatrix sub3 = m.getSubMatrix(30, 34, 0, 5); + for (int i = 0; i < sub3.getRowDimension(); ++i) { + for (int j = 0; j < sub3.getColumnDimension(); ++j) { + assertEquals(new Fraction((i + 30) * 11 + (j + 0), 11), sub3.getEntry(i, j)); + } + } + + FieldMatrix sub4 = m.getSubMatrix(30, 32, 32, 35); + for (int i = 0; i < sub4.getRowDimension(); ++i) { + for (int j = 0; j < sub4.getColumnDimension(); ++j) { + assertEquals(new Fraction((i + 30) * 11 + (j + 32), 11), sub4.getEntry(i, j)); + } + } + + } + + //Additional Test for DenseFieldMatrixTest.testMultiply + + private Fraction[][] d3 = new Fraction[][] { + {new Fraction(1),new Fraction(2),new Fraction(3),new Fraction(4)}, + {new Fraction(5),new Fraction(6),new Fraction(7),new Fraction(8)} + }; + private Fraction[][] d4 = new Fraction[][] { + {new Fraction(1)}, + {new Fraction(2)}, + {new Fraction(3)}, + {new Fraction(4)} + }; + private Fraction[][] d5 = new Fraction[][] {{new Fraction(30)},{new Fraction(70)}}; + + public void testMultiply2() { + FieldMatrix m3 = new DenseFieldMatrix(d3); + FieldMatrix m4 = new DenseFieldMatrix(d4); + FieldMatrix m5 = new DenseFieldMatrix(d5); + TestUtils.assertEquals(m3.multiply(m4), m5); + } + + /** test trace */ + public void testTrace() { + FieldMatrix m = new DenseFieldMatrix(id); + assertEquals(new Fraction(3),m.getTrace()); + m = new DenseFieldMatrix(testData2); + try { + m.getTrace(); + fail("Expecting NonSquareMatrixException"); + } catch (NonSquareMatrixException ex) { + // ignored + } + } + + /** test scalarAdd */ + public void testScalarAdd() { + FieldMatrix m = new DenseFieldMatrix(testData); + TestUtils.assertEquals(new DenseFieldMatrix(testDataPlus2), + m.scalarAdd(new Fraction(2))); + } + + /** test operate */ + public void testOperate() { + FieldMatrix m = new DenseFieldMatrix(id); + TestUtils.assertEquals(testVector, m.operate(testVector)); + TestUtils.assertEquals(testVector, m.operate(new FieldVectorImpl(testVector)).getData()); + m = new DenseFieldMatrix(bigSingular); + try { + m.operate(testVector); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + + public void testOperateLarge() { + int p = (3 * DenseFieldMatrix.BLOCK_SIZE) / 2; + int q = (5 * DenseFieldMatrix.BLOCK_SIZE) / 2; + int r = 2 * DenseFieldMatrix.BLOCK_SIZE; + Random random = new Random(111007463902334l); + FieldMatrix m1 = createRandomMatrix(random, p, q); + FieldMatrix m2 = createRandomMatrix(random, q, r); + FieldMatrix m1m2 = m1.multiply(m2); + for (int i = 0; i < r; ++i) { + TestUtils.assertEquals(m1m2.getColumn(i), m1.operate(m2.getColumn(i))); + } + } + + public void testOperatePremultiplyLarge() { + int p = (3 * DenseFieldMatrix.BLOCK_SIZE) / 2; + int q = (5 * DenseFieldMatrix.BLOCK_SIZE) / 2; + int r = 2 * DenseFieldMatrix.BLOCK_SIZE; + Random random = new Random(111007463902334l); + FieldMatrix m1 = createRandomMatrix(random, p, q); + FieldMatrix m2 = createRandomMatrix(random, q, r); + FieldMatrix m1m2 = m1.multiply(m2); + for (int i = 0; i < p; ++i) { + TestUtils.assertEquals(m1m2.getRow(i), m2.preMultiply(m1.getRow(i))); + } + } + + /** test issue MATH-209 */ + public void testMath209() { + FieldMatrix a = new DenseFieldMatrix(new Fraction[][] { + { new Fraction(1), new Fraction(2) }, + { new Fraction(3), new Fraction(4) }, + { new Fraction(5), new Fraction(6) } + }); + Fraction[] b = a.operate(new Fraction[] { new Fraction(1), new Fraction(1) }); + assertEquals(a.getRowDimension(), b.length); + assertEquals( new Fraction(3), b[0]); + assertEquals( new Fraction(7), b[1]); + assertEquals(new Fraction(11), b[2]); + } + + /** test transpose */ + public void testTranspose() { + FieldMatrix m = new DenseFieldMatrix(testData); + FieldMatrix mIT = new FieldLUDecompositionImpl(m).getSolver().getInverse().transpose(); + FieldMatrix mTI = new FieldLUDecompositionImpl(m.transpose()).getSolver().getInverse(); + TestUtils.assertEquals(mIT, mTI); + m = new DenseFieldMatrix(testData2); + FieldMatrix mt = new DenseFieldMatrix(testData2T); + TestUtils.assertEquals(mt, m.transpose()); + } + + /** test preMultiply by vector */ + public void testPremultiplyVector() { + FieldMatrix m = new DenseFieldMatrix(testData); + TestUtils.assertEquals(m.preMultiply(testVector), preMultTest); + TestUtils.assertEquals(m.preMultiply(new FieldVectorImpl(testVector).getData()), + preMultTest); + m = new DenseFieldMatrix(bigSingular); + try { + m.preMultiply(testVector); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + + public void testPremultiply() { + FieldMatrix m3 = new DenseFieldMatrix(d3); + FieldMatrix m4 = new DenseFieldMatrix(d4); + FieldMatrix m5 = new DenseFieldMatrix(d5); + TestUtils.assertEquals(m4.preMultiply(m3), m5); + + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix mInv = new DenseFieldMatrix(testDataInv); + DenseFieldMatrix identity = new DenseFieldMatrix(id); + TestUtils.assertEquals(m.preMultiply(mInv), identity); + TestUtils.assertEquals(mInv.preMultiply(m), identity); + TestUtils.assertEquals(m.preMultiply(identity), m); + TestUtils.assertEquals(identity.preMultiply(mInv), mInv); + try { + m.preMultiply(new DenseFieldMatrix(bigSingular)); + fail("Expecting illegalArgumentException"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + + public void testGetVectors() { + FieldMatrix m = new DenseFieldMatrix(testData); + TestUtils.assertEquals(m.getRow(0), testDataRow1); + TestUtils.assertEquals(m.getColumn(2), testDataCol3); + try { + m.getRow(10); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // ignored + } + try { + m.getColumn(-1); + fail("expecting MatrixIndexException"); + } catch (MatrixIndexException ex) { + // ignored + } + } + + public void testGetEntry() { + FieldMatrix m = new DenseFieldMatrix(testData); + assertEquals(m.getEntry(0,1),new Fraction(2)); + 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 + Fraction[][] matrixData = { + {new Fraction(1),new Fraction(2),new Fraction(3)}, + {new Fraction(2),new Fraction(5),new Fraction(3)} + }; + FieldMatrix m = new DenseFieldMatrix(matrixData); + // One more with three rows, two columns + Fraction[][] matrixData2 = { + {new Fraction(1),new Fraction(2)}, + {new Fraction(2),new Fraction(5)}, + {new Fraction(1), new Fraction(7)} + }; + FieldMatrix n = new DenseFieldMatrix(matrixData2); + // Now multiply m by n + FieldMatrix p = m.multiply(n); + assertEquals(2, p.getRowDimension()); + assertEquals(2, p.getColumnDimension()); + // Invert p + FieldMatrix pInverse = new FieldLUDecompositionImpl(p).getSolver().getInverse(); + assertEquals(2, pInverse.getRowDimension()); + assertEquals(2, pInverse.getColumnDimension()); + + // Solve example + Fraction[][] coefficientsData = { + {new Fraction(2), new Fraction(3), new Fraction(-2)}, + {new Fraction(-1), new Fraction(7), new Fraction(6)}, + {new Fraction(4), new Fraction(-3), new Fraction(-5)} + }; + FieldMatrix coefficients = new DenseFieldMatrix(coefficientsData); + Fraction[] constants = {new Fraction(1), new Fraction(-2), new Fraction(1)}; + Fraction[] solution = new FieldLUDecompositionImpl(coefficients).getSolver().solve(constants); + assertEquals(new Fraction(2).multiply(solution[0]). + add(new Fraction(3).multiply(solution[1])). + subtract(new Fraction(2).multiply(solution[2])), + constants[0]); + assertEquals(new Fraction(-1).multiply(solution[0]). + add(new Fraction(7).multiply(solution[1])). + add(new Fraction(6).multiply(solution[2])), + constants[1]); + assertEquals(new Fraction(4).multiply(solution[0]). + subtract(new Fraction(3).multiply(solution[1])). + subtract(new Fraction(5).multiply(solution[2])), + constants[2]); + + } + + // test submatrix accessors + public void testGetSubMatrix() { + FieldMatrix m = new DenseFieldMatrix(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(FieldMatrix m, Fraction[][] reference, + int startRow, int endRow, int startColumn, int endColumn, + boolean mustFail) { + try { + FieldMatrix sub = m.getSubMatrix(startRow, endRow, startColumn, endColumn); + assertEquals(new DenseFieldMatrix(reference), sub); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + private void checkGetSubMatrix(FieldMatrix m, Fraction[][] reference, + int[] selectedRows, int[] selectedColumns, + boolean mustFail) { + try { + FieldMatrix sub = m.getSubMatrix(selectedRows, selectedColumns); + assertEquals(new DenseFieldMatrix(reference), sub); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + public void testGetSetMatrixLarge() { + int n = 3 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = + new DenseFieldMatrix(FractionField.getInstance(), n, n); + FieldMatrix sub = + new DenseFieldMatrix(FractionField.getInstance(), n - 4, n - 4).scalarAdd(new Fraction(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(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + assertEquals(sub, m.getSubMatrix(2, n - 3, 2, n - 3)); + + } + + public void testCopySubMatrix() { + FieldMatrix m = new DenseFieldMatrix(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(FieldMatrix m, Fraction[][] reference, + int startRow, int endRow, int startColumn, int endColumn, + boolean mustFail) { + try { + Fraction[][] sub = (reference == null) ? + new Fraction[1][1] : + new Fraction[reference.length][reference[0].length]; + m.copySubMatrix(startRow, endRow, startColumn, endColumn, sub); + assertEquals(new DenseFieldMatrix(reference), new DenseFieldMatrix(sub)); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + private void checkCopy(FieldMatrix m, Fraction[][] reference, + int[] selectedRows, int[] selectedColumns, + boolean mustFail) { + try { + Fraction[][] sub = (reference == null) ? + new Fraction[1][1] : + new Fraction[reference.length][reference[0].length]; + m.copySubMatrix(selectedRows, selectedColumns, sub); + assertEquals(new DenseFieldMatrix(reference), new DenseFieldMatrix(sub)); + if (mustFail) { + fail("Expecting MatrixIndexException"); + } + } catch (MatrixIndexException e) { + if (!mustFail) { + throw e; + } + } + } + + public void testGetRowMatrix() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldMatrix mRow0 = new DenseFieldMatrix(subRow0); + FieldMatrix mRow3 = new DenseFieldMatrix(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() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldMatrix mRow3 = new DenseFieldMatrix(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 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = + new DenseFieldMatrix(FractionField.getInstance(), n, n); + FieldMatrix sub = + new DenseFieldMatrix(FractionField.getInstance(), 1, n).scalarAdd(new Fraction(1)); + + m.setRowMatrix(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + assertEquals(sub, m.getRowMatrix(2)); + + } + + public void testGetColumnMatrix() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldMatrix mColumn1 = new DenseFieldMatrix(subColumn1); + FieldMatrix mColumn3 = new DenseFieldMatrix(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() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldMatrix mColumn3 = new DenseFieldMatrix(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 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = + new DenseFieldMatrix(FractionField.getInstance(), n, n); + FieldMatrix sub = + new DenseFieldMatrix(FractionField.getInstance(), n, 1).scalarAdd(new Fraction(1)); + + m.setColumnMatrix(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + assertEquals(sub, m.getColumnMatrix(2)); + + } + + public void testGetRowVector() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldVector mRow0 = new FieldVectorImpl(subRow0[0]); + FieldVector mRow3 = new FieldVectorImpl(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() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldVector mRow3 = new FieldVectorImpl(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 FieldVectorImpl(FractionField.getInstance(), 5)); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetRowVectorLarge() { + int n = 3 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = new DenseFieldMatrix(FractionField.getInstance(), n, n); + FieldVector sub = new FieldVectorImpl(n, new Fraction(1)); + + m.setRowVector(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + assertEquals(sub, m.getRowVector(2)); + + } + + public void testGetColumnVector() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldVector mColumn1 = columnToVector(subColumn1); + FieldVector 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() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + FieldVector 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 FieldVectorImpl(FractionField.getInstance(), 5)); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetColumnVectorLarge() { + int n = 3 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = new DenseFieldMatrix(FractionField.getInstance(), n, n); + FieldVector sub = new FieldVectorImpl(n, new Fraction(1)); + + m.setColumnVector(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + assertEquals(sub, m.getColumnVector(2)); + + } + + private FieldVector columnToVector(Fraction[][] column) { + Fraction[] data = new Fraction[column.length]; + for (int i = 0; i < data.length; ++i) { + data[i] = column[i][0]; + } + return new FieldVectorImpl(data, false); + } + + public void testGetRow() { + FieldMatrix m = new DenseFieldMatrix(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() { + FieldMatrix m = new DenseFieldMatrix(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 Fraction[5]); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetRowLarge() { + int n = 3 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = new DenseFieldMatrix(FractionField.getInstance(), n, n); + Fraction[] sub = new Fraction[n]; + Arrays.fill(sub, new Fraction(1)); + + m.setRow(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (i != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + checkArrays(sub, m.getRow(2)); + + } + + public void testGetColumn() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + Fraction[] mColumn1 = columnToArray(subColumn1); + Fraction[] 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() { + FieldMatrix m = new DenseFieldMatrix(subTestData); + Fraction[] 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 Fraction[5]); + fail("Expecting InvalidMatrixException"); + } catch (InvalidMatrixException ex) { + // expected + } + } + + public void testGetSetColumnLarge() { + int n = 3 * DenseFieldMatrix.BLOCK_SIZE; + FieldMatrix m = new DenseFieldMatrix(FractionField.getInstance(), n, n); + Fraction[] sub = new Fraction[n]; + Arrays.fill(sub, new Fraction(1)); + + m.setColumn(2, sub); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (j != 2) { + assertEquals(new Fraction(0), m.getEntry(i, j)); + } else { + assertEquals(new Fraction(1), m.getEntry(i, j)); + } + } + } + checkArrays(sub, m.getColumn(2)); + + } + + private Fraction[] columnToArray(Fraction[][] column) { + Fraction[] data = new Fraction[column.length]; + for (int i = 0; i < data.length; ++i) { + data[i] = column[i][0]; + } + return data; + } + + private void checkArrays(Fraction[] expected, Fraction[] actual) { + assertEquals(expected.length, actual.length); + for (int i = 0; i < expected.length; ++i) { + assertEquals(expected[i], actual[i]); + } + } + + public void testEqualsAndHashCode() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + DenseFieldMatrix m1 = (DenseFieldMatrix) m.copy(); + DenseFieldMatrix mt = (DenseFieldMatrix) 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 DenseFieldMatrix(bigSingular))); + } + + public void testToString() { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + assertEquals("DenseFieldMatrix{{1,2,3},{2,5,3},{1,0,8}}", m.toString()); + } + + public void testSetSubMatrix() throws Exception { + DenseFieldMatrix m = new DenseFieldMatrix(testData); + m.setSubMatrix(detData2,1,1); + FieldMatrix expected = new DenseFieldMatrix + (new Fraction[][] {{new Fraction(1),new Fraction(2),new Fraction(3)},{new Fraction(2),new Fraction(1),new Fraction(3)},{new Fraction(1),new Fraction(2),new Fraction(4)}}); + assertEquals(expected, m); + + m.setSubMatrix(detData2,0,0); + expected = new DenseFieldMatrix + (new Fraction[][] {{new Fraction(1),new Fraction(3),new Fraction(3)},{new Fraction(2),new Fraction(4),new Fraction(3)},{new Fraction(1),new Fraction(2),new Fraction(4)}}); + assertEquals(expected, m); + + m.setSubMatrix(testDataPlus2,0,0); + expected = new DenseFieldMatrix + (new Fraction[][] {{new Fraction(3),new Fraction(4),new Fraction(5)},{new Fraction(4),new Fraction(7),new Fraction(5)},{new Fraction(3),new Fraction(2),new Fraction(10)}}); + assertEquals(expected, m); + + // javadoc example + DenseFieldMatrix matrix = + new DenseFieldMatrix(new Fraction[][] { + {new Fraction(1), new Fraction(2), new Fraction(3), new Fraction(4)}, + {new Fraction(5), new Fraction(6), new Fraction(7), new Fraction(8)}, + {new Fraction(9), new Fraction(0), new Fraction(1) , new Fraction(2)} + }); + matrix.setSubMatrix(new Fraction[][] { + {new Fraction(3), new Fraction(4)}, + {new Fraction(5), new Fraction(6)} + }, 1, 1); + expected = + new DenseFieldMatrix(new Fraction[][] { + {new Fraction(1), new Fraction(2), new Fraction(3),new Fraction(4)}, + {new Fraction(5), new Fraction(3), new Fraction(4), new Fraction(8)}, + {new Fraction(9), new Fraction(5) ,new Fraction(6), new Fraction(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 Fraction[][] {{new Fraction(1)}, {new Fraction(2), new Fraction(3)}}, 0, 0); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException e) { + // expected + } + + // empty + try { + m.setSubMatrix(new Fraction[][] {{}}, 0, 0); + fail("expecting IllegalArgumentException"); + } catch (IllegalArgumentException e) { + // expected + } + + } + + public void testWalk() { + int rows = 150; + int columns = 75; + + FieldMatrix m = new DenseFieldMatrix(FractionField.getInstance(), rows, columns); + m.walkInRowOrder(new SetVisitor()); + GetVisitor getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new DenseFieldMatrix(FractionField.getInstance(), 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(new Fraction(0), m.getEntry(i, 0)); + assertEquals(new Fraction(0), m.getEntry(i, columns - 1)); + } + for (int j = 0; j < columns; ++j) { + assertEquals(new Fraction(0), m.getEntry(0, j)); + assertEquals(new Fraction(0), m.getEntry(rows - 1, j)); + } + + m = new DenseFieldMatrix(FractionField.getInstance(), rows, columns); + m.walkInColumnOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInOptimizedOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new DenseFieldMatrix(FractionField.getInstance(), 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(new Fraction(0), m.getEntry(i, 0)); + assertEquals(new Fraction(0), m.getEntry(i, columns - 1)); + } + for (int j = 0; j < columns; ++j) { + assertEquals(new Fraction(0), m.getEntry(0, j)); + assertEquals(new Fraction(0), m.getEntry(rows - 1, j)); + } + + m = new DenseFieldMatrix(FractionField.getInstance(), rows, columns); + m.walkInOptimizedOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInRowOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new DenseFieldMatrix(FractionField.getInstance(), 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(new Fraction(0), m.getEntry(i, 0)); + assertEquals(new Fraction(0), m.getEntry(i, columns - 1)); + } + for (int j = 0; j < columns; ++j) { + assertEquals(new Fraction(0), m.getEntry(0, j)); + assertEquals(new Fraction(0), m.getEntry(rows - 1, j)); + } + + m = new DenseFieldMatrix(FractionField.getInstance(), rows, columns); + m.walkInOptimizedOrder(new SetVisitor()); + getVisitor = new GetVisitor(); + m.walkInColumnOrder(getVisitor); + assertEquals(rows * columns, getVisitor.getCount()); + + m = new DenseFieldMatrix(FractionField.getInstance(), 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(new Fraction(0), m.getEntry(i, 0)); + assertEquals(new Fraction(0), m.getEntry(i, columns - 1)); + } + for (int j = 0; j < columns; ++j) { + assertEquals(new Fraction(0), m.getEntry(0, j)); + assertEquals(new Fraction(0), m.getEntry(rows - 1, j)); + } + + } + + private static class SetVisitor extends DefaultFieldMatrixChangingVisitor { + private static final long serialVersionUID = 6586716383170553060L; + public SetVisitor() { + super(Fraction.ZERO); + } + @Override + public Fraction visit(int i, int j, Fraction value) { + return new Fraction(i * 11 + j, 11); + } + } + + private static class GetVisitor extends DefaultFieldMatrixPreservingVisitor { + private static final long serialVersionUID = 8363001284977267825L; + private int count; + public GetVisitor() { + super(Fraction.ZERO); + count = 0; + } + @Override + public void visit(int i, int j, Fraction value) { + ++count; + assertEquals(new Fraction(i * 11 + j, 11), value); + } + public int getCount() { + return count; + } + } + + private DenseFieldMatrix createRandomMatrix(Random r, int rows, int columns) { + DenseFieldMatrix m = + new DenseFieldMatrix(FractionField.getInstance(), rows, columns); + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < columns; ++j) { + int p = r.nextInt(20) - 10; + int q = r.nextInt(20) - 10; + if (q == 0) { + q = 1; + } + m.setEntry(i, j, new Fraction(p, q)); + } + } + return m; + } + +} +