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