001    /*
002     * Copyright (C) 2010 by Frode Carlsen, Holger Arndt
003     *
004     * This file is part of the Universal Java Matrix Package (UJMP).
005     * See the NOTICE file distributed with this work for additional
006     * information regarding copyright ownership and licensing.
007     *
008     * UJMP is free software; you can redistribute it and/or modify
009     * it under the terms of the GNU Lesser General Public License as
010     * published by the Free Software Foundation; either version 2
011     * of the License, or (at your option) any later version.
012     *
013     * UJMP is distributed in the hope that it will be useful,
014     * but WITHOUT ANY WARRANTY; without even the implied warranty of
015     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
016     * GNU Lesser General Public License for more details.
017     *
018     * You should have received a copy of the GNU Lesser General Public
019     * License along with UJMP; if not, write to the
020     * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
021     * Boston, MA  02110-1301  USA
022     */
023    package org.ujmp.core.doublematrix.impl;
024    
025    import static org.ujmp.core.util.VerifyUtil.assertTrue;
026    
027    import java.util.Arrays;
028    
029    import org.ujmp.core.Coordinates;
030    import org.ujmp.core.Matrix;
031    import org.ujmp.core.annotation.Annotation;
032    import org.ujmp.core.calculation.Mtimes;
033    import org.ujmp.core.calculation.Calculation.Ret;
034    import org.ujmp.core.doublematrix.DenseDoubleMatrix2D;
035    import org.ujmp.core.doublematrix.impl.BlockMatrixLayout.BlockOrder;
036    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
037    import org.ujmp.core.exceptions.MatrixException;
038    import org.ujmp.core.interfaces.HasBlockDoubleArray2D;
039    import org.ujmp.core.objectmatrix.calculation.Transpose;
040    import org.ujmp.core.util.UJMPSettings;
041    import org.ujmp.core.util.concurrent.PFor;
042    
043    /**
044     * A dense 2D matrix with square block layout. The data in the matrix is
045     * re-organised as a block, tiled layout, which on modern CPUs with multiple
046     * caches reduces the number of cache misses, by providing better cache locality
047     * and cache temporality.
048     * 
049     * <h4>Block Layout (example)</h4>
050     * <p>
051     * Example: a 4x4 matrix with block size 2x2 is internally re-organised into 4
052     * separate blocks :
053     * <p>
054     * <code>
055     * |(0,0) , (0, 1) | (0, 2) , (0, 3) |<br>
056     * |(1,0) , (1, 1) | (1, 2) , (1, 3) |<br>
057     * |----------------------------------<br>
058     * |(2,0) , (2, 1) | (2, 2) , (2, 3) |<br>
059     * |(3,0) , (3, 1) | (3, 2) , (3, 3) |<br>
060     * </code>
061     * <p>
062     * This can be described as a matrix of the following blocks: <br>
063     * <p>
064     * <code>
065     * | A | B |<br>
066     * | C | D |<br>
067     * </code>
068     * <p>
069     * Each block is mapped to a separate one-dimensional array. For example block
070     * A:
071     * <p>
072     * <code>
073     * [ (0,0), (1,0), (0, 1), (1, 1)]
074     * </code>
075     * <p>
076     * <p>
077     * This layout is similar to what's described in [ref: II fig. 4b].
078     * 
079     * <h4>Choice of block size:</h4>
080     * <p>
081     * The blocks may be square and can be configured by the user at runtime.
082     * However, matrix multiplication performance will be sensitive to the choice of
083     * block size, depending on the amount of CPU cache available on the system.
084     * <p>
085     * 
086     * @see I. 'Efficient Matrix Multiplication Using Cache Conscious Data Layouts';
087     *      Neungsoo Park, Wenheng Liuy, Viktor K. Prasanna, Cauligi Raghavendra;
088     *      Department of Electrical Engineering{Systems), University of Southern
089     *      California
090     * @see II. 'Tiling, Block Data Layout, and Memory Hierarchy Performance';
091     *      Neungsoo Park, Bo Hong, and Viktor K. Prasanna
092     * 
093     * @author Frode Carlsen
094     * @author Holger Arndt
095     */
096    public class BlockDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D implements
097                    HasBlockDoubleArray2D {
098    
099            private static final long serialVersionUID = -5131649082019624021L;
100    
101            private static int deriveDefaultBlockStripeSize(int rows, int cols) {
102                    // TODO pick a suitable size
103                    return (rows < UJMPSettings.getDefaultBlockSize() && cols < UJMPSettings
104                                    .getDefaultBlockSize()) ? 50 : UJMPSettings.getDefaultBlockSize();
105            }
106    
107            /** Matrix data by block number. */
108            private double[][] data;
109    
110            /** Layout of matrix and blocks. */
111            protected BlockMatrixLayout layout;
112    
113            /** Dimensions of the matrix. */
114            private final long[] size;
115    
116            public BlockDenseDoubleMatrix2D(final double[][] values) {
117                    this(values.length, values[0].length, BlockOrder.ROWMAJOR);
118                    fill(values);
119            }
120    
121            /**
122             * Create a block matrix from a jagged array. <br>
123             * All rows of the values array must have the same length.
124             * 
125             * @param values
126             *            - the data to populate the matrix with.
127             * @param blockStripeSize
128             *            - length of one side of a block
129             * @param blockOrder
130             *            - see {@link BlockOrder}.
131             * @throws NullPointerException
132             *             if values is null, or values[0] is null.
133             * @throws ArrayIndexOutOfBoundsException
134             *             if any row is shorter than the first row.
135             */
136            public BlockDenseDoubleMatrix2D(final double[][] values, final int blockStripeSize,
137                            final BlockOrder blockOrder) {
138                    this(values.length, values[0].length, blockStripeSize, blockOrder);
139                    fill(values);
140            }
141    
142            /**
143             * Create a new matrix with the specified size, and specified block stripe
144             * size.
145             * 
146             * @param rows
147             *            - number of rows of the matrix.
148             * @param cols
149             *            - number of columns of the matrix.
150             * @param blockStripeSize
151             *            - length of one side of a square block.
152             * @throws IllegalArgumentException
153             *             if rows, cols or blockStripeSize are 0 or less, or blockOrder
154             *             is null.
155             */
156            public BlockDenseDoubleMatrix2D(final int rows, final int cols, final int blockStripeSize,
157                            final BlockOrder blockOrder) {
158                    assertTrue(rows > 0, "rows<=0");
159                    assertTrue(cols > 0, "cols<=0");
160                    assertTrue(blockStripeSize > 0, "blockStripeSize<=0");
161                    assertTrue(blockOrder != null, "blockOrder == null");
162    
163                    this.size = new long[] { rows, cols };
164    
165                    // layout structure for the blocks in the matrix
166                    this.layout = new BlockMatrixLayout(rows, cols, blockStripeSize, blockOrder);
167    
168                    this.data = new double[this.layout.numberOfBlocks][];
169            }
170    
171            /**
172             * Create a new matrix with the given size (rows, cols) and block layout.
173             * 
174             * @see #BlockMatrix(int, int)
175             */
176            public BlockDenseDoubleMatrix2D(final long rows, final long cols, final BlockOrder blockOrder) {
177                    this((int) rows, (int) cols, deriveDefaultBlockStripeSize((int) rows, (int) cols),
178                                    blockOrder);
179            }
180    
181            /**
182             * Create a new matrix with the specified size, and specified block stripe
183             * size.
184             * 
185             * @param rows
186             *            - number of rows of the matrix.
187             * @param cols
188             *            - number of columns of the matrix.
189             * @throws IllegalArgumentException
190             *             if rows, cols are 0 or less.
191             */
192            public BlockDenseDoubleMatrix2D(final int rows, final int cols) {
193                    this(rows, cols, deriveDefaultBlockStripeSize(rows, cols), BlockOrder.ROWMAJOR);
194            }
195    
196            /**
197             * Create a new matrix with the given size (rows, cols) and default block
198             * layout.
199             */
200            public BlockDenseDoubleMatrix2D(final long... size) {
201                    this(size[ROW], size[COLUMN], BlockOrder.ROWMAJOR);
202            }
203    
204            /**
205             * Constructor which takes an existing Matrix to copy data and structure
206             * from. <br>
207             * Block stripe size will be defaulted internally.
208             * 
209             * @param m
210             *            - matrix to copy data from.
211             */
212            public BlockDenseDoubleMatrix2D(Matrix m) {
213                    this(m, deriveDefaultBlockStripeSize((int) m.getRowCount(), (int) m.getColumnCount()));
214            }
215    
216            /**
217             * Constructor which takes an existing BlockMatrix to copy data and
218             * structure from. <br>
219             */
220            public BlockDenseDoubleMatrix2D(final BlockDenseDoubleMatrix2D m) {
221                    this((int) m.size[ROW], (int) m.size[COLUMN], m.layout.blockStripe, m.layout.blockOrder);
222                    for (int i = m.layout.numberOfBlocks; --i != -1;) {
223                            final double[] block = m.data[i];
224                            if (block != null) {
225                                    // cannot use Arrays.copyOf(): not supported in Java 5
226                                    this.data[i] = new double[block.length];
227                                    System.arraycopy(block, 0, this.data[i], 0, block.length);
228                            }
229                    }
230                    Annotation a = m.getAnnotation();
231                    if (a != null) {
232                            setAnnotation(a.clone());
233                    }
234            }
235    
236            /**
237             * Constructor which takes a Matrix and a proposed default block stripe
238             * size.
239             * 
240             * @param m
241             *            - matrix containing existing values.
242             * @param blockStripeSize
243             *            - proposed default block size.
244             */
245            public BlockDenseDoubleMatrix2D(Matrix m, int blockStripeSize) {
246                    this(m, blockStripeSize, BlockOrder.ROWMAJOR);
247            }
248    
249            /**
250             * Constructor which takes a Matrix and a proposed default block stripe
251             * size.
252             * 
253             * @param m
254             *            - matrix containing existing values.
255             * @param blockStripeSize
256             *            - proposed default block size.
257             * @param blockOrder
258             *            row major or column major
259             */
260            public BlockDenseDoubleMatrix2D(Matrix m, int blockStripeSize, BlockOrder blockOrder) {
261                    this((int) m.getRowCount(), (int) m.getColumnCount(), blockStripeSize, blockOrder);
262    
263                    if (m instanceof DenseDoubleMatrix2D) {
264                            final DenseDoubleMatrix2D mDense = (DenseDoubleMatrix2D) m;
265                            final int mRows = (int) mDense.getRowCount(), mColumns = (int) mDense.getColumnCount();
266                            for (int j = 0; j < mColumns; j++) {
267                                    for (int i = 0; i < mRows; i++) {
268                                            setDouble(mDense.getDouble(i, j), i, j);
269                                    }
270                            }
271                    } else {
272                            for (long[] c : m.availableCoordinates()) {
273                                    setDouble(m.getAsDouble(c), c);
274                            }
275                    }
276            }
277    
278            protected void addBlockData(final int row, final int column, final double[] newData) {
279                    int blockNumber = layout.getBlockNumber(row, column);
280    
281                    if (null == data[blockNumber]) {
282                            // init first block
283                            synchronized (data) {
284                                    data[blockNumber] = newData;
285                                    return;
286                            }
287                    }
288    
289                    final double[] block = data[blockNumber];
290                    synchronized (block) {
291                            for (int i = newData.length; --i >= 0;) {
292                                    block[i] += newData[i];
293                            }
294                    }
295            }
296    
297            /**
298             * @see #fill(double[][], int, int)
299             */
300            public void fill(final double[][] data) {
301                    fill(data, 0, 0);
302            }
303    
304            /**
305             * Populate matrix with given data.
306             * 
307             * @param data
308             *            - to fill into matrix
309             * @param startRow
310             *            - row to start filling in data at
311             * @param startCol
312             *            - col to start at
313             */
314            public void fill(final double[][] data, final int startRow, final int startCol) {
315                    final int rows = data.length;
316                    final int cols = data[0].length;
317    
318                    assertTrue(startRow < rows && startRow < getRowCount(), "illegal startRow: %s", startRow);
319                    assertTrue(startCol < cols && startCol < getColumnCount(), "illegal startCol: %s", startCol);
320                    assertTrue(rows <= getRowCount(), "too many rows in input: %s: max allowed = %s", rows,
321                                    getRowCount());
322                    assertTrue(cols <= getColumnCount(), "too many columns in input: %s: max allowed = %s",
323                                    cols, getColumnCount());
324    
325                    for (int i = startRow; i < rows; i++) {
326                            for (int j = startCol; j < cols; j++) {
327                                    setDouble(data[i][j], i, j);
328                            }
329                    }
330            }
331    
332            /**
333             * Get block holding the specified row and column. If none exist, create
334             * one.
335             * 
336             * @param row
337             *            - in matrix
338             * @param column
339             *            - in matrix
340             * @return double[] block where the given row,column is held.
341             */
342            double[] getBlockData(int row, int column) {
343                    int blockNumber = layout.getBlockNumber(row, column);
344                    double[] block = data[blockNumber];
345    
346                    if (null == block) {
347                            block = new double[layout.getBlockSize(row, column)];
348                            data[blockNumber] = block;
349                    }
350                    return data[blockNumber];
351            }
352    
353            /**
354             * @return {@link BlockMatrixLayout} of this matrix.
355             */
356            public final BlockMatrixLayout getBlockLayout() {
357                    return this.layout;
358            }
359    
360            /** @return blockSize of this matrix. */
361            public final int getBlockStripeSize() {
362                    return layout.blockStripe;
363            }
364    
365            public double getDouble(final int row, final int col) {
366                    double[] block = data[layout.getBlockNumber(row, col)];
367                    if (null == block) {
368                            return 0.0d;
369                    }
370                    return block[layout.getIndexInBlock(row, col)];
371            }
372    
373            public double getDouble(long row, long column) {
374                    return this.getDouble((int) row, (int) column);
375            }
376    
377            public long[] getSize() {
378                    return this.size;
379            }
380    
381            /**
382             * Shortcut to create a BlockMatrix for target
383             */
384            public Matrix mtimes(Matrix m2) {
385                    if (m2 instanceof DenseDoubleMatrix2D) {
386                            final DenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D((int) getRowCount(),
387                                            (int) m2.getColumnCount(), layout.blockStripe, BlockOrder.ROWMAJOR);
388                            Mtimes.DENSEDOUBLEMATRIX2D.calc(this, (DenseDoubleMatrix2D) m2, result);
389                            return result;
390                    } else {
391                            return super.mtimes(m2);
392                    }
393            }
394    
395            public void setDouble(double value, int row, int column) {
396                    final double[] block = getBlockData(row, column);
397                    block[layout.getIndexInBlock(row, column)] = value;
398            }
399    
400            public void setDouble(double value, long row, long column) {
401                    setDouble(value, (int) row, (int) column);
402            }
403    
404            /**
405             * Returns the transpose of the current matrix.
406             * 
407             * @return transpose of this matrix.
408             */
409            @Override
410            public Matrix transpose() throws MatrixException {
411                    return transpose(Ret.NEW);
412            }
413    
414            @Override
415            public BlockDenseDoubleMatrix2D clone() {
416                    return new BlockDenseDoubleMatrix2D(this);
417            }
418    
419            @Override
420            public Matrix transpose(final Ret returnType) {
421                    // swap rows, column dimensions and block order
422                    final BlockOrder transOrder = (BlockOrder.ROWMAJOR == layout.blockOrder ? BlockOrder.COLUMNMAJOR
423                                    : BlockOrder.ROWMAJOR);
424    
425                    final int step = layout.blockStripe;
426                    final BlockDenseDoubleMatrix2D transMat = new BlockDenseDoubleMatrix2D(
427                                    (int) getColumnCount(), (int) getRowCount(), step, transOrder);
428    
429                    double[] block;
430    
431                    for (int i = 0; i < layout.rows; i += step) {
432                            for (int j = 0; j < layout.columns; j += step) {
433                                    // shuffle blocks to new position
434                                    final int blockNumberA = layout.getBlockNumber(i, j);
435                                    block = this.data[blockNumberA];
436    
437                                    if (returnType == Ret.NEW && null != block) {
438                                            // cannot use Arrays.copyOf(): not supported in Java 5
439                                            final double[] newBlock = new double[block.length];
440                                            System.arraycopy(block, 0, newBlock, 0, block.length);
441                                            block = newBlock;
442                                    }
443    
444                                    final int blockNumberB = transMat.layout.getBlockNumber(j, i);
445                                    transMat.data[blockNumberB] = block;
446                            }
447                    }
448    
449                    if (returnType == Ret.ORIG) {
450                            this.layout = transMat.layout;
451                            this.data = transMat.data;
452                            System.arraycopy(transMat.size, 0, this.size, 0, transMat.size.length);
453                            if (getAnnotation() != null) {
454                                    setAnnotation(Transpose.transposeAnnotation(getAnnotation(), Coordinates
455                                                    .transpose(getSize())));
456                            }
457                            return this;
458                    } else if (returnType == Ret.LINK) {
459                            return super.transpose(Ret.LINK);
460                    } else {
461                            if (getAnnotation() != null) {
462                                    transMat.setAnnotation(Transpose.transposeAnnotation(getAnnotation(), Coordinates
463                                                    .transpose(getSize())));
464                            }
465                            return transMat;
466                    }
467            }
468    
469            public final double[][] getBlockDoubleArray2D() {
470                    return data;
471            }
472    
473            public Matrix plus(final double value) {
474                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
475                    if (result.data.length < 100) {
476                            for (int i = result.data.length; --i != -1;) {
477                                    double[] block = result.data[i];
478                                    if (block == null) {
479                                            block = new double[layout.blockArea];
480                                    }
481                                    for (int j = block.length; --j != -1;) {
482                                            block[j] += value;
483                                    }
484                            }
485                    } else {
486                            new PFor(0, result.data.length) {
487                                    public void step(int i) {
488                                            double[] block = result.data[i];
489                                            if (block == null) {
490                                                    block = new double[layout.blockArea];
491                                            }
492                                            for (int j = block.length; --j != -1;) {
493                                                    block[j] += value;
494                                            }
495                                    }
496                            };
497                    }
498                    return result;
499            }
500    
501            public Matrix plus(Matrix value) {
502                    if (value instanceof BlockDenseDoubleMatrix2D) {
503                            final BlockDenseDoubleMatrix2D b = (BlockDenseDoubleMatrix2D) value;
504                            if (b.layout.rows == layout.rows && b.layout.columns == layout.columns
505                                            && b.layout.blockOrder == layout.blockOrder
506                                            && b.layout.blockStripe == layout.blockStripe) {
507                                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
508                                    if (result.data.length < 100) {
509                                            for (int i = result.data.length; --i != -1;) {
510                                                    final double[] block2 = b.data[i];
511                                                    if (block2 == null) {
512                                                            continue;
513                                                    }
514                                                    if (result.data[i] == null) {
515                                                            result.data[i] = new double[b.data[i].length];
516                                                    }
517                                                    final double[] block = result.data[i];
518                                                    for (int j = block.length; --j != -1;) {
519                                                            block[j] += block2[j];
520                                                    }
521                                            }
522                                    } else {
523                                            new PFor(0, result.data.length - 1) {
524                                                    public void step(int i) {
525                                                            final double[] block2 = b.data[i];
526                                                            if (block2 == null) {
527                                                                    return;
528                                                            }
529                                                            if (result.data[i] == null) {
530                                                                    result.data[i] = new double[b.data[i].length];
531                                                            }
532                                                            final double[] block = result.data[i];
533                                                            for (int j = block.length; --j != -1;) {
534                                                                    block[j] += block2[j];
535                                                            }
536                                                    }
537                                            };
538                                    }
539                                    return result;
540                            }
541                    }
542                    return super.plus(value);
543            }
544    
545            public Matrix minus(Matrix value) {
546                    if (value instanceof BlockDenseDoubleMatrix2D) {
547                            final BlockDenseDoubleMatrix2D b = (BlockDenseDoubleMatrix2D) value;
548                            if (b.layout.rows == layout.rows && b.layout.columns == layout.columns
549                                            && b.layout.blockOrder == layout.blockOrder
550                                            && b.layout.blockStripe == layout.blockStripe) {
551                                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
552                                    if (result.data.length < 100) {
553                                            for (int i = data.length; --i != -1;) {
554                                                    final double[] block2 = b.data[i];
555                                                    if (block2 == null) {
556                                                            continue;
557                                                    }
558                                                    if (result.data[i] == null) {
559                                                            result.data[i] = new double[b.data[i].length];
560                                                    }
561                                                    final double[] block = result.data[i];
562                                                    for (int j = block.length; --j != -1;) {
563                                                            block[j] -= block2[j];
564                                                    }
565                                            }
566                                    } else {
567                                            new PFor(0, result.data.length - 1) {
568                                                    public void step(int i) {
569                                                            final double[] block2 = b.data[i];
570                                                            if (block2 == null) {
571                                                                    return;
572                                                            }
573                                                            if (result.data[i] == null) {
574                                                                    result.data[i] = new double[b.data[i].length];
575                                                            }
576                                                            final double[] block = result.data[i];
577                                                            for (int j = block.length; --j != -1;) {
578                                                                    block[j] -= block2[j];
579                                                            }
580                                                    }
581                                            };
582                                    }
583                                    return result;
584                            }
585                    }
586                    return super.minus(value);
587            }
588    
589            public Matrix times(Matrix value) {
590                    if (value instanceof BlockDenseDoubleMatrix2D) {
591                            final BlockDenseDoubleMatrix2D b = (BlockDenseDoubleMatrix2D) value;
592                            if (b.layout.rows == layout.rows && b.layout.columns == layout.columns
593                                            && b.layout.blockOrder == layout.blockOrder
594                                            && b.layout.blockStripe == layout.blockStripe) {
595                                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
596                                    if (result.data.length < 100) {
597                                            for (int i = result.data.length; --i != -1;) {
598                                                    final double[] block2 = b.data[i];
599                                                    if (block2 == null) {
600                                                            // multiply with 0.0, clear block in result
601                                                            result.data[i] = null;
602                                                    } else {
603                                                            final double[] block = result.data[i];
604                                                            if (block == null) {
605                                                                    // multiply with 0.0
606                                                                    continue;
607                                                            }
608                                                            for (int j = block.length; --j != -1;) {
609                                                                    block[j] *= block2[j];
610                                                            }
611                                                    }
612                                            }
613                                    } else {
614                                            new PFor(0, result.data.length - 1) {
615                                                    public void step(int i) {
616                                                            final double[] block2 = b.data[i];
617                                                            if (block2 == null) {
618                                                                    // multiply with 0.0, clear block in result
619                                                                    result.data[i] = null;
620                                                            } else {
621                                                                    final double[] block = result.data[i];
622                                                                    if (block == null) {
623                                                                            // multiply with 0.0
624                                                                            return;
625                                                                    }
626                                                                    for (int j = block.length; --j != -1;) {
627                                                                            block[j] *= block2[j];
628                                                                    }
629                                                            }
630                                                    }
631                                            };
632                                    }
633                                    return result;
634                            }
635                    }
636                    return super.times(value);
637            }
638    
639            public Matrix divide(Matrix value) {
640                    if (value instanceof BlockDenseDoubleMatrix2D) {
641                            final BlockDenseDoubleMatrix2D b = (BlockDenseDoubleMatrix2D) value;
642                            if (b.layout.rows == layout.rows && b.layout.columns == layout.columns
643                                            && b.layout.blockOrder == layout.blockOrder
644                                            && b.layout.blockStripe == layout.blockStripe) {
645                                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
646                                    if (result.data.length < 100) {
647                                            for (int i = result.data.length; --i != -1;) {
648                                                    final double[] block2 = b.data[i];
649                                                    if (block2 == null) {
650                                                            // divide by 0.0, fill block with NaN
651                                                            if (result.data[i] == null) {
652                                                                    result.data[i] = new double[layout.blockArea];
653                                                            }
654                                                            Arrays.fill(result.data[i], Double.NaN);
655                                                    } else {
656                                                            final double[] block = result.data[i];
657                                                            if (block == null) {
658                                                                    // divide 0.0 by x = 0.0: nothing to do
659                                                                    continue;
660                                                            }
661                                                            for (int j = block.length; --j != -1;) {
662                                                                    block[j] /= block2[j];
663                                                            }
664                                                    }
665                                            }
666                                    } else {
667                                            new PFor(0, result.data.length - 1) {
668                                                    public void step(int i) {
669                                                            final double[] block2 = b.data[i];
670                                                            if (block2 == null) {
671                                                                    // divide by 0.0, fill block with NaN
672                                                                    if (result.data[i] == null) {
673                                                                            result.data[i] = new double[layout.blockArea];
674                                                                    }
675                                                                    Arrays.fill(result.data[i], Double.NaN);
676                                                            } else {
677                                                                    final double[] block = result.data[i];
678                                                                    if (block == null) {
679                                                                            // divide 0.0 by x = 0.0: nothing to do
680                                                                            return;
681                                                                    }
682                                                                    for (int j = block.length; --j != -1;) {
683                                                                            block[j] /= block2[j];
684                                                                    }
685                                                            }
686                                                    }
687                                            };
688                                    }
689    
690                                    return result;
691                            }
692                    }
693                    return super.times(value);
694            }
695    
696            public Matrix minus(final double value) {
697                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
698                    if (result.data.length < 100) {
699                            for (int i = result.data.length; --i != -1;) {
700                                    double[] block = result.data[i];
701                                    if (block == null) {
702                                            block = new double[layout.blockArea];
703                                    }
704                                    for (int j = block.length; --j != -1;) {
705                                            block[j] -= value;
706                                    }
707                            }
708                    } else {
709                            new PFor(0, result.data.length - 1) {
710                                    public void step(int i) {
711                                            double[] block = result.data[i];
712                                            if (block == null) {
713                                                    block = new double[layout.blockArea];
714                                            }
715                                            for (int j = block.length; --j != -1;) {
716                                                    block[j] -= value;
717                                            }
718                                    }
719                            };
720                    }
721                    return result;
722            }
723    
724            public Matrix times(final double value) {
725                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
726                    if (result.data.length < 100) {
727                            for (int i = result.data.length; --i != -1;) {
728                                    final double[] block = result.data[i];
729                                    if (block != null) {
730                                            for (int j = block.length; --j != -1;) {
731                                                    block[j] *= value;
732                                            }
733                                    }
734                            }
735                    } else {
736                            new PFor(0, result.data.length - 1) {
737                                    public void step(int i) {
738                                            final double[] block = result.data[i];
739                                            if (block != null) {
740                                                    for (int j = block.length; --j != -1;) {
741                                                            block[j] *= value;
742                                                    }
743                                            }
744                                    }
745                            };
746                    }
747                    return result;
748            }
749    
750            public Matrix divide(final double value) {
751                    final BlockDenseDoubleMatrix2D result = new BlockDenseDoubleMatrix2D(this);
752                    if (result.data.length < 100) {
753                            for (int i = result.data.length; --i != -1;) {
754                                    double[] block = result.data[i];
755                                    if (block == null) {
756                                            block = new double[layout.blockArea];
757                                    }
758                                    for (int j = block.length; --j != -1;) {
759                                            block[j] /= value;
760                                    }
761                            }
762                    } else {
763                            new PFor(0, result.data.length - 1) {
764                                    public void step(int i) {
765                                            double[] block = result.data[i];
766                                            if (block == null) {
767                                                    block = new double[layout.blockArea];
768                                            }
769                                            for (int j = block.length; --j != -1;) {
770                                                    block[j] /= value;
771                                            }
772                                    }
773                            };
774                    }
775                    return result;
776            }
777    
778            /**
779             * Change layout of blocks in this matrix (e.g. switch form rowmajor to
780             * columnmajor).
781             * 
782             * @param order
783             *            - new block layout order.
784             * @return old BlockOrder.
785             */
786            public BlockOrder setBlockOrder(BlockOrder order) {
787                    assertTrue(order != null, "block order cannot be null");
788    
789                    if (order == layout.blockOrder) {
790                            return order; // quick exit, already same
791                    }
792    
793                    BlockMatrixLayout newLayout = new BlockMatrixLayout(layout.rows, layout.columns,
794                                    layout.blockStripe, order);
795    
796                    // swap order of each block
797                    for (int r = 0; r < layout.rows; r += layout.blockStripe) {
798                            for (int c = 0; c < layout.columns; c += layout.blockStripe) {
799    
800                                    int blockNumber = layout.getBlockNumber(r, c);
801                                    if (data[blockNumber] == null) {
802                                            continue;
803                                    } else if (order == BlockOrder.ROWMAJOR) {
804                                            data[blockNumber] = layout.toRowMajorBlock(data[blockNumber], r, c);
805                                    } else {
806                                            data[blockNumber] = layout.toColMajorBlock(data[blockNumber], r, c);
807                                    }
808    
809                            }
810                    }
811    
812                    BlockOrder oldLayout = this.layout.blockOrder;
813                    this.layout = newLayout;
814                    return oldLayout;
815    
816            }
817    
818    };