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