    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math.linear;
    019
    020    import java.io.Serializable;
    021    import java.util.Arrays;
    022
    023    import org.apache.commons.math.MathRuntimeException;
    024    import org.apache.commons.math.linear.MatrixVisitorException;
    025    import org.apache.commons.math.exception.util.LocalizedFormats;
    026    import org.apache.commons.math.util.FastMath;
    027
    028    /**
    029     * Cache-friendly implementation of RealMatrix using a flat arrays to store
    030     * square blocks of the matrix.
    031     * <p>
    032     * This implementation is specially designed to be cache-friendly. Square blocks are
    033     * stored as small arrays and allow efficient traversal of data both in row major direction
    034     * and columns major direction, one block at a time. This greatly increases performances
    035     * for algorithms that use crossed directions loops like multiplication or transposition.
    036     * </p>
    037     * <p>
    038     * The size of square blocks is a static parameter. It may be tuned according to the cache
    039     * size of the target computer processor. As a rule of thumbs, it should be the largest
    040     * value that allows three blocks to be simultaneously cached (this is necessary for example
    041     * for matrix multiplication). The default value is to use 52x52 blocks which is well suited
    042     * for processors with 64k L1 cache (one block holds 2704 values or 21632 bytes). This value
    043     * could be lowered to 36x36 for processors with 32k L1 cache.
    044     * </p>
    045     * <p>
    046     * The regular blocks represent {@link #BLOCK_SIZE} x {@link #BLOCK_SIZE} squares. Blocks
    047     * at right hand side and bottom side which may be smaller to fit matrix dimensions. The square
    048     * blocks are flattened in row major order in single dimension arrays which are therefore
    049     * {@link #BLOCK_SIZE}<sup>2</sup> elements long for regular blocks. The blocks are themselves
    050     * organized in row major order.
    051     * </p>
    052     * <p>
    053     * As an example, for a block size of 52x52, a 100x60 matrix would be stored in 4 blocks.
    054     * Block 0 would be a double[2704] array holding the upper left 52x52 square, block 1 would be
    055     * a double[416] array holding the upper right 52x8 rectangle, block 2 would be a double[2496]
    056     * array holding the lower left 48x52 rectangle and block 3 would be a double[384] array
    057     * holding the lower right 48x8 rectangle.
    058     * </p>
    059     * <p>
    060     * The layout complexity overhead versus simple mapping of matrices to java
    061     * arrays is negligible for small matrices (about 1%). The gain from cache efficiency leads
    062     * to up to 3-fold improvements for matrices of moderate to large size.
    063     * </p>
    064     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    065     * @since 2.0
    066     */
    067    public class BlockRealMatrix extends AbstractRealMatrix implements Serializable {
    068
    069        /** Block size. */
    070        public static final int BLOCK_SIZE = 52;
    071
    072        /** Serializable version identifier */
    073        private static final long serialVersionUID = 4991895511313664478L;
    074
    075        /** Blocks of matrix entries. */
    076        private final double blocks[][];
    077
    078        /** Number of rows of the matrix. */
    079        private final int rows;
    080
    081        /** Number of columns of the matrix. */
    082        private final int columns;
    083
    084        /** Number of block rows of the matrix. */
    085        private final int blockRows;
    086
    087        /** Number of block columns of the matrix. */
    088        private final int blockColumns;
    089
    090        /**
    091         * Create a new matrix with the supplied row and column dimensions.
    092         *
    093         * @param rows  the number of rows in the new matrix
    094         * @param columns  the number of columns in the new matrix
    095         * @throws IllegalArgumentException if row or column dimension is not
    096         *  positive
    097         */
    098        public BlockRealMatrix(final int rows, final int columns)
    099            throws IllegalArgumentException {
    100
    101            super(rows, columns);
    102            this.rows    = rows;
    103            this.columns = columns;
    104
    105            // number of blocks
    106            blockRows    = (rows    + BLOCK_SIZE - 1) / BLOCK_SIZE;
    107            blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
    108
    109            // allocate storage blocks, taking care of smaller ones at right and bottom
    110            blocks = createBlocksLayout(rows, columns);
    111
    112        }
    113
    114        /**
    115         * Create a new dense matrix copying entries from raw layout data.
    116         * <p>The input array <em>must</em> already be in raw layout.</p>
    117         * <p>Calling this constructor is equivalent to call:
    118         * <pre>matrix = new BlockRealMatrix(rawData.length, rawData[0].length,
    119         *                                   toBlocksLayout(rawData), false);</pre>
    120         * </p>
    121         * @param rawData data for new matrix, in raw layout
    122         *
    123         * @exception IllegalArgumentException if <code>blockData</code> shape is
    124         * inconsistent with block layout
    125         * @see #BlockRealMatrix(int, int, double[][], boolean)
    126         */
    127        public BlockRealMatrix(final double[][] rawData)
    128            throws IllegalArgumentException {
    129            this(rawData.length, rawData[0].length, toBlocksLayout(rawData), false);
    130        }
    131
    132        /**
    133         * Create a new dense matrix copying entries from block layout data.
    134         * <p>The input array <em>must</em> already be in blocks layout.</p>
    135         * @param rows  the number of rows in the new matrix
    136         * @param columns  the number of columns in the new matrix
    137         * @param blockData data for new matrix
    138         * @param copyArray if true, the input array will be copied, otherwise
    139         * it will be referenced
    140         *
    141         * @exception IllegalArgumentException if <code>blockData</code> shape is
    142         * inconsistent with block layout
    143         * @see #createBlocksLayout(int, int)
    144         * @see #toBlocksLayout(double[][])
    145         * @see #BlockRealMatrix(double[][])
    146         */
    147        public BlockRealMatrix(final int rows, final int columns,
    148                               final double[][] blockData, final boolean copyArray)
    149            throws IllegalArgumentException {
    150
    151            super(rows, columns);
    152            this.rows    = rows;
    153            this.columns = columns;
    154
    155            // number of blocks
    156            blockRows    = (rows    + BLOCK_SIZE - 1) / BLOCK_SIZE;
    157            blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
    158
    159            if (copyArray) {
    160                // allocate storage blocks, taking care of smaller ones at right and bottom
    161                blocks = new double[blockRows * blockColumns][];
    162            } else {
    163                // reference existing array
    164                blocks = blockData;
    165            }
    166
    167            int index = 0;
    168            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    169                final int iHeight = blockHeight(iBlock);
    170                for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++index) {
    171                    if (blockData[index].length != iHeight * blockWidth(jBlock)) {
    172                        throw MathRuntimeException.createIllegalArgumentException(
    173                                LocalizedFormats.WRONG_BLOCK_LENGTH,
    174                                blockData[index].length, iHeight * blockWidth(jBlock));
    175                    }
    176                    if (copyArray) {
    177                        blocks[index] = blockData[index].clone();
    178                    }
    179                }
    180            }
    181
    182        }
    183
    184        /**
    185         * Convert a data array from raw layout to blocks layout.
    186         * <p>
    187         * Raw layout is the straightforward layout where element at row i and
    188         * column j is in array element <code>rawData[i][j]</code>. Blocks layout
    189         * is the layout used in {@link BlockRealMatrix} instances, where the matrix
    190         * is split in square blocks (except at right and bottom side where blocks may
    191         * be rectangular to fit matrix size) and each block is stored in a flattened
    192         * one-dimensional array.
    193         * </p>
    194         * <p>
    195         * This method creates an array in blocks layout from an input array in raw layout.
    196         * It can be used to provide the array argument of the {@link
    197         * #BlockRealMatrix(int, int, double[][], boolean)} constructor.
    198         * </p>
    199         * @param rawData data array in raw layout
    200         * @return a new data array containing the same entries but in blocks layout
    201         * @exception IllegalArgumentException if <code>rawData</code> is not rectangular
    202         *  (not all rows have the same length)
    203         * @see #createBlocksLayout(int, int)
    204         * @see #BlockRealMatrix(int, int, double[][], boolean)
    205         */
    206        public static double[][] toBlocksLayout(final double[][] rawData)
    207            throws IllegalArgumentException {
    208
    209            final int rows         = rawData.length;
    210            final int columns      = rawData[0].length;
    211            final int blockRows    = (rows    + BLOCK_SIZE - 1) / BLOCK_SIZE;
    212            final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
    213
    214            // safety checks
    215            for (int i = 0; i < rawData.length; ++i) {
    216                final int length = rawData[i].length;
    217                if (length != columns) {
    218                    throw MathRuntimeException.createIllegalArgumentException(
    219                            LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    220                            columns, length);
    221                }
    222            }
    223
    224            // convert array
    225            final double[][] blocks = new double[blockRows * blockColumns][];
    226            int blockIndex = 0;
    227            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    228                final int pStart  = iBlock * BLOCK_SIZE;
    229                final int pEnd    = FastMath.min(pStart + BLOCK_SIZE, rows);
    230                final int iHeight = pEnd - pStart;
    231                for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    232                    final int qStart = jBlock * BLOCK_SIZE;
    233                    final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    234                    final int jWidth = qEnd - qStart;
    235
    236                    // allocate new block
    237                    final double[] block = new double[iHeight * jWidth];
    238                    blocks[blockIndex] = block;
    239
    240                    // copy data
    241                    int index = 0;
    242                    for (int p = pStart; p < pEnd; ++p) {
    243                        System.arraycopy(rawData[p], qStart, block, index, jWidth);
    244                        index += jWidth;
    245                    }
    246
    247                    ++blockIndex;
    248
    249                }
    250            }
    251
    252            return blocks;
    253
    254        }
    255
    256        /**
    257         * Create a data array in blocks layout.
    258         * <p>
    259         * This method can be used to create the array argument of the {@link
    260         * #BlockRealMatrix(int, int, double[][], boolean)} constructor.
    261         * </p>
    262         * @param rows  the number of rows in the new matrix
    263         * @param columns  the number of columns in the new matrix
    264         * @return a new data array in blocks layout
    265         * @see #toBlocksLayout(double[][])
    266         * @see #BlockRealMatrix(int, int, double[][], boolean)
    267         */
    268        public static double[][] createBlocksLayout(final int rows, final int columns) {
    269
    270            final int blockRows    = (rows    + BLOCK_SIZE - 1) / BLOCK_SIZE;
    271            final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
    272
    273            final double[][] blocks = new double[blockRows * blockColumns][];
    274            int blockIndex = 0;
    275            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    276                final int pStart  = iBlock * BLOCK_SIZE;
    277                final int pEnd    = FastMath.min(pStart + BLOCK_SIZE, rows);
    278                final int iHeight = pEnd - pStart;
    279                for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    280                    final int qStart = jBlock * BLOCK_SIZE;
    281                    final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    282                    final int jWidth = qEnd - qStart;
    283                    blocks[blockIndex] = new double[iHeight * jWidth];
    284                    ++blockIndex;
    285                }
    286            }
    287
    288            return blocks;
    289
    290        }
    291
    292        /** {@inheritDoc} */
    293        @Override
    294        public BlockRealMatrix createMatrix(final int rowDimension, final int columnDimension)
    295            throws IllegalArgumentException {
    296            return new BlockRealMatrix(rowDimension, columnDimension);
    297        }
    298
    299        /** {@inheritDoc} */
    300        @Override
    301        public BlockRealMatrix copy() {
    302
    303            // create an empty matrix
    304            BlockRealMatrix copied = new BlockRealMatrix(rows, columns);
    305
    306            // copy the blocks
    307            for (int i = 0; i < blocks.length; ++i) {
    308                System.arraycopy(blocks[i], 0, copied.blocks[i], 0, blocks[i].length);
    309            }
    310
    311            return copied;
    312
    313        }
    314
    315        /** {@inheritDoc} */
    316        @Override
    317        public BlockRealMatrix add(final RealMatrix m)
    318            throws IllegalArgumentException {
    319            try {
    320                return add((BlockRealMatrix) m);
    321            } catch (ClassCastException cce) {
    322
    323                // safety check
    324                MatrixUtils.checkAdditionCompatible(this, m);
    325
    326                final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    327
    328                // perform addition block-wise, to ensure good cache behavior
    329                int blockIndex = 0;
    330                for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
    331                    for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
    332
    333                        // perform addition on the current block
    334                        final double[] outBlock = out.blocks[blockIndex];
    335                        final double[] tBlock   = blocks[blockIndex];
    336                        final int      pStart   = iBlock * BLOCK_SIZE;
    337                        final int      pEnd     = FastMath.min(pStart + BLOCK_SIZE, rows);
    338                        final int      qStart   = jBlock * BLOCK_SIZE;
    339                        final int      qEnd     = FastMath.min(qStart + BLOCK_SIZE, columns);
    340                        int k = 0;
    341                        for (int p = pStart; p < pEnd; ++p) {
    342                            for (int q = qStart; q < qEnd; ++q) {
    343                                outBlock[k] = tBlock[k] + m.getEntry(p, q);
    344                                ++k;
    345                            }
    346                        }
    347
    348                        // go to next block
    349                        ++blockIndex;
    350
    351                    }
    352                }
    353
    354                return out;
    355
    356            }
    357        }
    358
    359        /**
    360         * Compute the sum of this and <code>m</code>.
    361         *
    362         * @param m    matrix to be added
    363         * @return     this + m
    364         * @throws  IllegalArgumentException if m is not the same size as this
    365         */
    366        public BlockRealMatrix add(final BlockRealMatrix m)
    367            throws IllegalArgumentException {
    368
    369            // safety check
    370            MatrixUtils.checkAdditionCompatible(this, m);
    371
    372            final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    373
    374            // perform addition block-wise, to ensure good cache behavior
    375            for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
    376                final double[] outBlock = out.blocks[blockIndex];
    377                final double[] tBlock   = blocks[blockIndex];
    378                final double[] mBlock   = m.blocks[blockIndex];
    379                for (int k = 0; k < outBlock.length; ++k) {
    380                    outBlock[k] = tBlock[k] + mBlock[k];
    381                }
    382            }
    383
    384            return out;
    385
    386        }
    387
    388        /** {@inheritDoc} */
    389        @Override
    390        public BlockRealMatrix subtract(final RealMatrix m)
    391            throws IllegalArgumentException {
    392            try {
    393                return subtract((BlockRealMatrix) m);
    394            } catch (ClassCastException cce) {
    395
    396                // safety check
    397                MatrixUtils.checkSubtractionCompatible(this, m);
    398
    399                final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    400
    401                // perform subtraction block-wise, to ensure good cache behavior
    402                int blockIndex = 0;
    403                for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
    404                    for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
    405
    406                        // perform subtraction on the current block
    407                        final double[] outBlock = out.blocks[blockIndex];
    408                        final double[] tBlock   = blocks[blockIndex];
    409                        final int      pStart   = iBlock * BLOCK_SIZE;
    410                        final int      pEnd     = FastMath.min(pStart + BLOCK_SIZE, rows);
    411                        final int      qStart   = jBlock * BLOCK_SIZE;
    412                        final int      qEnd     = FastMath.min(qStart + BLOCK_SIZE, columns);
    413                        int k = 0;
    414                        for (int p = pStart; p < pEnd; ++p) {
    415                            for (int q = qStart; q < qEnd; ++q) {
    416                                outBlock[k] = tBlock[k] - m.getEntry(p, q);
    417                                ++k;
    418                            }
    419                        }
    420
    421                        // go to next block
    422                        ++blockIndex;
    423
    424                    }
    425                }
    426
    427                return out;
    428
    429            }
    430        }
    431
    432        /**
    433         * Compute this minus <code>m</code>.
    434         *
    435         * @param m    matrix to be subtracted
    436         * @return     this - m
    437         * @throws  IllegalArgumentException if m is not the same size as this
    438         */
    439        public BlockRealMatrix subtract(final BlockRealMatrix m)
    440            throws IllegalArgumentException {
    441
    442            // safety check
    443            MatrixUtils.checkSubtractionCompatible(this, m);
    444
    445            final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    446
    447            // perform subtraction block-wise, to ensure good cache behavior
    448            for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
    449                final double[] outBlock = out.blocks[blockIndex];
    450                final double[] tBlock   = blocks[blockIndex];
    451                final double[] mBlock   = m.blocks[blockIndex];
    452                for (int k = 0; k < outBlock.length; ++k) {
    453                    outBlock[k] = tBlock[k] - mBlock[k];
    454                }
    455            }
    456
    457            return out;
    458
    459        }
    460
    461        /** {@inheritDoc} */
    462        @Override
    463        public BlockRealMatrix scalarAdd(final double d)
    464            throws IllegalArgumentException {
    465
    466            final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    467
    468            // perform subtraction block-wise, to ensure good cache behavior
    469            for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
    470                final double[] outBlock = out.blocks[blockIndex];
    471                final double[] tBlock   = blocks[blockIndex];
    472                for (int k = 0; k < outBlock.length; ++k) {
    473                    outBlock[k] = tBlock[k] + d;
    474                }
    475            }
    476
    477            return out;
    478
    479        }
    480
    481        /** {@inheritDoc} */
    482        @Override
    483        public RealMatrix scalarMultiply(final double d)
    484            throws IllegalArgumentException {
    485
    486            final BlockRealMatrix out = new BlockRealMatrix(rows, columns);
    487
    488            // perform subtraction block-wise, to ensure good cache behavior
    489            for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
    490                final double[] outBlock = out.blocks[blockIndex];
    491                final double[] tBlock   = blocks[blockIndex];
    492                for (int k = 0; k < outBlock.length; ++k) {
    493                    outBlock[k] = tBlock[k] * d;
    494                }
    495            }
    496
    497            return out;
    498
    499        }
    500
    501        /** {@inheritDoc} */
    502        @Override
    503        public BlockRealMatrix multiply(final RealMatrix m)
    504            throws IllegalArgumentException {
    505            try {
    506                return multiply((BlockRealMatrix) m);
    507            } catch (ClassCastException cce) {
    508
    509                // safety check
    510                MatrixUtils.checkMultiplicationCompatible(this, m);
    511
    512                final BlockRealMatrix out = new BlockRealMatrix(rows, m.getColumnDimension());
    513
    514                // perform multiplication block-wise, to ensure good cache behavior
    515                int blockIndex = 0;
    516                for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
    517
    518                    final int pStart = iBlock * BLOCK_SIZE;
    519                    final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    520
    521                    for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
    522
    523                        final int qStart = jBlock * BLOCK_SIZE;
    524                        final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, m.getColumnDimension());
    525
    526                        // select current block
    527                        final double[] outBlock = out.blocks[blockIndex];
    528
    529                        // perform multiplication on current block
    530                        for (int kBlock = 0; kBlock < blockColumns; ++kBlock) {
    531                            final int kWidth      = blockWidth(kBlock);
    532                            final double[] tBlock = blocks[iBlock * blockColumns + kBlock];
    533                            final int rStart      = kBlock * BLOCK_SIZE;
    534                            int k = 0;
    535                            for (int p = pStart; p < pEnd; ++p) {
    536                                final int lStart = (p - pStart) * kWidth;
    537                                final int lEnd   = lStart + kWidth;
    538                                for (int q = qStart; q < qEnd; ++q) {
    539                                    double sum = 0;
    540                                    int r = rStart;
    541                                    for (int l = lStart; l < lEnd; ++l) {
    542                                        sum += tBlock[l] * m.getEntry(r, q);
    543                                        ++r;
    544                                    }
    545                                    outBlock[k] += sum;
    546                                    ++k;
    547                                }
    548                            }
    549                        }
    550
    551                        // go to next block
    552                        ++blockIndex;
    553
    554                    }
    555                }
    556
    557                return out;
    558
    559            }
    560        }
    561
    562        /**
    563         * Returns the result of postmultiplying this by m.
    564         *
    565         * @param m    matrix to postmultiply by
    566         * @return     this * m
    567         * @throws     IllegalArgumentException
    568         *             if columnDimension(this) != rowDimension(m)
    569         */
    570        public BlockRealMatrix multiply(BlockRealMatrix m) throws IllegalArgumentException {
    571
    572            // safety check
    573            MatrixUtils.checkMultiplicationCompatible(this, m);
    574
    575            final BlockRealMatrix out = new BlockRealMatrix(rows, m.columns);
    576
    577            // perform multiplication block-wise, to ensure good cache behavior
    578            int blockIndex = 0;
    579            for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
    580
    581                final int pStart = iBlock * BLOCK_SIZE;
    582                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    583
    584                for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
    585                    final int jWidth = out.blockWidth(jBlock);
    586                    final int jWidth2 = jWidth  + jWidth;
    587                    final int jWidth3 = jWidth2 + jWidth;
    588                    final int jWidth4 = jWidth3 + jWidth;
    589
    590                    // select current block
    591                    final double[] outBlock = out.blocks[blockIndex];
    592
    593                    // perform multiplication on current block
    594                    for (int kBlock = 0; kBlock < blockColumns; ++kBlock) {
    595                        final int kWidth = blockWidth(kBlock);
    596                        final double[] tBlock = blocks[iBlock * blockColumns + kBlock];
    597                        final double[] mBlock = m.blocks[kBlock * m.blockColumns + jBlock];
    598                        int k = 0;
    599                        for (int p = pStart; p < pEnd; ++p) {
    600                            final int lStart = (p - pStart) * kWidth;
    601                            final int lEnd   = lStart + kWidth;
    602                            for (int nStart = 0; nStart < jWidth; ++nStart) {
    603                                double sum = 0;
    604                                int l = lStart;
    605                                int n = nStart;
    606                                while (l < lEnd - 3) {
    607                                    sum += tBlock[l] * mBlock[n] +
    608                                           tBlock[l + 1] * mBlock[n + jWidth] +
    609                                           tBlock[l + 2] * mBlock[n + jWidth2] +
    610                                           tBlock[l + 3] * mBlock[n + jWidth3];
    611                                    l += 4;
    612                                    n += jWidth4;
    613                                }
    614                                while (l < lEnd) {
    615                                    sum += tBlock[l++] * mBlock[n];
    616                                    n += jWidth;
    617                                }
    618                                outBlock[k] += sum;
    619                                ++k;
    620                            }
    621                        }
    622                    }
    623
    624                    // go to next block
    625                    ++blockIndex;
    626
    627                }
    628            }
    629
    630            return out;
    631
    632        }
    633
    634        /** {@inheritDoc} */
    635        @Override
    636        public double[][] getData() {
    637
    638            final double[][] data = new double[getRowDimension()][getColumnDimension()];
    639            final int lastColumns = columns - (blockColumns - 1) * BLOCK_SIZE;
    640
    641            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    642                final int pStart = iBlock * BLOCK_SIZE;
    643                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    644                int regularPos   = 0;
    645                int lastPos      = 0;
    646                for (int p = pStart; p < pEnd; ++p) {
    647                    final double[] dataP = data[p];
    648                    int blockIndex = iBlock * blockColumns;
    649                    int dataPos    = 0;
    650                    for (int jBlock = 0; jBlock < blockColumns - 1; ++jBlock) {
    651                        System.arraycopy(blocks[blockIndex++], regularPos, dataP, dataPos, BLOCK_SIZE);
    652                        dataPos += BLOCK_SIZE;
    653                    }
    654                    System.arraycopy(blocks[blockIndex], lastPos, dataP, dataPos, lastColumns);
    655                    regularPos += BLOCK_SIZE;
    656                    lastPos    += lastColumns;
    657                }
    658            }
    659
    660            return data;
    661
    662        }
    663
    664        /** {@inheritDoc} */
    665        @Override
    666        public double getNorm() {
    667            final double[] colSums = new double[BLOCK_SIZE];
    668            double maxColSum = 0;
    669            for (int jBlock = 0; jBlock < blockColumns; jBlock++) {
    670                final int jWidth = blockWidth(jBlock);
    671                Arrays.fill(colSums, 0, jWidth, 0.0);
    672                for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    673                    final int iHeight = blockHeight(iBlock);
    674                    final double[] block = blocks[iBlock * blockColumns + jBlock];
    675                    for (int j = 0; j < jWidth; ++j) {
    676                        double sum = 0;
    677                        for (int i = 0; i < iHeight; ++i) {
    678                            sum += FastMath.abs(block[i * jWidth + j]);
    679                        }
    680                        colSums[j] += sum;
    681                    }
    682                }
    683                for (int j = 0; j < jWidth; ++j) {
    684                    maxColSum = FastMath.max(maxColSum, colSums[j]);
    685                }
    686            }
    687            return maxColSum;
    688        }
    689
    690        /** {@inheritDoc} */
    691        @Override
    692        public double getFrobeniusNorm() {
    693            double sum2 = 0;
    694            for (int blockIndex = 0; blockIndex < blocks.length; ++blockIndex) {
    695                for (final double entry : blocks[blockIndex]) {
    696                    sum2 += entry * entry;
    697                }
    698            }
    699            return FastMath.sqrt(sum2);
    700        }
    701
    702        /** {@inheritDoc} */
    703        @Override
    704        public BlockRealMatrix getSubMatrix(final int startRow, final int endRow,
    705                                       final int startColumn, final int endColumn)
    706            throws MatrixIndexException {
    707
    708            // safety checks
    709            MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    710
    711            // create the output matrix
    712            final BlockRealMatrix out =
    713                new BlockRealMatrix(endRow - startRow + 1, endColumn - startColumn + 1);
    714
    715            // compute blocks shifts
    716            final int blockStartRow    = startRow    / BLOCK_SIZE;
    717            final int rowsShift        = startRow    % BLOCK_SIZE;
    718            final int blockStartColumn = startColumn / BLOCK_SIZE;
    719            final int columnsShift     = startColumn % BLOCK_SIZE;
    720
    721            // perform extraction block-wise, to ensure good cache behavior
    722            int pBlock = blockStartRow;
    723            for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
    724                final int iHeight = out.blockHeight(iBlock);
    725                int qBlock = blockStartColumn;
    726                for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
    727                    final int jWidth = out.blockWidth(jBlock);
    728
    729                    // handle one block of the output matrix
    730                    final int      outIndex = iBlock * out.blockColumns + jBlock;
    731                    final double[] outBlock = out.blocks[outIndex];
    732                    final int      index    = pBlock * blockColumns + qBlock;
    733                    final int      width    = blockWidth(qBlock);
    734
    735                    final int heightExcess = iHeight + rowsShift - BLOCK_SIZE;
    736                    final int widthExcess  = jWidth + columnsShift - BLOCK_SIZE;
    737                    if (heightExcess > 0) {
    738                        // the submatrix block spans on two blocks rows from the original matrix
    739                        if (widthExcess > 0) {
    740                            // the submatrix block spans on two blocks columns from the original matrix
    741                            final int width2 = blockWidth(qBlock + 1);
    742                            copyBlockPart(blocks[index], width,
    743                                          rowsShift, BLOCK_SIZE,
    744                                          columnsShift, BLOCK_SIZE,
    745                                          outBlock, jWidth, 0, 0);
    746                            copyBlockPart(blocks[index + 1], width2,
    747                                          rowsShift, BLOCK_SIZE,
    748                                          0, widthExcess,
    749                                          outBlock, jWidth, 0, jWidth - widthExcess);
    750                            copyBlockPart(blocks[index + blockColumns], width,
    751                                          0, heightExcess,
    752                                          columnsShift, BLOCK_SIZE,
    753                                          outBlock, jWidth, iHeight - heightExcess, 0);
    754                            copyBlockPart(blocks[index + blockColumns + 1], width2,
    755                                          0, heightExcess,
    756                                          0, widthExcess,
    757                                          outBlock, jWidth, iHeight - heightExcess, jWidth - widthExcess);
    758                        } else {
    759                            // the submatrix block spans on one block column from the original matrix
    760                            copyBlockPart(blocks[index], width,
    761                                          rowsShift, BLOCK_SIZE,
    762                                          columnsShift, jWidth + columnsShift,
    763                                          outBlock, jWidth, 0, 0);
    764                            copyBlockPart(blocks[index + blockColumns], width,
    765                                          0, heightExcess,
    766                                          columnsShift, jWidth + columnsShift,
    767                                          outBlock, jWidth, iHeight - heightExcess, 0);
    768                        }
    769                    } else {
    770                        // the submatrix block spans on one block row from the original matrix
    771                        if (widthExcess > 0) {
    772                            // the submatrix block spans on two blocks columns from the original matrix
    773                            final int width2 = blockWidth(qBlock + 1);
    774                            copyBlockPart(blocks[index], width,
    775                                          rowsShift, iHeight + rowsShift,
    776                                          columnsShift, BLOCK_SIZE,
    777                                          outBlock, jWidth, 0, 0);
    778                            copyBlockPart(blocks[index + 1], width2,
    779                                          rowsShift, iHeight + rowsShift,
    780                                          0, widthExcess,
    781                                          outBlock, jWidth, 0, jWidth - widthExcess);
    782                        } else {
    783                            // the submatrix block spans on one block column from the original matrix
    784                            copyBlockPart(blocks[index], width,
    785                                          rowsShift, iHeight + rowsShift,
    786                                          columnsShift, jWidth + columnsShift,
    787                                          outBlock, jWidth, 0, 0);
    788                        }
    789                   }
    790
    791                    ++qBlock;
    792
    793                }
    794
    795                ++pBlock;
    796
    797            }
    798
    799            return out;
    800
    801        }
    802
    803        /**
    804         * Copy a part of a block into another one
    805         * <p>This method can be called only when the specified part fits in both
    806         * blocks, no verification is done here.</p>
    807         * @param srcBlock source block
    808         * @param srcWidth source block width ({@link #BLOCK_SIZE} or smaller)
    809         * @param srcStartRow start row in the source block
    810         * @param srcEndRow end row (exclusive) in the source block
    811         * @param srcStartColumn start column in the source block
    812         * @param srcEndColumn end column (exclusive) in the source block
    813         * @param dstBlock destination block
    814         * @param dstWidth destination block width ({@link #BLOCK_SIZE} or smaller)
    815         * @param dstStartRow start row in the destination block
    816         * @param dstStartColumn start column in the destination block
    817         */
    818        private void copyBlockPart(final double[] srcBlock, final int srcWidth,
    819                                   final int srcStartRow, final int srcEndRow,
    820                                   final int srcStartColumn, final int srcEndColumn,
    821                                   final double[] dstBlock, final int dstWidth,
    822                                   final int dstStartRow, final int dstStartColumn) {
    823            final int length = srcEndColumn - srcStartColumn;
    824            int srcPos = srcStartRow * srcWidth + srcStartColumn;
    825            int dstPos = dstStartRow * dstWidth + dstStartColumn;
    826            for (int srcRow = srcStartRow; srcRow < srcEndRow; ++srcRow) {
    827                System.arraycopy(srcBlock, srcPos, dstBlock, dstPos, length);
    828                srcPos += srcWidth;
    829                dstPos += dstWidth;
    830            }
    831        }
    832
    833        /** {@inheritDoc} */
    834        @Override
    835        public void setSubMatrix(final double[][] subMatrix, final int row, final int column)
    836            throws MatrixIndexException {
    837
    838            // safety checks
    839            final int refLength = subMatrix[0].length;
    840            if (refLength < 1) {
    841                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
    842            }
    843            final int endRow    = row + subMatrix.length - 1;
    844            final int endColumn = column + refLength - 1;
    845            MatrixUtils.checkSubMatrixIndex(this, row, endRow, column, endColumn);
    846            for (final double[] subRow : subMatrix) {
    847                if (subRow.length != refLength) {
    848                    throw MathRuntimeException.createIllegalArgumentException(
    849                            LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    850                            refLength, subRow.length);
    851                }
    852            }
    853
    854            // compute blocks bounds
    855            final int blockStartRow    = row / BLOCK_SIZE;
    856            final int blockEndRow      = (endRow + BLOCK_SIZE) / BLOCK_SIZE;
    857            final int blockStartColumn = column / BLOCK_SIZE;
    858            final int blockEndColumn   = (endColumn + BLOCK_SIZE) / BLOCK_SIZE;
    859
    860            // perform copy block-wise, to ensure good cache behavior
    861            for (int iBlock = blockStartRow; iBlock < blockEndRow; ++iBlock) {
    862                final int iHeight  = blockHeight(iBlock);
    863                final int firstRow = iBlock * BLOCK_SIZE;
    864                final int iStart   = FastMath.max(row,    firstRow);
    865                final int iEnd     = FastMath.min(endRow + 1, firstRow + iHeight);
    866
    867                for (int jBlock = blockStartColumn; jBlock < blockEndColumn; ++jBlock) {
    868                    final int jWidth      = blockWidth(jBlock);
    869                    final int firstColumn = jBlock * BLOCK_SIZE;
    870                    final int jStart      = FastMath.max(column,    firstColumn);
    871                    final int jEnd        = FastMath.min(endColumn + 1, firstColumn + jWidth);
    872                    final int jLength     = jEnd - jStart;
    873
    874                    // handle one block, row by row
    875                    final double[] block = blocks[iBlock * blockColumns + jBlock];
    876                    for (int i = iStart; i < iEnd; ++i) {
    877                        System.arraycopy(subMatrix[i - row], jStart - column,
    878                                         block, (i - firstRow) * jWidth + (jStart - firstColumn),
    879                                         jLength);
    880                    }
    881
    882                }
    883            }
    884        }
    885
    886        /** {@inheritDoc} */
    887        @Override
    888        public BlockRealMatrix getRowMatrix(final int row)
    889            throws MatrixIndexException {
    890
    891            MatrixUtils.checkRowIndex(this, row);
    892            final BlockRealMatrix out = new BlockRealMatrix(1, columns);
    893
    894            // perform copy block-wise, to ensure good cache behavior
    895            final int iBlock  = row / BLOCK_SIZE;
    896            final int iRow    = row - iBlock * BLOCK_SIZE;
    897            int outBlockIndex = 0;
    898            int outIndex      = 0;
    899            double[] outBlock = out.blocks[outBlockIndex];
    900            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    901                final int jWidth     = blockWidth(jBlock);
    902                final double[] block = blocks[iBlock * blockColumns + jBlock];
    903                final int available  = outBlock.length - outIndex;
    904                if (jWidth > available) {
    905                    System.arraycopy(block, iRow * jWidth, outBlock, outIndex, available);
    906                    outBlock = out.blocks[++outBlockIndex];
    907                    System.arraycopy(block, iRow * jWidth, outBlock, 0, jWidth - available);
    908                    outIndex = jWidth - available;
    909                } else {
    910                    System.arraycopy(block, iRow * jWidth, outBlock, outIndex, jWidth);
    911                    outIndex += jWidth;
    912                }
    913            }
    914
    915            return out;
    916
    917        }
    918
    919        /** {@inheritDoc} */
    920        @Override
    921        public void setRowMatrix(final int row, final RealMatrix matrix)
    922            throws MatrixIndexException, InvalidMatrixException {
    923            try {
    924                setRowMatrix(row, (BlockRealMatrix) matrix);
    925            } catch (ClassCastException cce) {
    926                super.setRowMatrix(row, matrix);
    927            }
    928        }
    929
    930        /**
    931         * Sets the entries in row number <code>row</code>
    932         * as a row matrix.  Row indices start at 0.
    933         *
    934         * @param row the row to be set
    935         * @param matrix row matrix (must have one row and the same number of columns
    936         * as the instance)
    937         * @throws MatrixIndexException if the specified row index is invalid
    938         * @throws InvalidMatrixException if the matrix dimensions do not match one
    939         * instance row
    940         */
    941        public void setRowMatrix(final int row, final BlockRealMatrix matrix)
    942            throws MatrixIndexException, InvalidMatrixException {
    943
    944            MatrixUtils.checkRowIndex(this, row);
    945            final int nCols = getColumnDimension();
    946            if ((matrix.getRowDimension() != 1) ||
    947                (matrix.getColumnDimension() != nCols)) {
    948                throw new InvalidMatrixException(
    949                        LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    950                        matrix.getRowDimension(), matrix.getColumnDimension(),
    951                        1, nCols);
    952            }
    953
    954            // perform copy block-wise, to ensure good cache behavior
    955            final int iBlock = row / BLOCK_SIZE;
    956            final int iRow   = row - iBlock * BLOCK_SIZE;
    957            int mBlockIndex  = 0;
    958            int mIndex       = 0;
    959            double[] mBlock  = matrix.blocks[mBlockIndex];
    960            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    961                final int jWidth     = blockWidth(jBlock);
    962                final double[] block = blocks[iBlock * blockColumns + jBlock];
    963                final int available  = mBlock.length - mIndex;
    964                if (jWidth > available) {
    965                    System.arraycopy(mBlock, mIndex, block, iRow * jWidth, available);
    966                    mBlock = matrix.blocks[++mBlockIndex];
    967                    System.arraycopy(mBlock, 0, block, iRow * jWidth, jWidth - available);
    968                    mIndex = jWidth - available;
    969                } else {
    970                    System.arraycopy(mBlock, mIndex, block, iRow * jWidth, jWidth);
    971                    mIndex += jWidth;
    972               }
    973            }
    974
    975        }
    976
    977        /** {@inheritDoc} */
    978        @Override
    979        public BlockRealMatrix getColumnMatrix(final int column)
    980            throws MatrixIndexException {
    981
    982            MatrixUtils.checkColumnIndex(this, column);
    983            final BlockRealMatrix out = new BlockRealMatrix(rows, 1);
    984
    985            // perform copy block-wise, to ensure good cache behavior
    986            final int jBlock  = column / BLOCK_SIZE;
    987            final int jColumn = column - jBlock * BLOCK_SIZE;
    988            final int jWidth  = blockWidth(jBlock);
    989            int outBlockIndex = 0;
    990            int outIndex      = 0;
    991            double[] outBlock = out.blocks[outBlockIndex];
    992            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    993                final int iHeight = blockHeight(iBlock);
    994                final double[] block = blocks[iBlock * blockColumns + jBlock];
    995                for (int i = 0; i < iHeight; ++i) {
    996                    if (outIndex >= outBlock.length) {
    997                        outBlock = out.blocks[++outBlockIndex];
    998                        outIndex = 0;
    999                    }
    1000                    outBlock[outIndex++] = block[i * jWidth + jColumn];
    1001                }
    1002            }
    1003
    1004            return out;
    1005
    1006        }
    1007
    1008        /** {@inheritDoc} */
    1009        @Override
    1010        public void setColumnMatrix(final int column, final RealMatrix matrix)
    1011            throws MatrixIndexException, InvalidMatrixException {
    1012            try {
    1013                setColumnMatrix(column, (BlockRealMatrix) matrix);
    1014            } catch (ClassCastException cce) {
    1015                super.setColumnMatrix(column, matrix);
    1016            }
    1017        }
    1018
    1019        /**
    1020         * Sets the entries in column number <code>column</code>
    1021         * as a column matrix.  Column indices start at 0.
    1022         *
    1023         * @param column the column to be set
    1024         * @param matrix column matrix (must have one column and the same number of rows
    1025         * as the instance)
    1026         * @throws MatrixIndexException if the specified column index is invalid
    1027         * @throws InvalidMatrixException if the matrix dimensions do not match one
    1028         * instance column
    1029         */
    1030        void setColumnMatrix(final int column, final BlockRealMatrix matrix)
    1031            throws MatrixIndexException, InvalidMatrixException {
    1032
    1033            MatrixUtils.checkColumnIndex(this, column);
    1034            final int nRows = getRowDimension();
    1035            if ((matrix.getRowDimension() != nRows) ||
    1036                (matrix.getColumnDimension() != 1)) {
    1037                throw new InvalidMatrixException(
    1038                        LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    1039                        matrix.getRowDimension(), matrix.getColumnDimension(),
    1040                        nRows, 1);
    1041            }
    1042
    1043            // perform copy block-wise, to ensure good cache behavior
    1044            final int jBlock  = column / BLOCK_SIZE;
    1045            final int jColumn = column - jBlock * BLOCK_SIZE;
    1046            final int jWidth  = blockWidth(jBlock);
    1047            int mBlockIndex = 0;
    1048            int mIndex      = 0;
    1049            double[] mBlock = matrix.blocks[mBlockIndex];
    1050            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1051                final int iHeight = blockHeight(iBlock);
    1052                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1053                for (int i = 0; i < iHeight; ++i) {
    1054                    if (mIndex >= mBlock.length) {
    1055                        mBlock = matrix.blocks[++mBlockIndex];
    1056                        mIndex = 0;
    1057                    }
    1058                    block[i * jWidth + jColumn] = mBlock[mIndex++];
    1059                }
    1060            }
    1061
    1062        }
    1063
    1064        /** {@inheritDoc} */
    1065        @Override
    1066        public RealVector getRowVector(final int row)
    1067            throws MatrixIndexException {
    1068
    1069            MatrixUtils.checkRowIndex(this, row);
    1070            final double[] outData = new double[columns];
    1071
    1072            // perform copy block-wise, to ensure good cache behavior
    1073            final int iBlock  = row / BLOCK_SIZE;
    1074            final int iRow    = row - iBlock * BLOCK_SIZE;
    1075            int outIndex      = 0;
    1076            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1077                final int jWidth     = blockWidth(jBlock);
    1078                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1079                System.arraycopy(block, iRow * jWidth, outData, outIndex, jWidth);
    1080                outIndex += jWidth;
    1081            }
    1082
    1083            return new ArrayRealVector(outData, false);
    1084
    1085        }
    1086
    1087        /** {@inheritDoc} */
    1088        @Override
    1089        public void setRowVector(final int row, final RealVector vector)
    1090            throws MatrixIndexException, InvalidMatrixException {
    1091            try {
    1092                setRow(row, ((ArrayRealVector) vector).getDataRef());
    1093            } catch (ClassCastException cce) {
    1094                super.setRowVector(row, vector);
    1095            }
    1096        }
    1097
    1098        /** {@inheritDoc} */
    1099        @Override
    1100        public RealVector getColumnVector(final int column)
    1101            throws MatrixIndexException {
    1102
    1103            MatrixUtils.checkColumnIndex(this, column);
    1104            final double[] outData = new double[rows];
    1105
    1106            // perform copy block-wise, to ensure good cache behavior
    1107            final int jBlock  = column / BLOCK_SIZE;
    1108            final int jColumn = column - jBlock * BLOCK_SIZE;
    1109            final int jWidth  = blockWidth(jBlock);
    1110            int outIndex      = 0;
    1111            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1112                final int iHeight = blockHeight(iBlock);
    1113                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1114                for (int i = 0; i < iHeight; ++i) {
    1115                    outData[outIndex++] = block[i * jWidth + jColumn];
    1116                }
    1117            }
    1118
    1119            return new ArrayRealVector(outData, false);
    1120
    1121        }
    1122
    1123        /** {@inheritDoc} */
    1124        @Override
    1125        public void setColumnVector(final int column, final RealVector vector)
    1126            throws MatrixIndexException, InvalidMatrixException {
    1127            try {
    1128                setColumn(column, ((ArrayRealVector) vector).getDataRef());
    1129            } catch (ClassCastException cce) {
    1130                super.setColumnVector(column, vector);
    1131            }
    1132        }
    1133
    1134        /** {@inheritDoc} */
    1135        @Override
    1136        public double[] getRow(final int row)
    1137            throws MatrixIndexException {
    1138
    1139            MatrixUtils.checkRowIndex(this, row);
    1140            final double[] out = new double[columns];
    1141
    1142            // perform copy block-wise, to ensure good cache behavior
    1143            final int iBlock  = row / BLOCK_SIZE;
    1144            final int iRow    = row - iBlock * BLOCK_SIZE;
    1145            int outIndex      = 0;
    1146            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1147                final int jWidth     = blockWidth(jBlock);
    1148                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1149                System.arraycopy(block, iRow * jWidth, out, outIndex, jWidth);
    1150                outIndex += jWidth;
    1151            }
    1152
    1153            return out;
    1154
    1155        }
    1156
    1157        /** {@inheritDoc} */
    1158        @Override
    1159        public void setRow(final int row, final double[] array)
    1160            throws MatrixIndexException, InvalidMatrixException {
    1161
    1162            MatrixUtils.checkRowIndex(this, row);
    1163            final int nCols = getColumnDimension();
    1164            if (array.length != nCols) {
    1165                throw new InvalidMatrixException(
    1166                        LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    1167                        1, array.length, 1, nCols);
    1168            }
    1169
    1170            // perform copy block-wise, to ensure good cache behavior
    1171            final int iBlock  = row / BLOCK_SIZE;
    1172            final int iRow    = row - iBlock * BLOCK_SIZE;
    1173            int outIndex      = 0;
    1174            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1175                final int jWidth     = blockWidth(jBlock);
    1176                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1177                System.arraycopy(array, outIndex, block, iRow * jWidth, jWidth);
    1178                outIndex += jWidth;
    1179            }
    1180
    1181        }
    1182
    1183        /** {@inheritDoc} */
    1184        @Override
    1185        public double[] getColumn(final int column)
    1186            throws MatrixIndexException {
    1187
    1188            MatrixUtils.checkColumnIndex(this, column);
    1189            final double[] out = new double[rows];
    1190
    1191            // perform copy block-wise, to ensure good cache behavior
    1192            final int jBlock  = column / BLOCK_SIZE;
    1193            final int jColumn = column - jBlock * BLOCK_SIZE;
    1194            final int jWidth  = blockWidth(jBlock);
    1195            int outIndex      = 0;
    1196            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1197                final int iHeight = blockHeight(iBlock);
    1198                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1199                for (int i = 0; i < iHeight; ++i) {
    1200                    out[outIndex++] = block[i * jWidth + jColumn];
    1201                }
    1202            }
    1203
    1204            return out;
    1205
    1206        }
    1207
    1208        /** {@inheritDoc} */
    1209        @Override
    1210        public void setColumn(final int column, final double[] array)
    1211            throws MatrixIndexException, InvalidMatrixException {
    1212
    1213            MatrixUtils.checkColumnIndex(this, column);
    1214            final int nRows = getRowDimension();
    1215            if (array.length != nRows) {
    1216                throw new InvalidMatrixException(
    1217                        LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    1218                        array.length, 1, nRows, 1);
    1219            }
    1220
    1221            // perform copy block-wise, to ensure good cache behavior
    1222            final int jBlock  = column / BLOCK_SIZE;
    1223            final int jColumn = column - jBlock * BLOCK_SIZE;
    1224            final int jWidth  = blockWidth(jBlock);
    1225            int outIndex      = 0;
    1226            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1227                final int iHeight = blockHeight(iBlock);
    1228                final double[] block = blocks[iBlock * blockColumns + jBlock];
    1229                for (int i = 0; i < iHeight; ++i) {
    1230                    block[i * jWidth + jColumn] = array[outIndex++];
    1231                }
    1232            }
    1233
    1234        }
    1235
    1236        /** {@inheritDoc} */
    1237        @Override
    1238        public double getEntry(final int row, final int column)
    1239            throws MatrixIndexException {
    1240            try {
    1241                final int iBlock = row    / BLOCK_SIZE;
    1242                final int jBlock = column / BLOCK_SIZE;
    1243                final int k      = (row    - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
    1244                                   (column - jBlock * BLOCK_SIZE);
    1245                return blocks[iBlock * blockColumns + jBlock][k];
    1246            } catch (ArrayIndexOutOfBoundsException e) {
    1247                throw new MatrixIndexException(
    1248                        LocalizedFormats.NO_SUCH_MATRIX_ENTRY,
    1249                        row, column, getRowDimension(), getColumnDimension());
    1250            }
    1251        }
    1252
    1253        /** {@inheritDoc} */
    1254        @Override
    1255        public void setEntry(final int row, final int column, final double value)
    1256            throws MatrixIndexException {
    1257            try {
    1258                final int iBlock = row    / BLOCK_SIZE;
    1259                final int jBlock = column / BLOCK_SIZE;
    1260                final int k      = (row    - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
    1261                                   (column - jBlock * BLOCK_SIZE);
    1262                blocks[iBlock * blockColumns + jBlock][k] = value;
    1263            } catch (ArrayIndexOutOfBoundsException e) {
    1264                throw new MatrixIndexException(
    1265                        LocalizedFormats.NO_SUCH_MATRIX_ENTRY,
    1266                        row, column, getRowDimension(), getColumnDimension());
    1267            }
    1268        }
    1269
    1270        /** {@inheritDoc} */
    1271        @Override
    1272        public void addToEntry(final int row, final int column, final double increment)
    1273            throws MatrixIndexException {
    1274            try {
    1275                final int iBlock = row    / BLOCK_SIZE;
    1276                final int jBlock = column / BLOCK_SIZE;
    1277                final int k      = (row    - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
    1278                                   (column - jBlock * BLOCK_SIZE);
    1279                blocks[iBlock * blockColumns + jBlock][k] += increment;
    1280            } catch (ArrayIndexOutOfBoundsException e) {
    1281                throw new MatrixIndexException(
    1282                        LocalizedFormats.NO_SUCH_MATRIX_ENTRY,
    1283                        row, column, getRowDimension(), getColumnDimension());
    1284            }
    1285        }
    1286
    1287        /** {@inheritDoc} */
    1288        @Override
    1289        public void multiplyEntry(final int row, final int column, final double factor)
    1290            throws MatrixIndexException {
    1291            try {
    1292                final int iBlock = row    / BLOCK_SIZE;
    1293                final int jBlock = column / BLOCK_SIZE;
    1294                final int k      = (row    - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
    1295                                   (column - jBlock * BLOCK_SIZE);
    1296                blocks[iBlock * blockColumns + jBlock][k] *= factor;
    1297            } catch (ArrayIndexOutOfBoundsException e) {
    1298                throw new MatrixIndexException(
    1299                        LocalizedFormats.NO_SUCH_MATRIX_ENTRY,
    1300                        row, column, getRowDimension(), getColumnDimension());
    1301            }
    1302        }
    1303
    1304        /** {@inheritDoc} */
    1305        @Override
    1306        public BlockRealMatrix transpose() {
    1307
    1308            final int nRows = getRowDimension();
    1309            final int nCols = getColumnDimension();
    1310            final BlockRealMatrix out = new BlockRealMatrix(nCols, nRows);
    1311
    1312            // perform transpose block-wise, to ensure good cache behavior
    1313            int blockIndex = 0;
    1314            for (int iBlock = 0; iBlock < blockColumns; ++iBlock) {
    1315                for (int jBlock = 0; jBlock < blockRows; ++jBlock) {
    1316
    1317                    // transpose current block
    1318                    final double[] outBlock = out.blocks[blockIndex];
    1319                    final double[] tBlock   = blocks[jBlock * blockColumns + iBlock];
    1320                    final int      pStart   = iBlock * BLOCK_SIZE;
    1321                    final int      pEnd     = FastMath.min(pStart + BLOCK_SIZE, columns);
    1322                    final int      qStart   = jBlock * BLOCK_SIZE;
    1323                    final int      qEnd     = FastMath.min(qStart + BLOCK_SIZE, rows);
    1324                    int k = 0;
    1325                    for (int p = pStart; p < pEnd; ++p) {
    1326                        final int lInc = pEnd - pStart;
    1327                        int l = p - pStart;
    1328                        for (int q = qStart; q < qEnd; ++q) {
    1329                            outBlock[k] = tBlock[l];
    1330                            ++k;
    1331                            l+= lInc;
    1332                        }
    1333                    }
    1334
    1335                    // go to next block
    1336                    ++blockIndex;
    1337
    1338                }
    1339            }
    1340
    1341            return out;
    1342
    1343        }
    1344
    1345        /** {@inheritDoc} */
    1346        @Override
    1347        public int getRowDimension() {
    1348            return rows;
    1349        }
    1350
    1351        /** {@inheritDoc} */
    1352        @Override
    1353        public int getColumnDimension() {
    1354            return columns;
    1355        }
    1356
    1357        /** {@inheritDoc} */
    1358        @Override
    1359        public double[] operate(final double[] v)
    1360            throws IllegalArgumentException {
    1361
    1362            if (v.length != columns) {
    1363                throw MathRuntimeException.createIllegalArgumentException(
    1364                        LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    1365                        v.length, columns);
    1366            }
    1367            final double[] out = new double[rows];
    1368
    1369            // perform multiplication block-wise, to ensure good cache behavior
    1370            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1371                final int pStart = iBlock * BLOCK_SIZE;
    1372                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1373                for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1374                    final double[] block  = blocks[iBlock * blockColumns + jBlock];
    1375                    final int      qStart = jBlock * BLOCK_SIZE;
    1376                    final int      qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1377                    int k = 0;
    1378                    for (int p = pStart; p < pEnd; ++p) {
    1379                        double sum = 0;
    1380                        int q = qStart;
    1381                        while (q < qEnd - 3) {
    1382                            sum += block[k]     * v[q]     +
    1383                                   block[k + 1] * v[q + 1] +
    1384                                   block[k + 2] * v[q + 2] +
    1385                                   block[k + 3] * v[q + 3];
    1386                            k += 4;
    1387                            q += 4;
    1388                        }
    1389                        while (q < qEnd) {
    1390                            sum += block[k++] * v[q++];
    1391                        }
    1392                        out[p] += sum;
    1393                    }
    1394                }
    1395            }
    1396
    1397            return out;
    1398
    1399        }
    1400
    1401        /** {@inheritDoc} */
    1402        @Override
    1403        public double[] preMultiply(final double[] v)
    1404            throws IllegalArgumentException {
    1405
    1406            if (v.length != rows) {
    1407                throw MathRuntimeException.createIllegalArgumentException(
    1408                        LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    1409                        v.length, rows);
    1410            }
    1411            final double[] out = new double[columns];
    1412
    1413            // perform multiplication block-wise, to ensure good cache behavior
    1414            for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1415                final int jWidth  = blockWidth(jBlock);
    1416                final int jWidth2 = jWidth  + jWidth;
    1417                final int jWidth3 = jWidth2 + jWidth;
    1418                final int jWidth4 = jWidth3 + jWidth;
    1419                final int qStart = jBlock * BLOCK_SIZE;
    1420                final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1421                for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1422                    final double[] block  = blocks[iBlock * blockColumns + jBlock];
    1423                    final int      pStart = iBlock * BLOCK_SIZE;
    1424                    final int      pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1425                    for (int q = qStart; q < qEnd; ++q) {
    1426                        int k = q - qStart;
    1427                        double sum = 0;
    1428                        int p = pStart;
    1429                        while (p < pEnd - 3) {
    1430                            sum += block[k]           * v[p]     +
    1431                                   block[k + jWidth]  * v[p + 1] +
    1432                                   block[k + jWidth2] * v[p + 2] +
    1433                                   block[k + jWidth3] * v[p + 3];
    1434                            k += jWidth4;
    1435                            p += 4;
    1436                        }
    1437                        while (p < pEnd) {
    1438                            sum += block[k] * v[p++];
    1439                            k += jWidth;
    1440                        }
    1441                        out[q] += sum;
    1442                    }
    1443                }
    1444            }
    1445
    1446            return out;
    1447
    1448        }
    1449
    1450        /** {@inheritDoc} */
    1451        @Override
    1452        public double walkInRowOrder(final RealMatrixChangingVisitor visitor)
    1453            throws MatrixVisitorException {
    1454            visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
    1455            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1456                final int pStart = iBlock * BLOCK_SIZE;
    1457                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1458                for (int p = pStart; p < pEnd; ++p) {
    1459                    for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1460                        final int jWidth = blockWidth(jBlock);
    1461                        final int qStart = jBlock * BLOCK_SIZE;
    1462                        final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1463                        final double[] block = blocks[iBlock * blockColumns + jBlock];
    1464                        int k = (p - pStart) * jWidth;
    1465                        for (int q = qStart; q < qEnd; ++q) {
    1466                            block[k] = visitor.visit(p, q, block[k]);
    1467                            ++k;
    1468                        }
    1469                    }
    1470                 }
    1471            }
    1472            return visitor.end();
    1473        }
    1474
    1475        /** {@inheritDoc} */
    1476        @Override
    1477        public double walkInRowOrder(final RealMatrixPreservingVisitor visitor)
    1478            throws MatrixVisitorException {
    1479            visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
    1480            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1481                final int pStart = iBlock * BLOCK_SIZE;
    1482                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1483                for (int p = pStart; p < pEnd; ++p) {
    1484                    for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1485                        final int jWidth = blockWidth(jBlock);
    1486                        final int qStart = jBlock * BLOCK_SIZE;
    1487                        final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1488                        final double[] block = blocks[iBlock * blockColumns + jBlock];
    1489                        int k = (p - pStart) * jWidth;
    1490                        for (int q = qStart; q < qEnd; ++q) {
    1491                            visitor.visit(p, q, block[k]);
    1492                            ++k;
    1493                        }
    1494                    }
    1495                 }
    1496            }
    1497            return visitor.end();
    1498        }
    1499
    1500        /** {@inheritDoc} */
    1501        @Override
    1502        public double walkInRowOrder(final RealMatrixChangingVisitor visitor,
    1503                                     final int startRow, final int endRow,
    1504                                     final int startColumn, final int endColumn)
    1505            throws MatrixIndexException, MatrixVisitorException {
    1506            MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    1507            visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    1508            for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
    1509                final int p0     = iBlock * BLOCK_SIZE;
    1510                final int pStart = FastMath.max(startRow, p0);
    1511                final int pEnd   = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
    1512                for (int p = pStart; p < pEnd; ++p) {
    1513                    for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
    1514                        final int jWidth = blockWidth(jBlock);
    1515                        final int q0     = jBlock * BLOCK_SIZE;
    1516                        final int qStart = FastMath.max(startColumn, q0);
    1517                        final int qEnd   = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
    1518                        final double[] block = blocks[iBlock * blockColumns + jBlock];
    1519                        int k = (p - p0) * jWidth + qStart - q0;
    1520                        for (int q = qStart; q < qEnd; ++q) {
    1521                            block[k] = visitor.visit(p, q, block[k]);
    1522                            ++k;
    1523                        }
    1524                    }
    1525                 }
    1526            }
    1527            return visitor.end();
    1528        }
    1529
    1530        /** {@inheritDoc} */
    1531        @Override
    1532        public double walkInRowOrder(final RealMatrixPreservingVisitor visitor,
    1533                                     final int startRow, final int endRow,
    1534                                     final int startColumn, final int endColumn)
    1535            throws MatrixIndexException, MatrixVisitorException {
    1536            MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    1537            visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    1538            for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
    1539                final int p0     = iBlock * BLOCK_SIZE;
    1540                final int pStart = FastMath.max(startRow, p0);
    1541                final int pEnd   = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
    1542                for (int p = pStart; p < pEnd; ++p) {
    1543                    for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
    1544                        final int jWidth = blockWidth(jBlock);
    1545                        final int q0     = jBlock * BLOCK_SIZE;
    1546                        final int qStart = FastMath.max(startColumn, q0);
    1547                        final int qEnd   = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
    1548                        final double[] block = blocks[iBlock * blockColumns + jBlock];
    1549                        int k = (p - p0) * jWidth + qStart - q0;
    1550                        for (int q = qStart; q < qEnd; ++q) {
    1551                            visitor.visit(p, q, block[k]);
    1552                            ++k;
    1553                        }
    1554                    }
    1555                 }
    1556            }
    1557            return visitor.end();
    1558        }
    1559
    1560        /** {@inheritDoc} */
    1561        @Override
    1562        public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor)
    1563            throws MatrixVisitorException {
    1564            visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
    1565            int blockIndex = 0;
    1566            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1567                final int pStart = iBlock * BLOCK_SIZE;
    1568                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1569                for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1570                    final int qStart = jBlock * BLOCK_SIZE;
    1571                    final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1572                    final double[] block = blocks[blockIndex];
    1573                    int k = 0;
    1574                    for (int p = pStart; p < pEnd; ++p) {
    1575                        for (int q = qStart; q < qEnd; ++q) {
    1576                            block[k] = visitor.visit(p, q, block[k]);
    1577                            ++k;
    1578                        }
    1579                    }
    1580                    ++blockIndex;
    1581                }
    1582            }
    1583            return visitor.end();
    1584        }
    1585
    1586        /** {@inheritDoc} */
    1587        @Override
    1588        public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor)
    1589            throws MatrixVisitorException {
    1590            visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
    1591            int blockIndex = 0;
    1592            for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
    1593                final int pStart = iBlock * BLOCK_SIZE;
    1594                final int pEnd   = FastMath.min(pStart + BLOCK_SIZE, rows);
    1595                for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
    1596                    final int qStart = jBlock * BLOCK_SIZE;
    1597                    final int qEnd   = FastMath.min(qStart + BLOCK_SIZE, columns);
    1598                    final double[] block = blocks[blockIndex];
    1599                    int k = 0;
    1600                    for (int p = pStart; p < pEnd; ++p) {
    1601                        for (int q = qStart; q < qEnd; ++q) {
    1602                            visitor.visit(p, q, block[k]);
    1603                            ++k;
    1604                        }
    1605                    }
    1606                    ++blockIndex;
    1607                }
    1608            }
    1609            return visitor.end();
    1610        }
    1611
    1612        /** {@inheritDoc} */
    1613        @Override
    1614        public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor,
    1615                                           final int startRow, final int endRow,
    1616                                           final int startColumn, final int endColumn)
    1617            throws MatrixIndexException, MatrixVisitorException {
    1618            MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    1619            visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    1620            for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
    1621                final int p0     = iBlock * BLOCK_SIZE;
    1622                final int pStart = FastMath.max(startRow, p0);
    1623                final int pEnd   = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
    1624                for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
    1625                    final int jWidth = blockWidth(jBlock);
    1626                    final int q0     = jBlock * BLOCK_SIZE;
    1627                    final int qStart = FastMath.max(startColumn, q0);
    1628                    final int qEnd   = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
    1629                    final double[] block = blocks[iBlock * blockColumns + jBlock];
    1630                    for (int p = pStart; p < pEnd; ++p) {
    1631                        int k = (p - p0) * jWidth + qStart - q0;
    1632                        for (int q = qStart; q < qEnd; ++q) {
    1633                            block[k] = visitor.visit(p, q, block[k]);
    1634                            ++k;
    1635                        }
    1636                    }
    1637                }
    1638            }
    1639            return visitor.end();
    1640        }
    1641
    1642        /** {@inheritDoc} */
    1643        @Override
    1644        public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor,
    1645                                           final int startRow, final int endRow,
    1646                                           final int startColumn, final int endColumn)
    1647            throws MatrixIndexException, MatrixVisitorException {
    1648            MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    1649            visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    1650            for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
    1651                final int p0     = iBlock * BLOCK_SIZE;
    1652                final int pStart = FastMath.max(startRow, p0);
    1653                final int pEnd   = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
    1654                for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
    1655                    final int jWidth = blockWidth(jBlock);
    1656                    final int q0     = jBlock * BLOCK_SIZE;
    1657                    final int qStart = FastMath.max(startColumn, q0);
    1658                    final int qEnd   = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
    1659                    final double[] block = blocks[iBlock * blockColumns + jBlock];
    1660                    for (int p = pStart; p < pEnd; ++p) {
    1661                        int k = (p - p0) * jWidth + qStart - q0;
    1662                        for (int q = qStart; q < qEnd; ++q) {
    1663                            visitor.visit(p, q, block[k]);
    1664                            ++k;
    1665                        }
    1666                    }
    1667                }
    1668            }
    1669            return visitor.end();
    1670        }
    1671
    1672        /**
    1673         * Get the height of a block.
    1674         * @param blockRow row index (in block sense) of the block
    1675         * @return height (number of rows) of the block
    1676         */
    1677        private int blockHeight(final int blockRow) {
    1678            return (blockRow == blockRows - 1) ? rows - blockRow * BLOCK_SIZE : BLOCK_SIZE;
    1679        }
    1680
    1681        /**
    1682         * Get the width of a block.
    1683         * @param blockColumn column index (in block sense) of the block
    1684         * @return width (number of columns) of the block
    1685         */
    1686        private int blockWidth(final int blockColumn) {
    1687            return (blockColumn == blockColumns - 1) ? columns - blockColumn * BLOCK_SIZE : BLOCK_SIZE;
    1688        }
    1689
    1690    }