You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/04/20 22:19:49 UTC

svn commit: r766852 [1/2] - in /commons/proper/math/trunk/src: java/org/apache/commons/math/linear/DenseFieldMatrix.java java/org/apache/commons/math/linear/MatrixUtils.java test/org/apache/commons/math/linear/DenseFieldMatrixTest.java

Author: luc
Date: Mon Apr 20 20:19:49 2009
New Revision: 766852

URL: http://svn.apache.org/viewvc?rev=766852&view=rev
Log:
added DenseFieldMatrix

Added:
    commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java   (with props)
    commons/proper/math/trunk/src/test/org/apache/commons/math/linear/DenseFieldMatrixTest.java   (with props)
Modified:
    commons/proper/math/trunk/src/java/org/apache/commons/math/linear/MatrixUtils.java

Added: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java?rev=766852&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java Mon Apr 20 20:19:49 2009
@@ -0,0 +1,1613 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.linear;
+
+import org.apache.commons.math.Field;
+import org.apache.commons.math.FieldElement;
+import org.apache.commons.math.MathRuntimeException;
+
+/**
+ * Cache-friendly implementation of FieldMatrix using a flat arrays to store
+ * square blocks of the matrix.
+ * <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>
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class DenseFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMatrix<T> {
+    
+    /** Serializable version identifier */
+    private static final long serialVersionUID = -4602336630143123183L;
+
+    /** Block size. */
+    public static final int BLOCK_SIZE = 36;
+
+    /** Blocks of matrix entries. */
+    private final T blocks[][];
+
+    /** Number of rows of the matrix. */
+    private final int rows;
+
+    /** Number of columns of the matrix. */
+    private final int columns;
+
+    /** Number of block rows of the matrix. */
+    private final int blockRows;
+
+    /** Number of block columns of the matrix. */
+    private final int blockColumns;
+
+    /**
+     * Create a new matrix with the supplied row and column dimensions.
+     *
+     * @param field field to which the elements belong
+     * @param rows  the number of rows in the new matrix
+     * @param columns  the number of columns in the new matrix
+     * @throws IllegalArgumentException if row or column dimension is not
+     *  positive
+     */
+    public DenseFieldMatrix(final Field<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 DenseFieldMatrix<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 #DenseFieldMatrix(int, int, T[][], boolean)
+     */
+    public DenseFieldMatrix(final T[][] rawData)
+        throws IllegalArgumentException {
+        this(rawData.length, rawData[0].length, toBlocksLayout(rawData), false);
+    }
+
+    /**
+     * Create a new dense matrix copying entries from block layout data.
+     * <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(int, int)
+     * @see #toBlocksLayout(T[][])
+     * @see #DenseFieldMatrix(T[][])
+     */
+    public DenseFieldMatrix(final int rows, final int columns,
+                           final T[][] blockData, final boolean copyArray)
+        throws IllegalArgumentException {
+
+        super(extractField(blockData), rows, columns);
+        this.rows    = rows;
+        this.columns = columns;
+
+        // number of blocks
+        blockRows    = (rows    + BLOCK_SIZE - 1) / BLOCK_SIZE;
+        blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
+
+        if (copyArray) {
+            // allocate storage blocks, taking care of smaller ones at right and bottom
+            blocks = buildArray(getField(), blockRows * blockColumns, -1);
+        } else {
+            // reference existing array
+            blocks = blockData;
+        }
+
+        int index = 0;
+        for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
+            final int iHeight = blockHeight(iBlock);
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++index) {
+                if (blockData[index].length != iHeight * blockWidth(jBlock)) {
+                    throw MathRuntimeException.createIllegalArgumentException(
+                            "wrong array shape (block length = {0}, expected {1})",
+                            blockData[index].length, iHeight * blockWidth(jBlock));
+                }
+                if (copyArray) {
+                    blocks[index] = blockData[index].clone();
+                }
+            }
+        }
+
+    }
+
+    /**
+     * Convert a data array from raw layout to blocks layout.
+     * <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 DenseFieldMatrix} instances, where the matrix
+     * is split in square blocks (except at right and bottom side where blocks may
+     * be rectangular to fit matrix size) and each block is stored in a flattened
+     * one-dimensional array.
+     * </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
+     * DenseFieldMatrix#DenseFieldMatrix(int, int, T[][], boolean)} constructor.
+     * </p>
+     * @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(int, int)
+     * @see #DenseFieldMatrix(int, int, T[][], 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);
+        for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart  = iBlock * BLOCK_SIZE;
+            final int pEnd    = Math.min(pStart + BLOCK_SIZE, rows);
+            final int iHeight = pEnd - pStart;
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
+                final int qStart = jBlock * BLOCK_SIZE;
+                final int qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+                final int jWidth = qEnd - qStart;
+
+                // allocate new block
+                final T[] block = buildArray(field, iHeight * jWidth);
+                blocks[blockIndex] = block;
+
+                // copy data
+                for (int p = pStart, index = 0; p < pEnd; ++p, index += jWidth) {
+                    System.arraycopy(rawData[p], qStart, block, index, jWidth);
+                }
+
+            }
+        }
+
+        return blocks;
+
+    }
+
+    /**
+     * Create a data array in blocks layout.
+     * <p>
+     * This method can be used to create the array argument of the {@link
+     * DenseFieldMatrix#DenseFieldMatrix(int, int, T[][], boolean)} constructor.
+     * </p>
+     * @param rows  the number of rows in the new matrix
+     * @param columns  the number of columns in the new matrix
+     * @return a new data array in blocks layout
+     * @see #toBlocksLayout(T[][])
+     * @see #DenseFieldMatrix(int, int, T[][], boolean)
+     */
+    public static <T 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);
+        for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart  = iBlock * BLOCK_SIZE;
+            final int pEnd    = Math.min(pStart + BLOCK_SIZE, rows);
+            final int iHeight = pEnd - pStart;
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
+                final int qStart = jBlock * BLOCK_SIZE;
+                final int qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+                final int jWidth = qEnd - qStart;
+                blocks[blockIndex] = buildArray(field, iHeight * jWidth);
+            }
+        }
+
+        return blocks;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public FieldMatrix<T> createMatrix(final int rowDimension, final int columnDimension)
+        throws IllegalArgumentException {
+        return new DenseFieldMatrix<T>(getField(), rowDimension, columnDimension);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public FieldMatrix<T> copy() {
+
+        // create an empty matrix
+        DenseFieldMatrix<T> copied = new DenseFieldMatrix<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((DenseFieldMatrix<T>) m);
+        } catch (ClassCastException cce) {
+
+            // safety check
+            checkAdditionCompatible(m);
+
+            final DenseFieldMatrix<T> out = new DenseFieldMatrix<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);
+                    for (int p = pStart, k = 0; p < pEnd; ++p) {
+                        for (int q = qStart; q < qEnd; ++q, ++k) {
+                            outBlock[k] = tBlock[k].add(m.getEntry(p, q));
+                        }
+                    }
+
+                    // go to next block
+                    ++blockIndex;
+
+                }
+            }
+
+            return out;
+
+        }
+    }
+
+    /**
+     * Compute the sum of this and <code>m</code>.
+     *
+     * @param m    matrix to be added
+     * @return     this + m
+     * @throws  IllegalArgumentException if m is not the same size as this
+     */
+    public DenseFieldMatrix<T> add(final DenseFieldMatrix<T> m)
+        throws IllegalArgumentException {
+
+        // safety check
+        checkAdditionCompatible(m);
+
+        final DenseFieldMatrix<T> out = new DenseFieldMatrix<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((DenseFieldMatrix<T>) m);
+        } catch (ClassCastException cce) {
+
+            // safety check
+            checkSubtractionCompatible(m);
+
+            final DenseFieldMatrix<T> out = new DenseFieldMatrix<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);
+                    for (int p = pStart, k = 0; p < pEnd; ++p) {
+                        for (int q = qStart; q < qEnd; ++q, ++k) {
+                            outBlock[k] = tBlock[k].subtract(m.getEntry(p, q));
+                        }
+                    }
+
+                    // go to next block
+                    ++blockIndex;
+
+                }
+            }
+
+            return out;
+
+        }
+    }
+
+    /**
+     * Compute this minus <code>m</code>.
+     *
+     * @param m    matrix to be subtracted
+     * @return     this - m
+     * @throws  IllegalArgumentException if m is not the same size as this
+     */
+    public DenseFieldMatrix<T> subtract(final DenseFieldMatrix<T> m)
+        throws IllegalArgumentException {
+
+        // safety check
+        checkSubtractionCompatible(m);
+
+        final DenseFieldMatrix<T> out = new DenseFieldMatrix<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 DenseFieldMatrix<T> out = new DenseFieldMatrix<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 DenseFieldMatrix<T> out = new DenseFieldMatrix<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((DenseFieldMatrix<T>) m);
+        } catch (ClassCastException cce) {
+
+            // safety check
+            checkMultiplicationCompatible(m);
+
+            final DenseFieldMatrix<T> out = new DenseFieldMatrix<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;
+                        for (int p = pStart, k = 0; p < pEnd; ++p) {
+                            final int lStart = (p - pStart) * kWidth;
+                            final int lEnd   = lStart + kWidth;
+                            for (int q = qStart; q < qEnd; ++q) {
+                                T sum = zero;
+                                for (int l = lStart, r = rStart; l < lEnd; ++l, ++r) {
+                                    sum = sum.add(tBlock[l].multiply(m.getEntry(r, q)));
+                                }
+                                outBlock[k] = outBlock[k].add(sum);
+                                ++k;
+                            }
+                        }
+                    }
+
+                    // go to next block
+                    ++blockIndex;
+
+                }
+            }
+
+            return out;
+
+        }
+    }
+
+    /**
+     * Returns the result of postmultiplying this by m.
+     *
+     * @param m    matrix to postmultiply by
+     * @return     this * m
+     * @throws     IllegalArgumentException
+     *             if columnDimension(this) != rowDimension(m)
+     */
+    public DenseFieldMatrix<T> multiply(DenseFieldMatrix<T> m) throws IllegalArgumentException {
+
+        // safety check
+        checkMultiplicationCompatible(m);
+
+        final DenseFieldMatrix<T> out = new DenseFieldMatrix<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];
+                    for (int p = pStart, k = 0; p < pEnd; ++p) {
+                        final int lStart = (p - pStart) * kWidth;
+                        final int lEnd   = lStart + kWidth;
+                        for (int nStart = 0; nStart < jWidth; ++nStart) {
+                            T sum = zero;
+                            int l = lStart;
+                            int n = nStart;
+                            while (l < lEnd - 3) {
+                                sum = sum.
+                                      add(tBlock[l].multiply(mBlock[n])).
+                                      add(tBlock[l + 1].multiply(mBlock[n + jWidth])).
+                                      add(tBlock[l + 2].multiply(mBlock[n + jWidth2])).
+                                      add(tBlock[l + 3].multiply(mBlock[n + jWidth3]));
+                                l += 4;
+                                n += jWidth4;
+                            }
+                            while (l < lEnd) {
+                                sum = sum.add(tBlock[l++].multiply(mBlock[n]));
+                                n += jWidth;
+                            }
+                            outBlock[k] = outBlock[k].add(sum);
+                            ++k;
+                        }
+                    }
+                }
+
+                // go to next block
+                ++blockIndex;
+
+            }
+        }
+
+        return out;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T[][] getData() {
+
+        final T[][] data = buildArray(getField(), getRowDimension(), getColumnDimension());
+        final int lastColumns = columns - (blockColumns - 1) * BLOCK_SIZE;
+
+        for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart = iBlock * BLOCK_SIZE;
+            final int pEnd   = Math.min(pStart + BLOCK_SIZE, rows);
+            int regularPos   = 0;
+            int lastPos      = 0;
+            for (int p = pStart; p < pEnd; ++p) {
+                final T[] dataP = data[p];
+                int blockIndex = iBlock * blockColumns;
+                int dataPos    = 0;
+                for (int jBlock = 0; jBlock < blockColumns - 1; ++jBlock) {
+                    System.arraycopy(blocks[blockIndex++], regularPos, dataP, dataPos, BLOCK_SIZE);
+                    dataPos += BLOCK_SIZE;
+                }
+                System.arraycopy(blocks[blockIndex], lastPos, dataP, dataPos, lastColumns);
+                regularPos += BLOCK_SIZE;
+                lastPos    += lastColumns;
+            }
+        }
+
+        return data;
+        
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public FieldMatrix<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 DenseFieldMatrix<T> out =
+            new DenseFieldMatrix<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
+        for (int iBlock = 0, pBlock = blockStartRow; iBlock < out.blockRows; ++iBlock, ++pBlock) {
+            final int iHeight = out.blockHeight(iBlock);
+            for (int jBlock = 0, qBlock = blockStartColumn; jBlock < out.blockColumns; ++jBlock, ++qBlock) {
+                final int jWidth = out.blockWidth(jBlock);
+
+                // handle one block of the output matrix
+                final int      outIndex = iBlock * out.blockColumns + jBlock;
+                final T[] outBlock = out.blocks[outIndex];
+                final int      index    = pBlock * blockColumns + qBlock;
+                final int      width    = blockWidth(qBlock);
+
+                final int heightExcess = iHeight + rowsShift - BLOCK_SIZE;
+                final int widthExcess  = jWidth + columnsShift - BLOCK_SIZE;
+                if (heightExcess > 0) {
+                    // the submatrix block spans on two blocks rows from the original matrix
+                    if (widthExcess > 0) {
+                        // the submatrix block spans on two blocks columns from the original matrix
+                        final int width2 = blockWidth(qBlock + 1);
+                        copyBlockPart(blocks[index], width,
+                                      rowsShift, BLOCK_SIZE,
+                                      columnsShift, BLOCK_SIZE,
+                                      outBlock, jWidth, 0, 0);
+                        copyBlockPart(blocks[index + 1], width2,
+                                      rowsShift, BLOCK_SIZE,
+                                      0, widthExcess,
+                                      outBlock, jWidth, 0, jWidth - widthExcess);
+                        copyBlockPart(blocks[index + blockColumns], width,
+                                      0, heightExcess,
+                                      columnsShift, BLOCK_SIZE,
+                                      outBlock, jWidth, iHeight - heightExcess, 0);
+                        copyBlockPart(blocks[index + blockColumns + 1], width2,
+                                      0, heightExcess,
+                                      0, widthExcess,
+                                      outBlock, jWidth, iHeight - heightExcess, jWidth - widthExcess);
+                    } else {
+                        // the submatrix block spans on one block column from the original matrix
+                        copyBlockPart(blocks[index], width,
+                                      rowsShift, BLOCK_SIZE,
+                                      columnsShift, jWidth + columnsShift,
+                                      outBlock, jWidth, 0, 0);
+                        copyBlockPart(blocks[index + blockColumns], width,
+                                      0, heightExcess,
+                                      columnsShift, jWidth + columnsShift,
+                                      outBlock, jWidth, iHeight - heightExcess, 0);
+                    }
+                } else {
+                    // the submatrix block spans on one block row from the original matrix
+                    if (widthExcess > 0) {
+                        // the submatrix block spans on two blocks columns from the original matrix
+                        final int width2 = blockWidth(qBlock + 1);
+                        copyBlockPart(blocks[index], width,
+                                      rowsShift, iHeight + rowsShift,
+                                      columnsShift, BLOCK_SIZE,
+                                      outBlock, jWidth, 0, 0);
+                        copyBlockPart(blocks[index + 1], width2,
+                                      rowsShift, iHeight + rowsShift,
+                                      0, widthExcess,
+                                      outBlock, jWidth, 0, jWidth - widthExcess);
+                    } else {
+                        // the submatrix block spans on one block column from the original matrix
+                        copyBlockPart(blocks[index], width,
+                                      rowsShift, iHeight + rowsShift,
+                                      columnsShift, jWidth + columnsShift,
+                                      outBlock, jWidth, 0, 0);
+                    }
+               }
+
+            }
+        }
+
+        return out;
+
+    }
+
+    /**
+     * Copy a part of a block into another one
+     * <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 DenseFieldMatrix<T> out = new DenseFieldMatrix<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, (DenseFieldMatrix<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 DenseFieldMatrix<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 DenseFieldMatrix<T> out = new DenseFieldMatrix<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, (DenseFieldMatrix<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 DenseFieldMatrix<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 FieldVectorImpl<T>(outData, false);
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void setRowVector(final int row, final FieldVector<T> vector)
+        throws MatrixIndexException, InvalidMatrixException {
+        try {
+            setRow(row, ((FieldVectorImpl<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 FieldVectorImpl<T>(outData, false);
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void setColumnVector(final int column, final FieldVector<T> vector)
+        throws MatrixIndexException, InvalidMatrixException {
+        try {
+            setColumn(column, ((FieldVectorImpl<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 DenseFieldMatrix<T> out = new DenseFieldMatrix<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);
+                for (int p = pStart, k = 0; p < pEnd; ++p) {
+                    final int lInc = pEnd - pStart;
+                    for (int q = qStart, l = p - pStart; q < qEnd; ++q, l+= lInc) {
+                        outBlock[k++] = tBlock[l];
+                    }
+                }
+
+                // go to next block
+                ++blockIndex;
+
+            }
+        }
+
+        return out;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getRowDimension() {
+        return rows;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getColumnDimension() {
+        return columns;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T[] operate(final T[] v)
+        throws IllegalArgumentException {
+
+        if (v.length != columns) {
+            throw MathRuntimeException.createIllegalArgumentException(
+                    "vector length mismatch: got {0} but expected {1}",
+                    v.length, columns);
+        }
+        final T[] out = buildArray(getField(), rows);
+        final T zero = getField().getZero();
+
+        // perform multiplication block-wise, to ensure good cache behavior
+        for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart = iBlock * BLOCK_SIZE;
+            final int pEnd   = Math.min(pStart + BLOCK_SIZE, rows);
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
+                final T[] block  = blocks[iBlock * blockColumns + jBlock];
+                final int      qStart = jBlock * BLOCK_SIZE;
+                final int      qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+                for (int p = pStart, k = 0; p < pEnd; ++p) {
+                    T sum = zero;
+                    int q = qStart;
+                    while (q < qEnd - 3) {
+                        sum = sum.
+                              add(block[k].multiply(v[q])).
+                              add(block[k + 1].multiply(v[q + 1])).
+                              add(block[k + 2].multiply(v[q + 2])).
+                              add(block[k + 3].multiply(v[q + 3]));
+                        k += 4;
+                        q += 4;
+                    }
+                    while (q < qEnd) {
+                        sum = sum.add(block[k++].multiply(v[q++]));
+                    }
+                    out[p] = out[p].add(sum);
+                }
+            }
+        }
+
+        return out;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T[] preMultiply(final T[] v)
+        throws IllegalArgumentException {
+
+        if (v.length != rows) {
+            throw MathRuntimeException.createIllegalArgumentException(
+                    "vector length mismatch: got {0} but expected {1}",
+                    v.length, rows);
+        }
+        final T[] out = buildArray(getField(), columns);
+        final T zero = getField().getZero();
+
+        // perform multiplication block-wise, to ensure good cache behavior
+        for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
+            final int jWidth  = blockWidth(jBlock);
+            final int jWidth2 = jWidth  + jWidth;
+            final int jWidth3 = jWidth2 + jWidth;
+            final int jWidth4 = jWidth3 + jWidth;
+            final int qStart = jBlock * BLOCK_SIZE;
+            final int qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
+                final T[] block  = blocks[iBlock * blockColumns + jBlock];
+                final int      pStart = iBlock * BLOCK_SIZE;
+                final int      pEnd   = Math.min(pStart + BLOCK_SIZE, rows);
+                for (int q = qStart; q < qEnd; ++q) {
+                    int k = q - qStart;
+                    T sum = zero;
+                    int p = pStart;
+                    while (p < pEnd - 3) {
+                        sum = sum.
+                              add(block[k].multiply(v[p])).
+                              add(block[k + jWidth].multiply(v[p + 1])).
+                              add(block[k + jWidth2].multiply(v[p + 2])).
+                              add(block[k + jWidth3].multiply(v[p + 3]));
+                        k += jWidth4;
+                        p += 4;
+                    }
+                    while (p < pEnd) {
+                        sum = sum.add(block[k].multiply(v[p++]));
+                        k += jWidth;
+                    }
+                    out[q] = out[q].add(sum);
+                }
+            }
+        }
+
+        return out;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInRowOrder(final FieldMatrixChangingVisitor<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];
+                    for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) {
+                        block[k] = visitor.visit(p, q, block[k]);
+                    }
+                }
+             }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInRowOrder(final FieldMatrixPreservingVisitor<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];
+                    for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) {
+                        visitor.visit(p, q, block[k]);
+                    }
+                }
+             }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInRowOrder(final FieldMatrixChangingVisitor<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];
+                    for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
+                        block[k] = visitor.visit(p, q, block[k]);
+                    }
+                }
+             }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInRowOrder(final FieldMatrixPreservingVisitor<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];
+                    for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
+                        visitor.visit(p, q, block[k]);
+                    }
+                }
+             }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> visitor)
+        throws MatrixVisitorException {
+        visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
+        for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart = iBlock * BLOCK_SIZE;
+            final int pEnd   = Math.min(pStart + BLOCK_SIZE, rows);
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
+                final int qStart = jBlock * BLOCK_SIZE;
+                final int qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+                final T[] block = blocks[blockIndex];
+                for (int p = pStart, k = 0; p < pEnd; ++p) {
+                    for (int q = qStart; q < qEnd; ++q, ++k) {
+                        block[k] = visitor.visit(p, q, block[k]);
+                    }
+                }
+            }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> visitor)
+        throws MatrixVisitorException {
+        visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
+        for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
+            final int pStart = iBlock * BLOCK_SIZE;
+            final int pEnd   = Math.min(pStart + BLOCK_SIZE, rows);
+            for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
+                final int qStart = jBlock * BLOCK_SIZE;
+                final int qEnd   = Math.min(qStart + BLOCK_SIZE, columns);
+                final T[] block = blocks[blockIndex];
+                for (int p = pStart, k = 0; p < pEnd; ++p) {
+                    for (int q = qStart; q < qEnd; ++q, ++k) {
+                        visitor.visit(p, q, block[k]);
+                    }
+                }
+            }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<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) {
+                    for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
+                        block[k] = visitor.visit(p, q, block[k]);
+                    }
+                }
+            }
+        }
+        return visitor.end();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<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) {
+                    for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
+                        visitor.visit(p, q, block[k]);
+                    }
+                }
+            }
+        }
+        return visitor.end();
+    }
+
+    /**
+     * Get the height of a block.
+     * @param blockRow row index (in block sense) of the block
+     * @return height (number of rows) of the block
+     */
+    private int blockHeight(final int blockRow) {
+        return (blockRow == blockRows - 1) ? rows - blockRow * BLOCK_SIZE : BLOCK_SIZE;
+    }
+
+    /**
+     * Get the width of a block.
+     * @param blockColumn column index (in block sense) of the block
+     * @return width (number of columns) of the block
+     */
+    private int blockWidth(final int blockColumn) {
+        return (blockColumn == blockColumns - 1) ? columns - blockColumn * BLOCK_SIZE : BLOCK_SIZE;
+    }
+
+}

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/DenseFieldMatrix.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/MatrixUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/MatrixUtils.java?rev=766852&r1=766851&r2=766852&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/linear/MatrixUtils.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/linear/MatrixUtils.java Mon Apr 20 20:19:49 2009
@@ -17,9 +17,13 @@
 
 package org.apache.commons.math.linear;
 
