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