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