+import java.lang.reflect.Array;
 import java.math.BigDecimal;
 import java.util.Arrays;
 
+import org.apache.commons.math.Field;
+import org.apache.commons.math.FieldElement;
+
 /**
  * A collection of static methods that operate on or return matrices.
  * 
@@ -47,6 +51,22 @@
     }
 
     /**
+     * Returns a {@link FieldMatrix} with specified dimensions.
+     * <p>The matrix elements are all set to field.getZero().</p>
+     * @param field field to which the matrix elements belong
+     * @param rows number of rows of the matrix
+     * @param columns number of columns of the matrix
+     * @return  FieldMatrix with specified dimensions
+     * @see #createFieldMatrix(FieldElement[][])
+     * @since 2.0
+     */
+    public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(final Field<T> field,
+                                                                               final int rows,
+                                                                               final int columns) {
+        return new DenseFieldMatrix<T>(field, rows, columns);
+    }
+
+    /**
      * Returns a {@link RealMatrix} whose entries are the the values in the
      * the input array.  The input array is copied, not referenced.
      * 
@@ -62,6 +82,24 @@
     }
 
     /**
+     * Returns a {@link FieldMatrix} whose entries are the the values in the
+     * the input array.
+     * <p>
+     * The input array is copied, not referenced.
+     * </p>
+     * @param data input array
+     * @return  RealMatrix containing the values of the array
+     * @throws IllegalArgumentException if <code>data</code> is not rectangular
+     *  (not all rows have the same length) or empty
+     * @throws NullPointerException if <code>data</code> is null
+     * @see #createFieldMatrix(Field, int, int)
+     * @since 2.0
+     */
+    public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(T[][] data) {
+        return new DenseFieldMatrix<T>(data);
+    }
+
+    /**
      * Returns <code>dimension x dimension</code> identity matrix.
      *
      * @param dimension dimension of identity matrix to generate
@@ -76,6 +114,48 @@
         }
         return m;
     }
+
+    /**
+     * Returns <code>dimension x dimension</code> identity matrix.
+     *
+     * @param dimension dimension of identity matrix to generate
+     * @return identity matrix
+     * @throws IllegalArgumentException if dimension is not positive
+     * @since 2.0
+     */
+    @SuppressWarnings("unchecked")
+    public static <T extends FieldElement<T>> FieldMatrix<T>
+        createFieldIdentityMatrix(final Field<T> field, final int dimension) {
+        final T zero = field.getZero();
+        final T one  = field.getOne();
+        final T[][] d = (T[][]) Array.newInstance(zero.getClass(), dimension, dimension);
+        for (int row = 0; row < dimension; row++) {
+            final T[] dRow = d[row];
+            Arrays.fill(dRow, zero);
+            dRow[row] = one;
+        }
+        return new FieldMatrixImpl<T>(d, false);
+    }
+
+    /**
+     * Returns <code>dimension x dimension</code> identity matrix.
+     *
+     * @param dimension dimension of identity matrix to generate
+     * @return identity matrix
+     * @throws IllegalArgumentException if dimension is not positive
+     * @since 1.1
+     * @deprecated since 2.0, replaced by {@link #createFieldIdentityMatrix(Field, int)}
+     */
+    @Deprecated
+    public static BigMatrix createBigIdentityMatrix(int dimension) {
+        final BigDecimal[][] d = new BigDecimal[dimension][dimension];
+        for (int row = 0; row < dimension; row++) {
+            final BigDecimal[] dRow = d[row];
+            Arrays.fill(dRow, BigMatrixImpl.ZERO);
+            dRow[row] = BigMatrixImpl.ONE;
+        }
+        return new BigMatrixImpl(d, false);
+    }
     
     /**
      * Returns a diagonal matrix with specified elements.
@@ -94,6 +174,24 @@
     }
     
     /**
+     * Returns a diagonal matrix with specified elements.
+     *
+     * @param diagonal diagonal elements of the matrix (the array elements
+     * will be copied)
+     * @return diagonal matrix
+     * @since 2.0
+     */
+    public static <T extends FieldElement<T>> FieldMatrix<T>
+        createFieldDiagonalMatrix(final T[] diagonal) {
+        final FieldMatrix<T> m =
+            createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length);
+        for (int i = 0; i < diagonal.length; ++i) {
+            m.setEntry(i, i, diagonal[i]);
+        }
+        return m;
+    }
+    
+    /**
      * Returns a {@link BigMatrix} whose entries are the the values in the
      * the input array.  The input array is copied, not referenced.
      * 
@@ -102,7 +200,9 @@
      * @throws IllegalArgumentException if <code>data</code> is not rectangular
      *  (not all rows have the same length) or empty
      * @throws NullPointerException if data is null
+     * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])}
      */
