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 }