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