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