/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  package org.apache.commons.math.linear;
  
  import java.io.Serializable;
  
  import org.apache.commons.math.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.
   * <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
   * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
   * @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  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 BlockFieldMatrix(final Field<T> 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.
      * <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 new matrix, in raw layout
      *
      * @exception IllegalArgumentException if <code>blockData</code> shape is
      * inconsistent with block layout
      * @see #BlockFieldMatrix(int, int, FieldElement[][], boolean)
      */
     public BlockFieldMatrix(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.
      * <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
      *
      * @exception IllegalArgumentException if <code>blockData</code> shape is
      * inconsistent with block layout
      * @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 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.
      * <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> the 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
      * @exception IllegalArgumentException if <code>rawData</code> 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 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<T> field = extractField(rawData);
         final T[][] blocks = buildArray(field, blockRows * blockColumns, -1);
         int blockIndex = 0;
         for (int iBlock = 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) {
                 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
                 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> the type of the field elements
      * @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
      * @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 = buildArray(field, blockRows * blockColumns, -1);
         int blockIndex = 0;
         for (int iBlock = 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) {
                 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);
                 ++blockIndex;
             }
         }
 
         return blocks;
 
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldMatrix<T> createMatrix(final int rowDimension, final int columnDimension)
         throws IllegalArgumentException {
         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 IllegalArgumentException {
         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     = Math.min(pStart + BLOCK_SIZE, rows);
                     final int      qStart   = jBlock * BLOCK_SIZE;
                     final int      qEnd     = Math.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 this and <code>m</code>.
      *
      * @param m    matrix to be added
      * @return     this + m
      * @throws  IllegalArgumentException if m is not the same size as this
      */
     public BlockFieldMatrix<T> add(final BlockFieldMatrix<T> m)
         throws IllegalArgumentException {
 
         // 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 IllegalArgumentException {
         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     = Math.min(pStart + BLOCK_SIZE, rows);
                     final int      qStart   = jBlock * BLOCK_SIZE;
                     final int      qEnd     = Math.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 this minus <code>m</code>.
      *
      * @param m    matrix to be subtracted
      * @return     this - m
      * @throws  IllegalArgumentException if m is not the same size as this
      */
     public BlockFieldMatrix<T> subtract(final BlockFieldMatrix<T> m)
         throws IllegalArgumentException {
 
         // 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)
         throws IllegalArgumentException {
 
         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)
         throws IllegalArgumentException {
 
         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 IllegalArgumentException {
         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   = 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;
                         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 this by m.
      *
      * @param m    matrix to postmultiply by
      * @return     this * m
      * @throws     IllegalArgumentException
      *             if columnDimension(this) != rowDimension(m)
      */
     public BlockFieldMatrix<T> multiply(BlockFieldMatrix<T> m) throws IllegalArgumentException {
 
         // 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   = 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];
                     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 = 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<T> 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 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 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<T> getRowMatrix(final int row)
         throws MatrixIndexException {
 
         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 MatrixIndexException, InvalidMatrixException {
         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 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 BlockFieldMatrix<T> 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<T> getColumnMatrix(final int column)
         throws MatrixIndexException {
 
         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 MatrixIndexException, InvalidMatrixException {
         try {
             setColumnMatrix(column, (BlockFieldMatrix<T>) matrix);
         } catch (ClassCastException cce) {
             super.setColumnMatrix(column, matrix);
         }
     }
 
     /**
      * Sets the entries in column number <code>column</code>
      * 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 BlockFieldMatrix<T> 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<T> 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 ArrayFieldVector<T>(outData, false);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public void setRowVector(final int row, final FieldVector<T> vector)
         throws MatrixIndexException, InvalidMatrixException {
         try {
             setRow(row, ((ArrayFieldVector<T>) vector).getDataRef());
         } catch (ClassCastException cce) {
             super.setRowVector(row, vector);
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldVector<T> 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 ArrayFieldVector<T>(outData, false);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public void setColumnVector(final int column, final FieldVector<T> vector)
         throws MatrixIndexException, InvalidMatrixException {
         try {
             setColumn(column, ((ArrayFieldVector<T>) 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<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     = Math.min(pStart + BLOCK_SIZE, columns);
                 final int      qStart   = jBlock * BLOCK_SIZE;
                 final int      qEnd     = Math.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 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);
                 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 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<T> 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];
                     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)
         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];
                     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 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];
                     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 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];
                     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)
         throws MatrixVisitorException {
         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   = Math.min(pStart + BLOCK_SIZE, rows);
             for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
                 final int qStart = jBlock * BLOCK_SIZE;
                 final int qEnd   = Math.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)
         throws MatrixVisitorException {
         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   = Math.min(pStart + BLOCK_SIZE, rows);
             for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
                 final int qStart = jBlock * BLOCK_SIZE;
                 final int qEnd   = Math.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 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) {
                     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 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) {
                     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;
     }
 
 }