http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java b/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java deleted file mode 100644 index cfbcb68..0000000 --- a/src/main/java/org/apache/commons/math3/linear/BiDiagonalTransformer.java +++ /dev/null @@ -1,380 +0,0 @@ -/* - * 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.math3.linear; - -import org.apache.commons.math3.util.FastMath; - - -/** - * Class transforming any matrix to bi-diagonal shape. - * <p>Any m × n matrix A can be written as the product of three matrices: - * A = U × B × V<sup>T</sup> with U an m × m orthogonal matrix, - * B an m × n bi-diagonal matrix (lower diagonal if m < n, upper diagonal - * otherwise), and V an n × n orthogonal matrix.</p> - * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is - * an intermediate step in more general decomposition algorithms like {@link - * SingularValueDecomposition Singular Value Decomposition}. This class is therefore - * intended for internal use by the library and is not public. As a consequence of - * this explicitly limited scope, many methods directly returns references to - * internal arrays, not copies.</p> - * @since 2.0 - */ -class BiDiagonalTransformer { - - /** Householder vectors. */ - private final double householderVectors[][]; - - /** Main diagonal. */ - private final double[] main; - - /** Secondary diagonal. */ - private final double[] secondary; - - /** Cached value of U. */ - private RealMatrix cachedU; - - /** Cached value of B. */ - private RealMatrix cachedB; - - /** Cached value of V. */ - private RealMatrix cachedV; - - /** - * Build the transformation to bi-diagonal shape of a matrix. - * @param matrix the matrix to transform. - */ - public BiDiagonalTransformer(RealMatrix matrix) { - - final int m = matrix.getRowDimension(); - final int n = matrix.getColumnDimension(); - final int p = FastMath.min(m, n); - householderVectors = matrix.getData(); - main = new double[p]; - secondary = new double[p - 1]; - cachedU = null; - cachedB = null; - cachedV = null; - - // transform matrix - if (m >= n) { - transformToUpperBiDiagonal(); - } else { - transformToLowerBiDiagonal(); - } - - } - - /** - * Returns the matrix U of the transform. - * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p> - * @return the U matrix - */ - public RealMatrix getU() { - - if (cachedU == null) { - - final int m = householderVectors.length; - final int n = householderVectors[0].length; - final int p = main.length; - final int diagOffset = (m >= n) ? 0 : 1; - final double[] diagonal = (m >= n) ? main : secondary; - double[][] ua = new double[m][m]; - - // fill up the part of the matrix not affected by Householder transforms - for (int k = m - 1; k >= p; --k) { - ua[k][k] = 1; - } - - // build up first part of the matrix by applying Householder transforms - for (int k = p - 1; k >= diagOffset; --k) { - final double[] hK = householderVectors[k]; - ua[k][k] = 1; - if (hK[k - diagOffset] != 0.0) { - for (int j = k; j < m; ++j) { - double alpha = 0; - for (int i = k; i < m; ++i) { - alpha -= ua[i][j] * householderVectors[i][k - diagOffset]; - } - alpha /= diagonal[k - diagOffset] * hK[k - diagOffset]; - - for (int i = k; i < m; ++i) { - ua[i][j] += -alpha * householderVectors[i][k - diagOffset]; - } - } - } - } - if (diagOffset > 0) { - ua[0][0] = 1; - } - cachedU = MatrixUtils.createRealMatrix(ua); - } - - // return the cached matrix - return cachedU; - - } - - /** - * Returns the bi-diagonal matrix B of the transform. - * @return the B matrix - */ - public RealMatrix getB() { - - if (cachedB == null) { - - final int m = householderVectors.length; - final int n = householderVectors[0].length; - double[][] ba = new double[m][n]; - for (int i = 0; i < main.length; ++i) { - ba[i][i] = main[i]; - if (m < n) { - if (i > 0) { - ba[i][i-1] = secondary[i - 1]; - } - } else { - if (i < main.length - 1) { - ba[i][i+1] = secondary[i]; - } - } - } - cachedB = MatrixUtils.createRealMatrix(ba); - } - - // return the cached matrix - return cachedB; - - } - - /** - * Returns the matrix V of the transform. - * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p> - * @return the V matrix - */ - public RealMatrix getV() { - - if (cachedV == null) { - - final int m = householderVectors.length; - final int n = householderVectors[0].length; - final int p = main.length; - final int diagOffset = (m >= n) ? 1 : 0; - final double[] diagonal = (m >= n) ? secondary : main; - double[][] va = new double[n][n]; - - // fill up the part of the matrix not affected by Householder transforms - for (int k = n - 1; k >= p; --k) { - va[k][k] = 1; - } - - // build up first part of the matrix by applying Householder transforms - for (int k = p - 1; k >= diagOffset; --k) { - final double[] hK = householderVectors[k - diagOffset]; - va[k][k] = 1; - if (hK[k] != 0.0) { - for (int j = k; j < n; ++j) { - double beta = 0; - for (int i = k; i < n; ++i) { - beta -= va[i][j] * hK[i]; - } - beta /= diagonal[k - diagOffset] * hK[k]; - - for (int i = k; i < n; ++i) { - va[i][j] += -beta * hK[i]; - } - } - } - } - if (diagOffset > 0) { - va[0][0] = 1; - } - cachedV = MatrixUtils.createRealMatrix(va); - } - - // return the cached matrix - return cachedV; - - } - - /** - * Get the Householder vectors of the transform. - * <p>Note that since this class is only intended for internal use, - * it returns directly a reference to its internal arrays, not a copy.</p> - * @return the main diagonal elements of the B matrix - */ - double[][] getHouseholderVectorsRef() { - return householderVectors; - } - - /** - * Get the main diagonal elements of the matrix B of the transform. - * <p>Note that since this class is only intended for internal use, - * it returns directly a reference to its internal arrays, not a copy.</p> - * @return the main diagonal elements of the B matrix - */ - double[] getMainDiagonalRef() { - return main; - } - - /** - * Get the secondary diagonal elements of the matrix B of the transform. - * <p>Note that since this class is only intended for internal use, - * it returns directly a reference to its internal arrays, not a copy.</p> - * @return the secondary diagonal elements of the B matrix - */ - double[] getSecondaryDiagonalRef() { - return secondary; - } - - /** - * Check if the matrix is transformed to upper bi-diagonal. - * @return true if the matrix is transformed to upper bi-diagonal - */ - boolean isUpperBiDiagonal() { - return householderVectors.length >= householderVectors[0].length; - } - - /** - * Transform original matrix to upper bi-diagonal form. - * <p>Transformation is done using alternate Householder transforms - * on columns and rows.</p> - */ - private void transformToUpperBiDiagonal() { - - final int m = householderVectors.length; - final int n = householderVectors[0].length; - for (int k = 0; k < n; k++) { - - //zero-out a column - double xNormSqr = 0; - for (int i = k; i < m; ++i) { - final double c = householderVectors[i][k]; - xNormSqr += c * c; - } - final double[] hK = householderVectors[k]; - final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); - main[k] = a; - if (a != 0.0) { - hK[k] -= a; - for (int j = k + 1; j < n; ++j) { - double alpha = 0; - for (int i = k; i < m; ++i) { - final double[] hI = householderVectors[i]; - alpha -= hI[j] * hI[k]; - } - alpha /= a * householderVectors[k][k]; - for (int i = k; i < m; ++i) { - final double[] hI = householderVectors[i]; - hI[j] -= alpha * hI[k]; - } - } - } - - if (k < n - 1) { - //zero-out a row - xNormSqr = 0; - for (int j = k + 1; j < n; ++j) { - final double c = hK[j]; - xNormSqr += c * c; - } - final double b = (hK[k + 1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); - secondary[k] = b; - if (b != 0.0) { - hK[k + 1] -= b; - for (int i = k + 1; i < m; ++i) { - final double[] hI = householderVectors[i]; - double beta = 0; - for (int j = k + 1; j < n; ++j) { - beta -= hI[j] * hK[j]; - } - beta /= b * hK[k + 1]; - for (int j = k + 1; j < n; ++j) { - hI[j] -= beta * hK[j]; - } - } - } - } - - } - } - - /** - * Transform original matrix to lower bi-diagonal form. - * <p>Transformation is done using alternate Householder transforms - * on rows and columns.</p> - */ - private void transformToLowerBiDiagonal() { - - final int m = householderVectors.length; - final int n = householderVectors[0].length; - for (int k = 0; k < m; k++) { - - //zero-out a row - final double[] hK = householderVectors[k]; - double xNormSqr = 0; - for (int j = k; j < n; ++j) { - final double c = hK[j]; - xNormSqr += c * c; - } - final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); - main[k] = a; - if (a != 0.0) { - hK[k] -= a; - for (int i = k + 1; i < m; ++i) { - final double[] hI = householderVectors[i]; - double alpha = 0; - for (int j = k; j < n; ++j) { - alpha -= hI[j] * hK[j]; - } - alpha /= a * householderVectors[k][k]; - for (int j = k; j < n; ++j) { - hI[j] -= alpha * hK[j]; - } - } - } - - if (k < m - 1) { - //zero-out a column - final double[] hKp1 = householderVectors[k + 1]; - xNormSqr = 0; - for (int i = k + 1; i < m; ++i) { - final double c = householderVectors[i][k]; - xNormSqr += c * c; - } - final double b = (hKp1[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); - secondary[k] = b; - if (b != 0.0) { - hKp1[k] -= b; - for (int j = k + 1; j < n; ++j) { - double beta = 0; - for (int i = k + 1; i < m; ++i) { - final double[] hI = householderVectors[i]; - beta -= hI[j] * hI[k]; - } - beta /= b * hKp1[k]; - for (int i = k + 1; i < m; ++i) { - final double[] hI = householderVectors[i]; - hI[j] -= beta * hI[k]; - } - } - } - } - - } - } - -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/BlockFieldMatrix.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/BlockFieldMatrix.java b/src/main/java/org/apache/commons/math3/linear/BlockFieldMatrix.java deleted file mode 100644 index 4a41f7f..0000000 --- a/src/main/java/org/apache/commons/math3/linear/BlockFieldMatrix.java +++ /dev/null @@ -1,1592 +0,0 @@ -/* - * 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.math3.linear; - -import java.io.Serializable; - -import org.apache.commons.math3.Field; -import org.apache.commons.math3.FieldElement; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.util.MathUtils; - -/** - * Cache-friendly implementation of FieldMatrix using a flat arrays to store - * square blocks of the matrix. - * <p> - * 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. - * </p> - * <p> - * 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. - * </p> - * <p> - * 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}<sup>2</sup> elements long for regular blocks. The blocks are themselves - * organized in row major order. - * </p> - * <p> - * 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. - * </p> - * <p> - * 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. - * </p> - * @param <T> the type of the field elements - * @since 2.0 - */ -public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMatrix<T> implements Serializable { - /** Block size. */ - public static final int BLOCK_SIZE = 36; - /** Serializable version identifier. */ - private static final long serialVersionUID = -4602336630143123183L; - /** 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 Number of rows in the new matrix. - * @param columns Number of columns in the new matrix. - * @throws NotStrictlyPositiveException if row or column dimension is not - * positive. - */ - public BlockFieldMatrix(final Field<T> field, final int rows, - final int columns) - throws NotStrictlyPositiveException { - 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. - * <p>The input array <em>must</em> already be in raw layout.</p> - * <p>Calling this constructor is equivalent to call: - * <pre>matrix = new BlockFieldMatrix<T>(getField(), rawData.length, rawData[0].length, - * toBlocksLayout(rawData), false);</pre> - * </p> - * - * @param rawData Data for the new matrix, in raw layout. - * @throws DimensionMismatchException if the {@code blockData} shape is - * inconsistent with block layout. - * @see #BlockFieldMatrix(int, int, FieldElement[][], boolean) - */ - public BlockFieldMatrix(final T[][] rawData) - throws DimensionMismatchException { - this(rawData.length, rawData[0].length, toBlocksLayout(rawData), false); - } - - /** - * Create a new dense matrix copying entries from block layout data. - * <p>The input array <em>must</em> already be in blocks layout.</p> - * @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 - * - * @throws DimensionMismatchException if the {@code blockData} shape is - * inconsistent with block layout. - * @throws NotStrictlyPositiveException if row or column dimension is not - * positive. - * @see #createBlocksLayout(Field, int, int) - * @see #toBlocksLayout(FieldElement[][]) - * @see #BlockFieldMatrix(FieldElement[][]) - */ - public BlockFieldMatrix(final int rows, final int columns, - final T[][] blockData, final boolean copyArray) - throws DimensionMismatchException, NotStrictlyPositiveException { - 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 = MathArrays.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 new DimensionMismatchException(blockData[index].length, - iHeight * blockWidth(jBlock)); - } - if (copyArray) { - blocks[index] = blockData[index].clone(); - } - } - } - } - - /** - * Convert a data array from raw layout to blocks layout. - * <p> - * Raw layout is the straightforward layout where element at row i and - * column j is in array element <code>rawData[i][j]</code>. Blocks layout - * is the layout used in {@link BlockFieldMatrix} 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. - * </p> - * <p> - * 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 - * #BlockFieldMatrix(int, int, FieldElement[][], boolean)} - * constructor. - * </p> - * @param <T> Type of the field elements. - * @param rawData Data array in raw layout. - * @return a new data array containing the same entries but in blocks layout - * @throws DimensionMismatchException if {@code rawData} is not rectangular - * (not all rows have the same length). - * @see #createBlocksLayout(Field, int, int) - * @see #BlockFieldMatrix(int, int, FieldElement[][], boolean) - */ - public static <T extends FieldElement<T>> T[][] toBlocksLayout(final T[][] rawData) - throws DimensionMismatchException { - - 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 new DimensionMismatchException(columns, length); - } - } - - // convert array - final Field<T> field = extractField(rawData); - final T[][] blocks = MathArrays.buildArray(field, blockRows * blockColumns, -1); - int blockIndex = 0; - for (int iBlock = 0; iBlock < blockRows; ++iBlock) { - final int pStart = iBlock * BLOCK_SIZE; - final int pEnd = FastMath.min(pStart + BLOCK_SIZE, rows); - final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - final int jWidth = qEnd - qStart; - - // allocate new block - final T[] block = MathArrays.buildArray(field, iHeight * jWidth); - blocks[blockIndex] = block; - - // copy data - int index = 0; - for (int p = pStart; p < pEnd; ++p) { - System.arraycopy(rawData[p], qStart, block, index, jWidth); - index += jWidth; - } - - ++blockIndex; - } - } - - return blocks; - } - - /** - * Create a data array in blocks layout. - * <p> - * This method can be used to create the array argument of the {@link - * #BlockFieldMatrix(int, int, FieldElement[][], boolean)} - * constructor. - * </p> - * @param <T> Type of the field elements. - * @param field Field to which the elements belong. - * @param rows Number of rows in the new matrix. - * @param columns Number of columns in the new matrix. - * @return a new data array in blocks layout. - * @see #toBlocksLayout(FieldElement[][]) - * @see #BlockFieldMatrix(int, int, FieldElement[][], boolean) - */ - public static <T extends FieldElement<T>> T[][] createBlocksLayout(final Field<T> 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 = MathArrays.buildArray(field, blockRows * blockColumns, -1); - int blockIndex = 0; - for (int iBlock = 0; iBlock < blockRows; ++iBlock) { - final int pStart = iBlock * BLOCK_SIZE; - final int pEnd = FastMath.min(pStart + BLOCK_SIZE, rows); - final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - final int jWidth = qEnd - qStart; - blocks[blockIndex] = MathArrays.buildArray(field, iHeight * jWidth); - ++blockIndex; - } - } - - return blocks; - } - - /** {@inheritDoc} */ - @Override - public FieldMatrix<T> createMatrix(final int rowDimension, - final int columnDimension) - throws NotStrictlyPositiveException { - return new BlockFieldMatrix<T>(getField(), rowDimension, - columnDimension); - } - - /** {@inheritDoc} */ - @Override - public FieldMatrix<T> copy() { - - // create an empty matrix - BlockFieldMatrix<T> copied = new BlockFieldMatrix<T>(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<T> add(final FieldMatrix<T> m) - throws MatrixDimensionMismatchException { - try { - return add((BlockFieldMatrix<T>) m); - } catch (ClassCastException cce) { - - // safety check - checkAdditionCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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 = FastMath.min(pStart + BLOCK_SIZE, rows); - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q) { - outBlock[k] = tBlock[k].add(m.getEntry(p, q)); - ++k; - } - } - - // go to next block - ++blockIndex; - - } - } - - return out; - } - } - - /** - * Compute the sum of {@code this} and {@code m}. - * - * @param m matrix to be added - * @return {@code this + m} - * @throws MatrixDimensionMismatchException if {@code m} is not the same - * size as {@code this} - */ - public BlockFieldMatrix<T> add(final BlockFieldMatrix<T> m) - throws MatrixDimensionMismatchException { - - // safety check - checkAdditionCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> subtract(final FieldMatrix<T> m) - throws MatrixDimensionMismatchException { - try { - return subtract((BlockFieldMatrix<T>) m); - } catch (ClassCastException cce) { - - // safety check - checkSubtractionCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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 = FastMath.min(pStart + BLOCK_SIZE, rows); - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q) { - outBlock[k] = tBlock[k].subtract(m.getEntry(p, q)); - ++k; - } - } - - // go to next block - ++blockIndex; - - } - } - - return out; - } - } - - /** - * Compute {@code this - m}. - * - * @param m matrix to be subtracted - * @return {@code this - m} - * @throws MatrixDimensionMismatchException if {@code m} is not the same - * size as {@code this} - */ - public BlockFieldMatrix<T> subtract(final BlockFieldMatrix<T> m) throws MatrixDimensionMismatchException { - // safety check - checkSubtractionCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> scalarAdd(final T d) { - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> scalarMultiply(final T d) { - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> multiply(final FieldMatrix<T> m) - throws DimensionMismatchException { - try { - return multiply((BlockFieldMatrix<T>) m); - } catch (ClassCastException cce) { - - // safety check - checkMultiplicationCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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 = FastMath.min(pStart + BLOCK_SIZE, rows); - - for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { - - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.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; - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - final int lStart = (p - pStart) * kWidth; - final int lEnd = lStart + kWidth; - for (int q = qStart; q < qEnd; ++q) { - T sum = zero; - int r = rStart; - for (int l = lStart; l < lEnd; ++l) { - sum = sum.add(tBlock[l].multiply(m.getEntry(r, q))); - ++r; - } - outBlock[k] = outBlock[k].add(sum); - ++k; - } - } - } - - // go to next block - ++blockIndex; - - } - } - - return out; - } - } - - /** - * Returns the result of postmultiplying {@code this} by {@code m}. - * - * @param m matrix to postmultiply by - * @return {@code this * m} - * @throws DimensionMismatchException if the matrices are not compatible. - */ - public BlockFieldMatrix<T> multiply(BlockFieldMatrix<T> m) - throws DimensionMismatchException { - - // safety check - checkMultiplicationCompatible(m); - - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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 = FastMath.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]; - int k = 0; - for (int p = pStart; 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 = MathArrays.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 = FastMath.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<T> getSubMatrix(final int startRow, final int endRow, - final int startColumn, - final int endColumn) - throws OutOfRangeException, NumberIsTooSmallException { - // safety checks - checkSubMatrixIndex(startRow, endRow, startColumn, endColumn); - - // create the output matrix - final BlockFieldMatrix<T> out = - new BlockFieldMatrix<T>(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 - int pBlock = blockStartRow; - for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { - final int iHeight = out.blockHeight(iBlock); - int qBlock = blockStartColumn; - for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { - 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); - } - } - ++qBlock; - } - ++pBlock; - } - - return out; - } - - /** - * Copy a part of a block into another one - * <p>This method can be called only when the specified part fits in both - * blocks, no verification is done here.</p> - * @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 DimensionMismatchException, OutOfRangeException, - NoDataException, NullArgumentException { - // safety checks - MathUtils.checkNotNull(subMatrix); - final int refLength = subMatrix[0].length; - if (refLength == 0) { - throw new NoDataException(LocalizedFormats.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 new DimensionMismatchException(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 = FastMath.max(row, firstRow); - final int iEnd = FastMath.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 = FastMath.max(column, firstColumn); - final int jEnd = FastMath.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<T> getRowMatrix(final int row) - throws OutOfRangeException { - checkRowIndex(row); - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> matrix) - throws MatrixDimensionMismatchException, OutOfRangeException { - try { - setRowMatrix(row, (BlockFieldMatrix<T>) matrix); - } catch (ClassCastException cce) { - super.setRowMatrix(row, matrix); - } - } - - /** - * Sets the entries in row number <code>row</code> - * 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 MatrixDimensionMismatchException if the matrix dimensions do - * not match one instance row. - * @throws OutOfRangeException if the specified row index is invalid. - */ - public void setRowMatrix(final int row, final BlockFieldMatrix<T> matrix) - throws MatrixDimensionMismatchException, OutOfRangeException { - checkRowIndex(row); - final int nCols = getColumnDimension(); - if ((matrix.getRowDimension() != 1) || - (matrix.getColumnDimension() != nCols)) { - throw new MatrixDimensionMismatchException(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<T> getColumnMatrix(final int column) - throws OutOfRangeException { - checkColumnIndex(column); - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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<T> matrix) - throws MatrixDimensionMismatchException, OutOfRangeException { - try { - setColumnMatrix(column, (BlockFieldMatrix<T>) matrix); - } catch (ClassCastException cce) { - super.setColumnMatrix(column, matrix); - } - } - - /** - * Sets the entries in column number {@code column} - * as a column matrix. Column indices start at 0. - * - * @param column Column to be set. - * @param matrix Column matrix (must have one column and the same number of rows - * as the instance). - * @throws MatrixDimensionMismatchException if the matrix dimensions do - * not match one instance column. - * @throws OutOfRangeException if the specified column index is invalid. - */ - void setColumnMatrix(final int column, final BlockFieldMatrix<T> matrix) - throws MatrixDimensionMismatchException, OutOfRangeException { - checkColumnIndex(column); - final int nRows = getRowDimension(); - if ((matrix.getRowDimension() != nRows) || - (matrix.getColumnDimension() != 1)) { - throw new MatrixDimensionMismatchException(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<T> getRowVector(final int row) - throws OutOfRangeException { - checkRowIndex(row); - final T[] outData = MathArrays.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 ArrayFieldVector<T>(getField(), outData, false); - } - - /** {@inheritDoc} */ - @Override - public void setRowVector(final int row, final FieldVector<T> vector) - throws MatrixDimensionMismatchException, OutOfRangeException { - try { - setRow(row, ((ArrayFieldVector<T>) vector).getDataRef()); - } catch (ClassCastException cce) { - super.setRowVector(row, vector); - } - } - - /** {@inheritDoc} */ - @Override - public FieldVector<T> getColumnVector(final int column) - throws OutOfRangeException { - checkColumnIndex(column); - final T[] outData = MathArrays.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 ArrayFieldVector<T>(getField(), outData, false); - } - - /** {@inheritDoc} */ - @Override - public void setColumnVector(final int column, final FieldVector<T> vector) - throws OutOfRangeException, MatrixDimensionMismatchException { - try { - setColumn(column, ((ArrayFieldVector<T>) vector).getDataRef()); - } catch (ClassCastException cce) { - super.setColumnVector(column, vector); - } - } - - /** {@inheritDoc} */ - @Override - public T[] getRow(final int row) throws OutOfRangeException { - checkRowIndex(row); - final T[] out = MathArrays.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 OutOfRangeException, MatrixDimensionMismatchException { - checkRowIndex(row); - final int nCols = getColumnDimension(); - if (array.length != nCols) { - throw new MatrixDimensionMismatchException(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 OutOfRangeException { - checkColumnIndex(column); - final T[] out = MathArrays.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 MatrixDimensionMismatchException, OutOfRangeException { - checkColumnIndex(column); - final int nRows = getRowDimension(); - if (array.length != nRows) { - throw new MatrixDimensionMismatchException(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 OutOfRangeException { - checkRowIndex(row); - checkColumnIndex(column); - - 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]; - } - - /** {@inheritDoc} */ - @Override - public void setEntry(final int row, final int column, final T value) - throws OutOfRangeException { - checkRowIndex(row); - checkColumnIndex(column); - - 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; - } - - /** {@inheritDoc} */ - @Override - public void addToEntry(final int row, final int column, final T increment) - throws OutOfRangeException { - checkRowIndex(row); - checkColumnIndex(column); - - 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); - } - - /** {@inheritDoc} */ - @Override - public void multiplyEntry(final int row, final int column, final T factor) - throws OutOfRangeException { - checkRowIndex(row); - checkColumnIndex(column); - - 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); - } - - /** {@inheritDoc} */ - @Override - public FieldMatrix<T> transpose() { - final int nRows = getRowDimension(); - final int nCols = getColumnDimension(); - final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(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 = FastMath.min(pStart + BLOCK_SIZE, columns); - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, rows); - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - final int lInc = pEnd - pStart; - int l = p - pStart; - for (int q = qStart; q < qEnd; ++q) { - outBlock[k] = tBlock[l]; - ++k; - l+= lInc; - } - } - - // 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 DimensionMismatchException { - if (v.length != columns) { - throw new DimensionMismatchException(v.length, columns); - } - final T[] out = MathArrays.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 = FastMath.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 = FastMath.min(qStart + BLOCK_SIZE, columns); - int k = 0; - for (int p = pStart; 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 DimensionMismatchException { - - if (v.length != rows) { - throw new DimensionMismatchException(v.length, rows); - } - final T[] out = MathArrays.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 = FastMath.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 = FastMath.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<T> visitor) { - 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 = FastMath.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 = FastMath.min(qStart + BLOCK_SIZE, columns); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - int k = (p - pStart) * jWidth; - for (int q = qStart; q < qEnd; ++q) { - block[k] = visitor.visit(p, q, block[k]); - ++k; - } - } - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInRowOrder(final FieldMatrixPreservingVisitor<T> visitor) { - 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 = FastMath.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 = FastMath.min(qStart + BLOCK_SIZE, columns); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - int k = (p - pStart) * jWidth; - for (int q = qStart; q < qEnd; ++q) { - visitor.visit(p, q, block[k]); - ++k; - } - } - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInRowOrder(final FieldMatrixChangingVisitor<T> visitor, - final int startRow, final int endRow, - final int startColumn, final int endColumn) - throws OutOfRangeException, NumberIsTooSmallException { - 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 = FastMath.max(startRow, p0); - final int pEnd = FastMath.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 = FastMath.max(startColumn, q0); - final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - int k = (p - p0) * jWidth + qStart - q0; - for (int q = qStart; q < qEnd; ++q) { - block[k] = visitor.visit(p, q, block[k]); - ++k; - } - } - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInRowOrder(final FieldMatrixPreservingVisitor<T> visitor, - final int startRow, final int endRow, - final int startColumn, final int endColumn) - throws OutOfRangeException, NumberIsTooSmallException { - 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 = FastMath.max(startRow, p0); - final int pEnd = FastMath.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 = FastMath.max(startColumn, q0); - final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - int k = (p - p0) * jWidth + qStart - q0; - for (int q = qStart; q < qEnd; ++q) { - visitor.visit(p, q, block[k]); - ++k; - } - } - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> visitor) { - visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - int blockIndex = 0; - for (int iBlock = 0; iBlock < blockRows; ++iBlock) { - final int pStart = iBlock * BLOCK_SIZE; - final int pEnd = FastMath.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - final T[] block = blocks[blockIndex]; - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q) { - block[k] = visitor.visit(p, q, block[k]); - ++k; - } - } - ++blockIndex; - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> visitor) { - visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - int blockIndex = 0; - for (int iBlock = 0; iBlock < blockRows; ++iBlock) { - final int pStart = iBlock * BLOCK_SIZE; - final int pEnd = FastMath.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { - final int qStart = jBlock * BLOCK_SIZE; - final int qEnd = FastMath.min(qStart + BLOCK_SIZE, columns); - final T[] block = blocks[blockIndex]; - int k = 0; - for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q) { - visitor.visit(p, q, block[k]); - ++k; - } - } - ++blockIndex; - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> visitor, - final int startRow, final int endRow, - final int startColumn, final int endColumn) - throws OutOfRangeException, NumberIsTooSmallException { - 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 = FastMath.max(startRow, p0); - final int pEnd = FastMath.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 = FastMath.max(startColumn, q0); - final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int p = pStart; p < pEnd; ++p) { - int k = (p - p0) * jWidth + qStart - q0; - for (int q = qStart; q < qEnd; ++q) { - block[k] = visitor.visit(p, q, block[k]); - ++k; - } - } - } - } - return visitor.end(); - } - - /** {@inheritDoc} */ - @Override - public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> visitor, - final int startRow, final int endRow, - final int startColumn, final int endColumn) - throws OutOfRangeException, NumberIsTooSmallException { - 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 = FastMath.max(startRow, p0); - final int pEnd = FastMath.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 = FastMath.max(startColumn, q0); - final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); - final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int p = pStart; p < pEnd; ++p) { - int k = (p - p0) * jWidth + qStart - q0; - for (int q = qStart; q < qEnd; ++q) { - visitor.visit(p, q, block[k]); - ++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; - } -}
