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