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