+    @Deprecated
     public static BigMatrix createBigMatrix(double[][] data) {
         return new BigMatrixImpl(data);
     }
@@ -116,7 +216,9 @@
      * @throws IllegalArgumentException if <code>data</code> is not rectangular
      *  (not all rows have the same length) or empty
      * @throws NullPointerException if data is null
+     * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])}
      */
+    @Deprecated
     public static BigMatrix createBigMatrix(BigDecimal[][] data) {
         return new BigMatrixImpl(data);
     }
@@ -136,7 +238,9 @@
      *  (not all rows have the same length) or empty
      * @throws NullPointerException if <code>data</code> is null
      * @see #createRealMatrix(double[][])
+     * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])}
      */
+    @Deprecated
     public static BigMatrix createBigMatrix(BigDecimal[][] data, boolean copyArray) {
         return new BigMatrixImpl(data, copyArray);
     }
@@ -150,7 +254,9 @@
      * @throws IllegalArgumentException if <code>data</code> is not rectangular
      *  (not all rows have the same length) or empty
      * @throws NullPointerException if data is null
+     * @deprecated since 2.0 replaced by {@link #createFieldMatrix(FieldElement[][])}
      */
+    @Deprecated
     public static BigMatrix createBigMatrix(String[][] data) {
         return new BigMatrixImpl(data);
     }
@@ -168,6 +274,18 @@
     }
     
     /**
+     * Creates a {@link FieldVector} using the data from the input array. 
+     * 
+     * @param data the input data
+     * @return a data.length FieldVector
+     * @throws IllegalArgumentException if <code>data</code> is empty
+     * @throws NullPointerException if <code>data</code>is null
+     */
+    public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final T[] data) {
+        return new FieldVectorImpl<T>(data, true);
+    }
+    
+    /**
      * Creates a row {@link RealMatrix} using the data from the input
      * array. 
      * 
@@ -186,6 +304,25 @@
     }
     
     /**
+     * Creates a row {@link FieldMatrix} using the data from the input
+     * array. 
+     * 
+     * @param rowData the input row data
+     * @return a 1 x rowData.length FieldMatrix
+     * @throws IllegalArgumentException if <code>rowData</code> is empty
+     * @throws NullPointerException if <code>rowData</code>is null
+     */
+    public static <T extends FieldElement<T>> FieldMatrix<T>
+        createRowFieldMatrix(final T[] rowData) {
+        final int nCols = rowData.length;
+        final FieldMatrix<T> m = createFieldMatrix(rowData[0].getField(), 1, nCols);
+        for (int i = 0; i < nCols; ++i) {
+            m.setEntry(0, i, rowData[i]);
+        }
+        return m;
+    }
+    
+    /**
      * Creates a row {@link BigMatrix} using the data from the input
      * array. 
      * 
@@ -193,7 +330,9 @@
      * @return a 1 x rowData.length BigMatrix
      * @throws IllegalArgumentException if <code>rowData</code> is empty
      * @throws NullPointerException if <code>rowData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createRowBigMatrix(double[] rowData) {
         final int nCols = rowData.length;
         final BigDecimal[][] data = new BigDecimal[1][nCols];
@@ -211,7 +350,9 @@
      * @return a 1 x rowData.length BigMatrix
      * @throws IllegalArgumentException if <code>rowData</code> is empty
      * @throws NullPointerException if <code>rowData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createRowBigMatrix(BigDecimal[] rowData) {
         final int nCols = rowData.length;
         final BigDecimal[][] data = new BigDecimal[1][nCols];
@@ -227,7 +368,9 @@
      * @return a 1 x rowData.length BigMatrix
      * @throws IllegalArgumentException if <code>rowData</code> is empty
      * @throws NullPointerException if <code>rowData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createRowFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createRowBigMatrix(String[] rowData) {
         final int nCols = rowData.length;
         final BigDecimal[][] data = new BigDecimal[1][nCols];
@@ -256,6 +399,25 @@
     }
     
     /**
+     * Creates a column {@link FieldMatrix} using the data from the input
+     * array.
+     * 
+     * @param columnData  the input column data
+     * @return a columnData x 1 FieldMatrix
+     * @throws IllegalArgumentException if <code>columnData</code> is empty
+     * @throws NullPointerException if <code>columnData</code>is null
+     */
+    public static <T extends FieldElement<T>> FieldMatrix<T>
+        createColumnFieldMatrix(final T[] columnData) {
+        final int nRows = columnData.length;
+        final FieldMatrix<T> m = createFieldMatrix(columnData[0].getField(), nRows, 1);
+        for (int i = 0; i < nRows; ++i) {
+            m.setEntry(i, 0, columnData[i]);
+        }
+        return m;
+    }
+    
+    /**
      * Creates a column {@link BigMatrix} using the data from the input
      * array.
      * 
@@ -263,7 +425,9 @@
      * @return a columnData x 1 BigMatrix
      * @throws IllegalArgumentException if <code>columnData</code> is empty
      * @throws NullPointerException if <code>columnData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createColumnBigMatrix(double[] columnData) {
         final int nRows = columnData.length;
         final BigDecimal[][] data = new BigDecimal[nRows][1];
@@ -281,7 +445,9 @@
      * @return a columnData x 1 BigMatrix
      * @throws IllegalArgumentException if <code>columnData</code> is empty
      * @throws NullPointerException if <code>columnData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createColumnBigMatrix(BigDecimal[] columnData) {
         final int nRows = columnData.length;
         final BigDecimal[][] data = new BigDecimal[nRows][1];
@@ -299,7 +465,9 @@
      * @return a columnData x 1 BigMatrix
      * @throws IllegalArgumentException if <code>columnData</code> is empty
      * @throws NullPointerException if <code>columnData</code>is null
+     * @deprecated since 2.0 replaced by {@link #createColumnFieldMatrix(FieldElement[])}
      */
+    @Deprecated
     public static BigMatrix createColumnBigMatrix(String[] columnData) {
         int nRows = columnData.length;
         final BigDecimal[][] data = new BigDecimal[nRows][1];
@@ -309,23 +477,5 @@
         return new BigMatrixImpl(data, false);
     }
     
-    /**
-     * Returns <code>dimension x dimension</code> identity matrix.
-     *
-     * @param dimension dimension of identity matrix to generate
-     * @return identity matrix
-     * @throws IllegalArgumentException if dimension is not positive
-     * @since 1.1
-     */
-    public static BigMatrix createBigIdentityMatrix(int dimension) {
-        final BigDecimal[][] d = new BigDecimal[dimension][dimension];
-        for (int row = 0; row < dimension; row++) {
-            final BigDecimal[] dRow = d[row];
-            Arrays.fill(dRow, BigMatrixImpl.ZERO);
-            dRow[row] = BigMatrixImpl.ONE;
-        }
-        return new BigMatrixImpl(d, false);
-    }
-    
 }