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: 783741 $ $Date: 2009-06-11 08:35:42 -0400 (Thu, 11 Jun 2009) $
064 * @since 2.0
065 */
066 public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMatrix<T> implements Serializable {
067
068 /** Serializable version identifier */
069 private static final long serialVersionUID = -4602336630143123183L;
070
071 /** Block size. */
072 public static final int BLOCK_SIZE = 36;
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 for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
230 final int pStart = iBlock * BLOCK_SIZE;
231 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
232 final int iHeight = pEnd - pStart;
233 for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
234 final int qStart = jBlock * BLOCK_SIZE;
235 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
236 final int jWidth = qEnd - qStart;
237
238 // allocate new block
239 final T[] block = buildArray(field, iHeight * jWidth);
240 blocks[blockIndex] = block;
241
242 // copy data
243 for (int p = pStart, index = 0; p < pEnd; ++p, index += jWidth) {
244 System.arraycopy(rawData[p], qStart, block, index, jWidth);
245 }
246
247 }
248 }
249
250 return blocks;
251
252 }
253
254 /**
255 * Create a data array in blocks layout.
256 * <p>
257 * This method can be used to create the array argument of the {@link
258 * #BlockFieldMatrix(int, int, FieldElement[][], boolean)}
259 * constructor.
260 * </p>
261 * @param <T> the type of the field elements
262 * @param field field to which the elements belong
263 * @param rows the number of rows in the new matrix
264 * @param columns the number of columns in the new matrix
265 * @return a new data array in blocks layout
266 * @see #toBlocksLayout(FieldElement[][])
267 * @see #BlockFieldMatrix(int, int, FieldElement[][], boolean)
268 */
269 public static <T extends FieldElement<T>> T[][] createBlocksLayout(final Field<T> field,
270 final int rows, final int columns) {
271
272 final int blockRows = (rows + BLOCK_SIZE - 1) / BLOCK_SIZE;
273 final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE;
274
275 final T[][] blocks = buildArray(field, blockRows * blockColumns, -1);
276 for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
277 final int pStart = iBlock * BLOCK_SIZE;
278 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
279 final int iHeight = pEnd - pStart;
280 for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
281 final int qStart = jBlock * BLOCK_SIZE;
282 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
283 final int jWidth = qEnd - qStart;
284 blocks[blockIndex] = buildArray(field, iHeight * jWidth);
285 }
286 }
287
288 return blocks;
289
290 }
291
292 /** {@inheritDoc} */
293 @Override
294 public FieldMatrix<T> createMatrix(final int rowDimension, final int columnDimension)
295 throws IllegalArgumentException {
296 return new BlockFieldMatrix<T>(getField(), rowDimension, columnDimension);
297 }
298
299 /** {@inheritDoc} */
300 @Override
301 public FieldMatrix<T> copy() {
302
303 // create an empty matrix
304 BlockFieldMatrix<T> copied = new BlockFieldMatrix<T>(getField(), rows, columns);
305
306 // copy the blocks
307 for (int i = 0; i < blocks.length; ++i) {
308 System.arraycopy(blocks[i], 0, copied.blocks[i], 0, blocks[i].length);
309 }
310
311 return copied;
312
313 }
314
315 /** {@inheritDoc} */
316 @Override
317 public FieldMatrix<T> add(final FieldMatrix<T> m)
318 throws IllegalArgumentException {
319 try {
320 return add((BlockFieldMatrix<T>) m);
321 } catch (ClassCastException cce) {
322
323 // safety check
324 checkAdditionCompatible(m);
325
326 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
327
328 // perform addition block-wise, to ensure good cache behavior
329 int blockIndex = 0;
330 for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
331 for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
332
333 // perform addition on the current block
334 final T[] outBlock = out.blocks[blockIndex];
335 final T[] tBlock = blocks[blockIndex];
336 final int pStart = iBlock * BLOCK_SIZE;
337 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
338 final int qStart = jBlock * BLOCK_SIZE;
339 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
340 for (int p = pStart, k = 0; p < pEnd; ++p) {
341 for (int q = qStart; q < qEnd; ++q, ++k) {
342 outBlock[k] = tBlock[k].add(m.getEntry(p, q));
343 }
344 }
345
346 // go to next block
347 ++blockIndex;
348
349 }
350 }
351
352 return out;
353
354 }
355 }
356
357 /**
358 * Compute the sum of this and <code>m</code>.
359 *
360 * @param m matrix to be added
361 * @return this + m
362 * @throws IllegalArgumentException if m is not the same size as this
363 */
364 public BlockFieldMatrix<T> add(final BlockFieldMatrix<T> m)
365 throws IllegalArgumentException {
366
367 // safety check
368 checkAdditionCompatible(m);
369
370 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
371
372 // perform addition block-wise, to ensure good cache behavior
373 for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
374 final T[] outBlock = out.blocks[blockIndex];
375 final T[] tBlock = blocks[blockIndex];
376 final T[] mBlock = m.blocks[blockIndex];
377 for (int k = 0; k < outBlock.length; ++k) {
378 outBlock[k] = tBlock[k].add(mBlock[k]);
379 }
380 }
381
382 return out;
383
384 }
385
386 /** {@inheritDoc} */
387 @Override
388 public FieldMatrix<T> subtract(final FieldMatrix<T> m)
389 throws IllegalArgumentException {
390 try {
391 return subtract((BlockFieldMatrix<T>) m);
392 } catch (ClassCastException cce) {
393
394 // safety check
395 checkSubtractionCompatible(m);
396
397 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
398
399 // perform subtraction block-wise, to ensure good cache behavior
400 int blockIndex = 0;
401 for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
402 for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
403
404 // perform subtraction on the current block
405 final T[] outBlock = out.blocks[blockIndex];
406 final T[] tBlock = blocks[blockIndex];
407 final int pStart = iBlock * BLOCK_SIZE;
408 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
409 final int qStart = jBlock * BLOCK_SIZE;
410 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
411 for (int p = pStart, k = 0; p < pEnd; ++p) {
412 for (int q = qStart; q < qEnd; ++q, ++k) {
413 outBlock[k] = tBlock[k].subtract(m.getEntry(p, q));
414 }
415 }
416
417 // go to next block
418 ++blockIndex;
419
420 }
421 }
422
423 return out;
424
425 }
426 }
427
428 /**
429 * Compute this minus <code>m</code>.
430 *
431 * @param m matrix to be subtracted
432 * @return this - m
433 * @throws IllegalArgumentException if m is not the same size as this
434 */
435 public BlockFieldMatrix<T> subtract(final BlockFieldMatrix<T> m)
436 throws IllegalArgumentException {
437
438 // safety check
439 checkSubtractionCompatible(m);
440
441 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
442
443 // perform subtraction block-wise, to ensure good cache behavior
444 for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
445 final T[] outBlock = out.blocks[blockIndex];
446 final T[] tBlock = blocks[blockIndex];
447 final T[] mBlock = m.blocks[blockIndex];
448 for (int k = 0; k < outBlock.length; ++k) {
449 outBlock[k] = tBlock[k].subtract(mBlock[k]);
450 }
451 }
452
453 return out;
454
455 }
456
457 /** {@inheritDoc} */
458 @Override
459 public FieldMatrix<T> scalarAdd(final T d)
460 throws IllegalArgumentException {
461
462 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
463
464 // perform subtraction block-wise, to ensure good cache behavior
465 for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
466 final T[] outBlock = out.blocks[blockIndex];
467 final T[] tBlock = blocks[blockIndex];
468 for (int k = 0; k < outBlock.length; ++k) {
469 outBlock[k] = tBlock[k].add(d);
470 }
471 }
472
473 return out;
474
475 }
476
477 /** {@inheritDoc} */
478 @Override
479 public FieldMatrix<T> scalarMultiply(final T d)
480 throws IllegalArgumentException {
481
482 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, columns);
483
484 // perform subtraction block-wise, to ensure good cache behavior
485 for (int blockIndex = 0; blockIndex < out.blocks.length; ++blockIndex) {
486 final T[] outBlock = out.blocks[blockIndex];
487 final T[] tBlock = blocks[blockIndex];
488 for (int k = 0; k < outBlock.length; ++k) {
489 outBlock[k] = tBlock[k].multiply(d);
490 }
491 }
492
493 return out;
494
495 }
496
497 /** {@inheritDoc} */
498 @Override
499 public FieldMatrix<T> multiply(final FieldMatrix<T> m)
500 throws IllegalArgumentException {
501 try {
502 return multiply((BlockFieldMatrix<T>) m);
503 } catch (ClassCastException cce) {
504
505 // safety check
506 checkMultiplicationCompatible(m);
507
508 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, m.getColumnDimension());
509 final T zero = getField().getZero();
510
511 // perform multiplication block-wise, to ensure good cache behavior
512 int blockIndex = 0;
513 for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
514
515 final int pStart = iBlock * BLOCK_SIZE;
516 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
517
518 for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
519
520 final int qStart = jBlock * BLOCK_SIZE;
521 final int qEnd = Math.min(qStart + BLOCK_SIZE, m.getColumnDimension());
522
523 // select current block
524 final T[] outBlock = out.blocks[blockIndex];
525
526 // perform multiplication on current block
527 for (int kBlock = 0; kBlock < blockColumns; ++kBlock) {
528 final int kWidth = blockWidth(kBlock);
529 final T[] tBlock = blocks[iBlock * blockColumns + kBlock];
530 final int rStart = kBlock * BLOCK_SIZE;
531 for (int p = pStart, k = 0; p < pEnd; ++p) {
532 final int lStart = (p - pStart) * kWidth;
533 final int lEnd = lStart + kWidth;
534 for (int q = qStart; q < qEnd; ++q) {
535 T sum = zero;
536 for (int l = lStart, r = rStart; l < lEnd; ++l, ++r) {
537 sum = sum.add(tBlock[l].multiply(m.getEntry(r, q)));
538 }
539 outBlock[k] = outBlock[k].add(sum);
540 ++k;
541 }
542 }
543 }
544
545 // go to next block
546 ++blockIndex;
547
548 }
549 }
550
551 return out;
552
553 }
554 }
555
556 /**
557 * Returns the result of postmultiplying this by m.
558 *
559 * @param m matrix to postmultiply by
560 * @return this * m
561 * @throws IllegalArgumentException
562 * if columnDimension(this) != rowDimension(m)
563 */
564 public BlockFieldMatrix<T> multiply(BlockFieldMatrix<T> m) throws IllegalArgumentException {
565
566 // safety check
567 checkMultiplicationCompatible(m);
568
569 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, m.columns);
570 final T zero = getField().getZero();
571
572 // perform multiplication block-wise, to ensure good cache behavior
573 int blockIndex = 0;
574 for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) {
575
576 final int pStart = iBlock * BLOCK_SIZE;
577 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
578
579 for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) {
580 final int jWidth = out.blockWidth(jBlock);
581 final int jWidth2 = jWidth + jWidth;
582 final int jWidth3 = jWidth2 + jWidth;
583 final int jWidth4 = jWidth3 + jWidth;
584
585 // select current block
586 final T[] outBlock = out.blocks[blockIndex];
587
588 // perform multiplication on current block
589 for (int kBlock = 0; kBlock < blockColumns; ++kBlock) {
590 final int kWidth = blockWidth(kBlock);
591 final T[] tBlock = blocks[iBlock * blockColumns + kBlock];
592 final T[] mBlock = m.blocks[kBlock * m.blockColumns + jBlock];
593 for (int p = pStart, k = 0; p < pEnd; ++p) {
594 final int lStart = (p - pStart) * kWidth;
595 final int lEnd = lStart + kWidth;
596 for (int nStart = 0; nStart < jWidth; ++nStart) {
597 T sum = zero;
598 int l = lStart;
599 int n = nStart;
600 while (l < lEnd - 3) {
601 sum = sum.
602 add(tBlock[l].multiply(mBlock[n])).
603 add(tBlock[l + 1].multiply(mBlock[n + jWidth])).
604 add(tBlock[l + 2].multiply(mBlock[n + jWidth2])).
605 add(tBlock[l + 3].multiply(mBlock[n + jWidth3]));
606 l += 4;
607 n += jWidth4;
608 }
609 while (l < lEnd) {
610 sum = sum.add(tBlock[l++].multiply(mBlock[n]));
611 n += jWidth;
612 }
613 outBlock[k] = outBlock[k].add(sum);
614 ++k;
615 }
616 }
617 }
618
619 // go to next block
620 ++blockIndex;
621
622 }
623 }
624
625 return out;
626
627 }
628
629 /** {@inheritDoc} */
630 @Override
631 public T[][] getData() {
632
633 final T[][] data = buildArray(getField(), getRowDimension(), getColumnDimension());
634 final int lastColumns = columns - (blockColumns - 1) * BLOCK_SIZE;
635
636 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
637 final int pStart = iBlock * BLOCK_SIZE;
638 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
639 int regularPos = 0;
640 int lastPos = 0;
641 for (int p = pStart; p < pEnd; ++p) {
642 final T[] dataP = data[p];
643 int blockIndex = iBlock * blockColumns;
644 int dataPos = 0;
645 for (int jBlock = 0; jBlock < blockColumns - 1; ++jBlock) {
646 System.arraycopy(blocks[blockIndex++], regularPos, dataP, dataPos, BLOCK_SIZE);
647 dataPos += BLOCK_SIZE;
648 }
649 System.arraycopy(blocks[blockIndex], lastPos, dataP, dataPos, lastColumns);
650 regularPos += BLOCK_SIZE;
651 lastPos += lastColumns;
652 }
653 }
654
655 return data;
656
657 }
658
659 /** {@inheritDoc} */
660 @Override
661 public FieldMatrix<T> getSubMatrix(final int startRow, final int endRow,
662 final int startColumn, final int endColumn)
663 throws MatrixIndexException {
664
665 // safety checks
666 checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
667
668 // create the output matrix
669 final BlockFieldMatrix<T> out =
670 new BlockFieldMatrix<T>(getField(), endRow - startRow + 1, endColumn - startColumn + 1);
671
672 // compute blocks shifts
673 final int blockStartRow = startRow / BLOCK_SIZE;
674 final int rowsShift = startRow % BLOCK_SIZE;
675 final int blockStartColumn = startColumn / BLOCK_SIZE;
676 final int columnsShift = startColumn % BLOCK_SIZE;
677
678 // perform extraction block-wise, to ensure good cache behavior
679 for (int iBlock = 0, pBlock = blockStartRow; iBlock < out.blockRows; ++iBlock, ++pBlock) {
680 final int iHeight = out.blockHeight(iBlock);
681 for (int jBlock = 0, qBlock = blockStartColumn; jBlock < out.blockColumns; ++jBlock, ++qBlock) {
682 final int jWidth = out.blockWidth(jBlock);
683
684 // handle one block of the output matrix
685 final int outIndex = iBlock * out.blockColumns + jBlock;
686 final T[] outBlock = out.blocks[outIndex];
687 final int index = pBlock * blockColumns + qBlock;
688 final int width = blockWidth(qBlock);
689
690 final int heightExcess = iHeight + rowsShift - BLOCK_SIZE;
691 final int widthExcess = jWidth + columnsShift - BLOCK_SIZE;
692 if (heightExcess > 0) {
693 // the submatrix block spans on two blocks rows from the original matrix
694 if (widthExcess > 0) {
695 // the submatrix block spans on two blocks columns from the original matrix
696 final int width2 = blockWidth(qBlock + 1);
697 copyBlockPart(blocks[index], width,
698 rowsShift, BLOCK_SIZE,
699 columnsShift, BLOCK_SIZE,
700 outBlock, jWidth, 0, 0);
701 copyBlockPart(blocks[index + 1], width2,
702 rowsShift, BLOCK_SIZE,
703 0, widthExcess,
704 outBlock, jWidth, 0, jWidth - widthExcess);
705 copyBlockPart(blocks[index + blockColumns], width,
706 0, heightExcess,
707 columnsShift, BLOCK_SIZE,
708 outBlock, jWidth, iHeight - heightExcess, 0);
709 copyBlockPart(blocks[index + blockColumns + 1], width2,
710 0, heightExcess,
711 0, widthExcess,
712 outBlock, jWidth, iHeight - heightExcess, jWidth - widthExcess);
713 } else {
714 // the submatrix block spans on one block column from the original matrix
715 copyBlockPart(blocks[index], width,
716 rowsShift, BLOCK_SIZE,
717 columnsShift, jWidth + columnsShift,
718 outBlock, jWidth, 0, 0);
719 copyBlockPart(blocks[index + blockColumns], width,
720 0, heightExcess,
721 columnsShift, jWidth + columnsShift,
722 outBlock, jWidth, iHeight - heightExcess, 0);
723 }
724 } else {
725 // the submatrix block spans on one block row from the original matrix
726 if (widthExcess > 0) {
727 // the submatrix block spans on two blocks columns from the original matrix
728 final int width2 = blockWidth(qBlock + 1);
729 copyBlockPart(blocks[index], width,
730 rowsShift, iHeight + rowsShift,
731 columnsShift, BLOCK_SIZE,
732 outBlock, jWidth, 0, 0);
733 copyBlockPart(blocks[index + 1], width2,
734 rowsShift, iHeight + rowsShift,
735 0, widthExcess,
736 outBlock, jWidth, 0, jWidth - widthExcess);
737 } else {
738 // the submatrix block spans on one block column from the original matrix
739 copyBlockPart(blocks[index], width,
740 rowsShift, iHeight + rowsShift,
741 columnsShift, jWidth + columnsShift,
742 outBlock, jWidth, 0, 0);
743 }
744 }
745
746 }
747 }
748
749 return out;
750
751 }
752
753 /**
754 * Copy a part of a block into another one
755 * <p>This method can be called only when the specified part fits in both
756 * blocks, no verification is done here.</p>
757 * @param srcBlock source block
758 * @param srcWidth source block width ({@link #BLOCK_SIZE} or smaller)
759 * @param srcStartRow start row in the source block
760 * @param srcEndRow end row (exclusive) in the source block
761 * @param srcStartColumn start column in the source block
762 * @param srcEndColumn end column (exclusive) in the source block
763 * @param dstBlock destination block
764 * @param dstWidth destination block width ({@link #BLOCK_SIZE} or smaller)
765 * @param dstStartRow start row in the destination block
766 * @param dstStartColumn start column in the destination block
767 */
768 private void copyBlockPart(final T[] srcBlock, final int srcWidth,
769 final int srcStartRow, final int srcEndRow,
770 final int srcStartColumn, final int srcEndColumn,
771 final T[] dstBlock, final int dstWidth,
772 final int dstStartRow, final int dstStartColumn) {
773 final int length = srcEndColumn - srcStartColumn;
774 int srcPos = srcStartRow * srcWidth + srcStartColumn;
775 int dstPos = dstStartRow * dstWidth + dstStartColumn;
776 for (int srcRow = srcStartRow; srcRow < srcEndRow; ++srcRow) {
777 System.arraycopy(srcBlock, srcPos, dstBlock, dstPos, length);
778 srcPos += srcWidth;
779 dstPos += dstWidth;
780 }
781 }
782
783 /** {@inheritDoc} */
784 @Override
785 public void setSubMatrix(final T[][] subMatrix, final int row, final int column)
786 throws MatrixIndexException {
787
788 // safety checks
789 final int refLength = subMatrix[0].length;
790 if (refLength < 1) {
791 throw MathRuntimeException.createIllegalArgumentException("matrix must have at least one column");
792 }
793 final int endRow = row + subMatrix.length - 1;
794 final int endColumn = column + refLength - 1;
795 checkSubMatrixIndex(row, endRow, column, endColumn);
796 for (final T[] subRow : subMatrix) {
797 if (subRow.length != refLength) {
798 throw MathRuntimeException.createIllegalArgumentException(
799 "some rows have length {0} while others have length {1}",
800 refLength, subRow.length);
801 }
802 }
803
804 // compute blocks bounds
805 final int blockStartRow = row / BLOCK_SIZE;
806 final int blockEndRow = (endRow + BLOCK_SIZE) / BLOCK_SIZE;
807 final int blockStartColumn = column / BLOCK_SIZE;
808 final int blockEndColumn = (endColumn + BLOCK_SIZE) / BLOCK_SIZE;
809
810 // perform copy block-wise, to ensure good cache behavior
811 for (int iBlock = blockStartRow; iBlock < blockEndRow; ++iBlock) {
812 final int iHeight = blockHeight(iBlock);
813 final int firstRow = iBlock * BLOCK_SIZE;
814 final int iStart = Math.max(row, firstRow);
815 final int iEnd = Math.min(endRow + 1, firstRow + iHeight);
816
817 for (int jBlock = blockStartColumn; jBlock < blockEndColumn; ++jBlock) {
818 final int jWidth = blockWidth(jBlock);
819 final int firstColumn = jBlock * BLOCK_SIZE;
820 final int jStart = Math.max(column, firstColumn);
821 final int jEnd = Math.min(endColumn + 1, firstColumn + jWidth);
822 final int jLength = jEnd - jStart;
823
824 // handle one block, row by row
825 final T[] block = blocks[iBlock * blockColumns + jBlock];
826 for (int i = iStart; i < iEnd; ++i) {
827 System.arraycopy(subMatrix[i - row], jStart - column,
828 block, (i - firstRow) * jWidth + (jStart - firstColumn),
829 jLength);
830 }
831
832 }
833 }
834 }
835
836 /** {@inheritDoc} */
837 @Override
838 public FieldMatrix<T> getRowMatrix(final int row)
839 throws MatrixIndexException {
840
841 checkRowIndex(row);
842 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), 1, columns);
843
844 // perform copy block-wise, to ensure good cache behavior
845 final int iBlock = row / BLOCK_SIZE;
846 final int iRow = row - iBlock * BLOCK_SIZE;
847 int outBlockIndex = 0;
848 int outIndex = 0;
849 T[] outBlock = out.blocks[outBlockIndex];
850 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
851 final int jWidth = blockWidth(jBlock);
852 final T[] block = blocks[iBlock * blockColumns + jBlock];
853 final int available = outBlock.length - outIndex;
854 if (jWidth > available) {
855 System.arraycopy(block, iRow * jWidth, outBlock, outIndex, available);
856 outBlock = out.blocks[++outBlockIndex];
857 System.arraycopy(block, iRow * jWidth, outBlock, 0, jWidth - available);
858 outIndex = jWidth - available;
859 } else {
860 System.arraycopy(block, iRow * jWidth, outBlock, outIndex, jWidth);
861 outIndex += jWidth;
862 }
863 }
864
865 return out;
866
867 }
868
869 /** {@inheritDoc} */
870 @Override
871 public void setRowMatrix(final int row, final FieldMatrix<T> matrix)
872 throws MatrixIndexException, InvalidMatrixException {
873 try {
874 setRowMatrix(row, (BlockFieldMatrix<T>) matrix);
875 } catch (ClassCastException cce) {
876 super.setRowMatrix(row, matrix);
877 }
878 }
879
880 /**
881 * Sets the entries in row number <code>row</code>
882 * as a row matrix. Row indices start at 0.
883 *
884 * @param row the row to be set
885 * @param matrix row matrix (must have one row and the same number of columns
886 * as the instance)
887 * @throws MatrixIndexException if the specified row index is invalid
888 * @throws InvalidMatrixException if the matrix dimensions do not match one
889 * instance row
890 */
891 public void setRowMatrix(final int row, final BlockFieldMatrix<T> matrix)
892 throws MatrixIndexException, InvalidMatrixException {
893
894 checkRowIndex(row);
895 final int nCols = getColumnDimension();
896 if ((matrix.getRowDimension() != 1) ||
897 (matrix.getColumnDimension() != nCols)) {
898 throw new InvalidMatrixException(
899 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
900 matrix.getRowDimension(), matrix.getColumnDimension(),
901 1, nCols);
902 }
903
904 // perform copy block-wise, to ensure good cache behavior
905 final int iBlock = row / BLOCK_SIZE;
906 final int iRow = row - iBlock * BLOCK_SIZE;
907 int mBlockIndex = 0;
908 int mIndex = 0;
909 T[] mBlock = matrix.blocks[mBlockIndex];
910 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
911 final int jWidth = blockWidth(jBlock);
912 final T[] block = blocks[iBlock * blockColumns + jBlock];
913 final int available = mBlock.length - mIndex;
914 if (jWidth > available) {
915 System.arraycopy(mBlock, mIndex, block, iRow * jWidth, available);
916 mBlock = matrix.blocks[++mBlockIndex];
917 System.arraycopy(mBlock, 0, block, iRow * jWidth, jWidth - available);
918 mIndex = jWidth - available;
919 } else {
920 System.arraycopy(mBlock, mIndex, block, iRow * jWidth, jWidth);
921 mIndex += jWidth;
922 }
923 }
924
925 }
926
927 /** {@inheritDoc} */
928 @Override
929 public FieldMatrix<T> getColumnMatrix(final int column)
930 throws MatrixIndexException {
931
932 checkColumnIndex(column);
933 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), rows, 1);
934
935 // perform copy block-wise, to ensure good cache behavior
936 final int jBlock = column / BLOCK_SIZE;
937 final int jColumn = column - jBlock * BLOCK_SIZE;
938 final int jWidth = blockWidth(jBlock);
939 int outBlockIndex = 0;
940 int outIndex = 0;
941 T[] outBlock = out.blocks[outBlockIndex];
942 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
943 final int iHeight = blockHeight(iBlock);
944 final T[] block = blocks[iBlock * blockColumns + jBlock];
945 for (int i = 0; i < iHeight; ++i) {
946 if (outIndex >= outBlock.length) {
947 outBlock = out.blocks[++outBlockIndex];
948 outIndex = 0;
949 }
950 outBlock[outIndex++] = block[i * jWidth + jColumn];
951 }
952 }
953
954 return out;
955
956 }
957
958 /** {@inheritDoc} */
959 @Override
960 public void setColumnMatrix(final int column, final FieldMatrix<T> matrix)
961 throws MatrixIndexException, InvalidMatrixException {
962 try {
963 setColumnMatrix(column, (BlockFieldMatrix<T>) matrix);
964 } catch (ClassCastException cce) {
965 super.setColumnMatrix(column, matrix);
966 }
967 }
968
969 /**
970 * Sets the entries in column number <code>column</code>
971 * as a column matrix. Column indices start at 0.
972 *
973 * @param column the column to be set
974 * @param matrix column matrix (must have one column and the same number of rows
975 * as the instance)
976 * @throws MatrixIndexException if the specified column index is invalid
977 * @throws InvalidMatrixException if the matrix dimensions do not match one
978 * instance column
979 */
980 void setColumnMatrix(final int column, final BlockFieldMatrix<T> matrix)
981 throws MatrixIndexException, InvalidMatrixException {
982
983 checkColumnIndex(column);
984 final int nRows = getRowDimension();
985 if ((matrix.getRowDimension() != nRows) ||
986 (matrix.getColumnDimension() != 1)) {
987 throw new InvalidMatrixException(
988 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
989 matrix.getRowDimension(), matrix.getColumnDimension(),
990 nRows, 1);
991 }
992
993 // perform copy block-wise, to ensure good cache behavior
994 final int jBlock = column / BLOCK_SIZE;
995 final int jColumn = column - jBlock * BLOCK_SIZE;
996 final int jWidth = blockWidth(jBlock);
997 int mBlockIndex = 0;
998 int mIndex = 0;
999 T[] mBlock = matrix.blocks[mBlockIndex];
1000 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1001 final int iHeight = blockHeight(iBlock);
1002 final T[] block = blocks[iBlock * blockColumns + jBlock];
1003 for (int i = 0; i < iHeight; ++i) {
1004 if (mIndex >= mBlock.length) {
1005 mBlock = matrix.blocks[++mBlockIndex];
1006 mIndex = 0;
1007 }
1008 block[i * jWidth + jColumn] = mBlock[mIndex++];
1009 }
1010 }
1011
1012 }
1013
1014 /** {@inheritDoc} */
1015 @Override
1016 public FieldVector<T> getRowVector(final int row)
1017 throws MatrixIndexException {
1018
1019 checkRowIndex(row);
1020 final T[] outData = buildArray(getField(), columns);
1021
1022 // perform copy block-wise, to ensure good cache behavior
1023 final int iBlock = row / BLOCK_SIZE;
1024 final int iRow = row - iBlock * BLOCK_SIZE;
1025 int outIndex = 0;
1026 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1027 final int jWidth = blockWidth(jBlock);
1028 final T[] block = blocks[iBlock * blockColumns + jBlock];
1029 System.arraycopy(block, iRow * jWidth, outData, outIndex, jWidth);
1030 outIndex += jWidth;
1031 }
1032
1033 return new ArrayFieldVector<T>(outData, false);
1034
1035 }
1036
1037 /** {@inheritDoc} */
1038 @Override
1039 public void setRowVector(final int row, final FieldVector<T> vector)
1040 throws MatrixIndexException, InvalidMatrixException {
1041 try {
1042 setRow(row, ((ArrayFieldVector<T>) vector).getDataRef());
1043 } catch (ClassCastException cce) {
1044 super.setRowVector(row, vector);
1045 }
1046 }
1047
1048 /** {@inheritDoc} */
1049 @Override
1050 public FieldVector<T> getColumnVector(final int column)
1051 throws MatrixIndexException {
1052
1053 checkColumnIndex(column);
1054 final T[] outData = buildArray(getField(), rows);
1055
1056 // perform copy block-wise, to ensure good cache behavior
1057 final int jBlock = column / BLOCK_SIZE;
1058 final int jColumn = column - jBlock * BLOCK_SIZE;
1059 final int jWidth = blockWidth(jBlock);
1060 int outIndex = 0;
1061 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1062 final int iHeight = blockHeight(iBlock);
1063 final T[] block = blocks[iBlock * blockColumns + jBlock];
1064 for (int i = 0; i < iHeight; ++i) {
1065 outData[outIndex++] = block[i * jWidth + jColumn];
1066 }
1067 }
1068
1069 return new ArrayFieldVector<T>(outData, false);
1070
1071 }
1072
1073 /** {@inheritDoc} */
1074 @Override
1075 public void setColumnVector(final int column, final FieldVector<T> vector)
1076 throws MatrixIndexException, InvalidMatrixException {
1077 try {
1078 setColumn(column, ((ArrayFieldVector<T>) vector).getDataRef());
1079 } catch (ClassCastException cce) {
1080 super.setColumnVector(column, vector);
1081 }
1082 }
1083
1084 /** {@inheritDoc} */
1085 @Override
1086 public T[] getRow(final int row)
1087 throws MatrixIndexException {
1088
1089 checkRowIndex(row);
1090 final T[] out = buildArray(getField(), columns);
1091
1092 // perform copy block-wise, to ensure good cache behavior
1093 final int iBlock = row / BLOCK_SIZE;
1094 final int iRow = row - iBlock * BLOCK_SIZE;
1095 int outIndex = 0;
1096 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1097 final int jWidth = blockWidth(jBlock);
1098 final T[] block = blocks[iBlock * blockColumns + jBlock];
1099 System.arraycopy(block, iRow * jWidth, out, outIndex, jWidth);
1100 outIndex += jWidth;
1101 }
1102
1103 return out;
1104
1105 }
1106
1107 /** {@inheritDoc} */
1108 @Override
1109 public void setRow(final int row, final T[] array)
1110 throws MatrixIndexException, InvalidMatrixException {
1111
1112 checkRowIndex(row);
1113 final int nCols = getColumnDimension();
1114 if (array.length != nCols) {
1115 throw new InvalidMatrixException(
1116 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
1117 1, array.length, 1, nCols);
1118 }
1119
1120 // perform copy block-wise, to ensure good cache behavior
1121 final int iBlock = row / BLOCK_SIZE;
1122 final int iRow = row - iBlock * BLOCK_SIZE;
1123 int outIndex = 0;
1124 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1125 final int jWidth = blockWidth(jBlock);
1126 final T[] block = blocks[iBlock * blockColumns + jBlock];
1127 System.arraycopy(array, outIndex, block, iRow * jWidth, jWidth);
1128 outIndex += jWidth;
1129 }
1130
1131 }
1132
1133 /** {@inheritDoc} */
1134 @Override
1135 public T[] getColumn(final int column)
1136 throws MatrixIndexException {
1137
1138 checkColumnIndex(column);
1139 final T[] out = buildArray(getField(), rows);
1140
1141 // perform copy block-wise, to ensure good cache behavior
1142 final int jBlock = column / BLOCK_SIZE;
1143 final int jColumn = column - jBlock * BLOCK_SIZE;
1144 final int jWidth = blockWidth(jBlock);
1145 int outIndex = 0;
1146 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1147 final int iHeight = blockHeight(iBlock);
1148 final T[] block = blocks[iBlock * blockColumns + jBlock];
1149 for (int i = 0; i < iHeight; ++i) {
1150 out[outIndex++] = block[i * jWidth + jColumn];
1151 }
1152 }
1153
1154 return out;
1155
1156 }
1157
1158 /** {@inheritDoc} */
1159 @Override
1160 public void setColumn(final int column, final T[] array)
1161 throws MatrixIndexException, InvalidMatrixException {
1162
1163 checkColumnIndex(column);
1164 final int nRows = getRowDimension();
1165 if (array.length != nRows) {
1166 throw new InvalidMatrixException(
1167 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
1168 array.length, 1, nRows, 1);
1169 }
1170
1171 // perform copy block-wise, to ensure good cache behavior
1172 final int jBlock = column / BLOCK_SIZE;
1173 final int jColumn = column - jBlock * BLOCK_SIZE;
1174 final int jWidth = blockWidth(jBlock);
1175 int outIndex = 0;
1176 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1177 final int iHeight = blockHeight(iBlock);
1178 final T[] block = blocks[iBlock * blockColumns + jBlock];
1179 for (int i = 0; i < iHeight; ++i) {
1180 block[i * jWidth + jColumn] = array[outIndex++];
1181 }
1182 }
1183
1184 }
1185
1186 /** {@inheritDoc} */
1187 @Override
1188 public T getEntry(final int row, final int column)
1189 throws MatrixIndexException {
1190 try {
1191 final int iBlock = row / BLOCK_SIZE;
1192 final int jBlock = column / BLOCK_SIZE;
1193 final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
1194 (column - jBlock * BLOCK_SIZE);
1195 return blocks[iBlock * blockColumns + jBlock][k];
1196 } catch (ArrayIndexOutOfBoundsException e) {
1197 throw new MatrixIndexException(
1198 "no entry at indices ({0}, {1}) in a {2}x{3} matrix",
1199 row, column, getRowDimension(), getColumnDimension());
1200 }
1201 }
1202
1203 /** {@inheritDoc} */
1204 @Override
1205 public void setEntry(final int row, final int column, final T value)
1206 throws MatrixIndexException {
1207 try {
1208 final int iBlock = row / BLOCK_SIZE;
1209 final int jBlock = column / BLOCK_SIZE;
1210 final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
1211 (column - jBlock * BLOCK_SIZE);
1212 blocks[iBlock * blockColumns + jBlock][k] = value;
1213 } catch (ArrayIndexOutOfBoundsException e) {
1214 throw new MatrixIndexException(
1215 "no entry at indices ({0}, {1}) in a {2}x{3} matrix",
1216 row, column, getRowDimension(), getColumnDimension());
1217 }
1218 }
1219
1220 /** {@inheritDoc} */
1221 @Override
1222 public void addToEntry(final int row, final int column, final T increment)
1223 throws MatrixIndexException {
1224 try {
1225 final int iBlock = row / BLOCK_SIZE;
1226 final int jBlock = column / BLOCK_SIZE;
1227 final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
1228 (column - jBlock * BLOCK_SIZE);
1229 final T[] blockIJ = blocks[iBlock * blockColumns + jBlock];
1230 blockIJ[k] = blockIJ[k].add(increment);
1231 } catch (ArrayIndexOutOfBoundsException e) {
1232 throw new MatrixIndexException(
1233 "no entry at indices ({0}, {1}) in a {2}x{3} matrix",
1234 row, column, getRowDimension(), getColumnDimension());
1235 }
1236 }
1237
1238 /** {@inheritDoc} */
1239 @Override
1240 public void multiplyEntry(final int row, final int column, final T factor)
1241 throws MatrixIndexException {
1242 try {
1243 final int iBlock = row / BLOCK_SIZE;
1244 final int jBlock = column / BLOCK_SIZE;
1245 final int k = (row - iBlock * BLOCK_SIZE) * blockWidth(jBlock) +
1246 (column - jBlock * BLOCK_SIZE);
1247 final T[] blockIJ = blocks[iBlock * blockColumns + jBlock];
1248 blockIJ[k] = blockIJ[k].multiply(factor);
1249 } catch (ArrayIndexOutOfBoundsException e) {
1250 throw new MatrixIndexException(
1251 "no entry at indices ({0}, {1}) in a {2}x{3} matrix",
1252 row, column, getRowDimension(), getColumnDimension());
1253 }
1254 }
1255
1256 /** {@inheritDoc} */
1257 @Override
1258 public FieldMatrix<T> transpose() {
1259
1260 final int nRows = getRowDimension();
1261 final int nCols = getColumnDimension();
1262 final BlockFieldMatrix<T> out = new BlockFieldMatrix<T>(getField(), nCols, nRows);
1263
1264 // perform transpose block-wise, to ensure good cache behavior
1265 int blockIndex = 0;
1266 for (int iBlock = 0; iBlock < blockColumns; ++iBlock) {
1267 for (int jBlock = 0; jBlock < blockRows; ++jBlock) {
1268
1269 // transpose current block
1270 final T[] outBlock = out.blocks[blockIndex];
1271 final T[] tBlock = blocks[jBlock * blockColumns + iBlock];
1272 final int pStart = iBlock * BLOCK_SIZE;
1273 final int pEnd = Math.min(pStart + BLOCK_SIZE, columns);
1274 final int qStart = jBlock * BLOCK_SIZE;
1275 final int qEnd = Math.min(qStart + BLOCK_SIZE, rows);
1276 for (int p = pStart, k = 0; p < pEnd; ++p) {
1277 final int lInc = pEnd - pStart;
1278 for (int q = qStart, l = p - pStart; q < qEnd; ++q, l+= lInc) {
1279 outBlock[k++] = tBlock[l];
1280 }
1281 }
1282
1283 // go to next block
1284 ++blockIndex;
1285
1286 }
1287 }
1288
1289 return out;
1290
1291 }
1292
1293 /** {@inheritDoc} */
1294 @Override
1295 public int getRowDimension() {
1296 return rows;
1297 }
1298
1299 /** {@inheritDoc} */
1300 @Override
1301 public int getColumnDimension() {
1302 return columns;
1303 }
1304
1305 /** {@inheritDoc} */
1306 @Override
1307 public T[] operate(final T[] v)
1308 throws IllegalArgumentException {
1309
1310 if (v.length != columns) {
1311 throw MathRuntimeException.createIllegalArgumentException(
1312 "vector length mismatch: got {0} but expected {1}",
1313 v.length, columns);
1314 }
1315 final T[] out = buildArray(getField(), rows);
1316 final T zero = getField().getZero();
1317
1318 // perform multiplication block-wise, to ensure good cache behavior
1319 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1320 final int pStart = iBlock * BLOCK_SIZE;
1321 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1322 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1323 final T[] block = blocks[iBlock * blockColumns + jBlock];
1324 final int qStart = jBlock * BLOCK_SIZE;
1325 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1326 for (int p = pStart, k = 0; p < pEnd; ++p) {
1327 T sum = zero;
1328 int q = qStart;
1329 while (q < qEnd - 3) {
1330 sum = sum.
1331 add(block[k].multiply(v[q])).
1332 add(block[k + 1].multiply(v[q + 1])).
1333 add(block[k + 2].multiply(v[q + 2])).
1334 add(block[k + 3].multiply(v[q + 3]));
1335 k += 4;
1336 q += 4;
1337 }
1338 while (q < qEnd) {
1339 sum = sum.add(block[k++].multiply(v[q++]));
1340 }
1341 out[p] = out[p].add(sum);
1342 }
1343 }
1344 }
1345
1346 return out;
1347
1348 }
1349
1350 /** {@inheritDoc} */
1351 @Override
1352 public T[] preMultiply(final T[] v)
1353 throws IllegalArgumentException {
1354
1355 if (v.length != rows) {
1356 throw MathRuntimeException.createIllegalArgumentException(
1357 "vector length mismatch: got {0} but expected {1}",
1358 v.length, rows);
1359 }
1360 final T[] out = buildArray(getField(), columns);
1361 final T zero = getField().getZero();
1362
1363 // perform multiplication block-wise, to ensure good cache behavior
1364 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1365 final int jWidth = blockWidth(jBlock);
1366 final int jWidth2 = jWidth + jWidth;
1367 final int jWidth3 = jWidth2 + jWidth;
1368 final int jWidth4 = jWidth3 + jWidth;
1369 final int qStart = jBlock * BLOCK_SIZE;
1370 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1371 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1372 final T[] block = blocks[iBlock * blockColumns + jBlock];
1373 final int pStart = iBlock * BLOCK_SIZE;
1374 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1375 for (int q = qStart; q < qEnd; ++q) {
1376 int k = q - qStart;
1377 T sum = zero;
1378 int p = pStart;
1379 while (p < pEnd - 3) {
1380 sum = sum.
1381 add(block[k].multiply(v[p])).
1382 add(block[k + jWidth].multiply(v[p + 1])).
1383 add(block[k + jWidth2].multiply(v[p + 2])).
1384 add(block[k + jWidth3].multiply(v[p + 3]));
1385 k += jWidth4;
1386 p += 4;
1387 }
1388 while (p < pEnd) {
1389 sum = sum.add(block[k].multiply(v[p++]));
1390 k += jWidth;
1391 }
1392 out[q] = out[q].add(sum);
1393 }
1394 }
1395 }
1396
1397 return out;
1398
1399 }
1400
1401 /** {@inheritDoc} */
1402 @Override
1403 public T walkInRowOrder(final FieldMatrixChangingVisitor<T> visitor)
1404 throws MatrixVisitorException {
1405 visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
1406 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1407 final int pStart = iBlock * BLOCK_SIZE;
1408 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1409 for (int p = pStart; p < pEnd; ++p) {
1410 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1411 final int jWidth = blockWidth(jBlock);
1412 final int qStart = jBlock * BLOCK_SIZE;
1413 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1414 final T[] block = blocks[iBlock * blockColumns + jBlock];
1415 for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) {
1416 block[k] = visitor.visit(p, q, block[k]);
1417 }
1418 }
1419 }
1420 }
1421 return visitor.end();
1422 }
1423
1424 /** {@inheritDoc} */
1425 @Override
1426 public T walkInRowOrder(final FieldMatrixPreservingVisitor<T> visitor)
1427 throws MatrixVisitorException {
1428 visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
1429 for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
1430 final int pStart = iBlock * BLOCK_SIZE;
1431 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1432 for (int p = pStart; p < pEnd; ++p) {
1433 for (int jBlock = 0; jBlock < blockColumns; ++jBlock) {
1434 final int jWidth = blockWidth(jBlock);
1435 final int qStart = jBlock * BLOCK_SIZE;
1436 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1437 final T[] block = blocks[iBlock * blockColumns + jBlock];
1438 for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) {
1439 visitor.visit(p, q, block[k]);
1440 }
1441 }
1442 }
1443 }
1444 return visitor.end();
1445 }
1446
1447 /** {@inheritDoc} */
1448 @Override
1449 public T walkInRowOrder(final FieldMatrixChangingVisitor<T> visitor,
1450 final int startRow, final int endRow,
1451 final int startColumn, final int endColumn)
1452 throws MatrixIndexException, MatrixVisitorException {
1453 checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
1454 visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
1455 for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
1456 final int p0 = iBlock * BLOCK_SIZE;
1457 final int pStart = Math.max(startRow, p0);
1458 final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
1459 for (int p = pStart; p < pEnd; ++p) {
1460 for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
1461 final int jWidth = blockWidth(jBlock);
1462 final int q0 = jBlock * BLOCK_SIZE;
1463 final int qStart = Math.max(startColumn, q0);
1464 final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
1465 final T[] block = blocks[iBlock * blockColumns + jBlock];
1466 for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
1467 block[k] = visitor.visit(p, q, block[k]);
1468 }
1469 }
1470 }
1471 }
1472 return visitor.end();
1473 }
1474
1475 /** {@inheritDoc} */
1476 @Override
1477 public T walkInRowOrder(final FieldMatrixPreservingVisitor<T> visitor,
1478 final int startRow, final int endRow,
1479 final int startColumn, final int endColumn)
1480 throws MatrixIndexException, MatrixVisitorException {
1481 checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
1482 visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
1483 for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
1484 final int p0 = iBlock * BLOCK_SIZE;
1485 final int pStart = Math.max(startRow, p0);
1486 final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
1487 for (int p = pStart; p < pEnd; ++p) {
1488 for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
1489 final int jWidth = blockWidth(jBlock);
1490 final int q0 = jBlock * BLOCK_SIZE;
1491 final int qStart = Math.max(startColumn, q0);
1492 final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
1493 final T[] block = blocks[iBlock * blockColumns + jBlock];
1494 for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
1495 visitor.visit(p, q, block[k]);
1496 }
1497 }
1498 }
1499 }
1500 return visitor.end();
1501 }
1502
1503 /** {@inheritDoc} */
1504 @Override
1505 public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> visitor)
1506 throws MatrixVisitorException {
1507 visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
1508 for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
1509 final int pStart = iBlock * BLOCK_SIZE;
1510 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1511 for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
1512 final int qStart = jBlock * BLOCK_SIZE;
1513 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1514 final T[] block = blocks[blockIndex];
1515 for (int p = pStart, k = 0; p < pEnd; ++p) {
1516 for (int q = qStart; q < qEnd; ++q, ++k) {
1517 block[k] = visitor.visit(p, q, block[k]);
1518 }
1519 }
1520 }
1521 }
1522 return visitor.end();
1523 }
1524
1525 /** {@inheritDoc} */
1526 @Override
1527 public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> visitor)
1528 throws MatrixVisitorException {
1529 visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
1530 for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) {
1531 final int pStart = iBlock * BLOCK_SIZE;
1532 final int pEnd = Math.min(pStart + BLOCK_SIZE, rows);
1533 for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) {
1534 final int qStart = jBlock * BLOCK_SIZE;
1535 final int qEnd = Math.min(qStart + BLOCK_SIZE, columns);
1536 final T[] block = blocks[blockIndex];
1537 for (int p = pStart, k = 0; p < pEnd; ++p) {
1538 for (int q = qStart; q < qEnd; ++q, ++k) {
1539 visitor.visit(p, q, block[k]);
1540 }
1541 }
1542 }
1543 }
1544 return visitor.end();
1545 }
1546
1547 /** {@inheritDoc} */
1548 @Override
1549 public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> visitor,
1550 final int startRow, final int endRow,
1551 final int startColumn, final int endColumn)
1552 throws MatrixIndexException, MatrixVisitorException {
1553 checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
1554 visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
1555 for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
1556 final int p0 = iBlock * BLOCK_SIZE;
1557 final int pStart = Math.max(startRow, p0);
1558 final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
1559 for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
1560 final int jWidth = blockWidth(jBlock);
1561 final int q0 = jBlock * BLOCK_SIZE;
1562 final int qStart = Math.max(startColumn, q0);
1563 final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
1564 final T[] block = blocks[iBlock * blockColumns + jBlock];
1565 for (int p = pStart; p < pEnd; ++p) {
1566 for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
1567 block[k] = visitor.visit(p, q, block[k]);
1568 }
1569 }
1570 }
1571 }
1572 return visitor.end();
1573 }
1574
1575 /** {@inheritDoc} */
1576 @Override
1577 public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> visitor,
1578 final int startRow, final int endRow,
1579 final int startColumn, final int endColumn)
1580 throws MatrixIndexException, MatrixVisitorException {
1581 checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
1582 visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
1583 for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
1584 final int p0 = iBlock * BLOCK_SIZE;
1585 final int pStart = Math.max(startRow, p0);
1586 final int pEnd = Math.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
1587 for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
1588 final int jWidth = blockWidth(jBlock);
1589 final int q0 = jBlock * BLOCK_SIZE;
1590 final int qStart = Math.max(startColumn, q0);
1591 final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
1592 final T[] block = blocks[iBlock * blockColumns + jBlock];
1593 for (int p = pStart; p < pEnd; ++p) {
1594 for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) {
1595 visitor.visit(p, q, block[k]);
1596 }
1597 }
1598 }
1599 }
1600 return visitor.end();
1601 }
1602
1603 /**
1604 * Get the height of a block.
1605 * @param blockRow row index (in block sense) of the block
1606 * @return height (number of rows) of the block
1607 */
1608 private int blockHeight(final int blockRow) {
1609 return (blockRow == blockRows - 1) ? rows - blockRow * BLOCK_SIZE : BLOCK_SIZE;
1610 }
1611
1612 /**
1613 * Get the width of a block.
1614 * @param blockColumn column index (in block sense) of the block
1615 * @return width (number of columns) of the block
1616 */
1617 private int blockWidth(final int blockColumn) {
1618 return (blockColumn == blockColumns - 1) ? columns - blockColumn * BLOCK_SIZE : BLOCK_SIZE;
1619 }
1620
1621 }