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