Java tutorial
/** * (C) Copyright IBM Corp. 2010, 2015 * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * */ package com.ibm.bi.dml.runtime.matrix.data; import java.io.DataInput; import java.io.DataOutput; import java.io.Externalizable; import java.io.IOException; import java.io.ObjectInput; import java.io.ObjectInputStream; import java.io.ObjectOutput; import java.io.ObjectOutputStream; import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; import java.util.Iterator; import java.util.Map; import org.apache.commons.math3.random.Well1024a; import org.apache.hadoop.io.DataInputBuffer; import com.ibm.bi.dml.conf.ConfigurationManager; import com.ibm.bi.dml.conf.DMLConfig; import com.ibm.bi.dml.hops.Hop.OpOp2; import com.ibm.bi.dml.hops.OptimizerUtils; import com.ibm.bi.dml.lops.MMTSJ.MMTSJType; import com.ibm.bi.dml.lops.MapMultChain.ChainType; import com.ibm.bi.dml.lops.PartialAggregate.CorrectionLocationType; import com.ibm.bi.dml.parser.DMLTranslator; import com.ibm.bi.dml.runtime.DMLRuntimeException; import com.ibm.bi.dml.runtime.DMLUnsupportedOperationException; import com.ibm.bi.dml.runtime.functionobjects.Builtin; import com.ibm.bi.dml.runtime.functionobjects.CM; import com.ibm.bi.dml.runtime.functionobjects.CTable; import com.ibm.bi.dml.runtime.functionobjects.DiagIndex; import com.ibm.bi.dml.runtime.functionobjects.Divide; import com.ibm.bi.dml.runtime.functionobjects.KahanFunction; import com.ibm.bi.dml.runtime.functionobjects.KahanPlus; import com.ibm.bi.dml.runtime.functionobjects.KahanPlusSq; import com.ibm.bi.dml.runtime.functionobjects.Multiply; import com.ibm.bi.dml.runtime.functionobjects.Plus; import com.ibm.bi.dml.runtime.functionobjects.ReduceAll; import com.ibm.bi.dml.runtime.functionobjects.ReduceCol; import com.ibm.bi.dml.runtime.functionobjects.ReduceRow; import com.ibm.bi.dml.runtime.functionobjects.SortIndex; import com.ibm.bi.dml.runtime.functionobjects.SwapIndex; import com.ibm.bi.dml.runtime.instructions.cp.CM_COV_Object; import com.ibm.bi.dml.runtime.instructions.cp.DoubleObject; import com.ibm.bi.dml.runtime.instructions.cp.KahanObject; import com.ibm.bi.dml.runtime.instructions.cp.ScalarObject; import com.ibm.bi.dml.runtime.matrix.MatrixCharacteristics; import com.ibm.bi.dml.runtime.matrix.data.LibMatrixBincell.BinaryAccessType; import com.ibm.bi.dml.runtime.matrix.mapred.IndexedMatrixValue; import com.ibm.bi.dml.runtime.matrix.mapred.MRJobConfiguration; import com.ibm.bi.dml.runtime.matrix.operators.AggregateBinaryOperator; import com.ibm.bi.dml.runtime.matrix.operators.AggregateOperator; import com.ibm.bi.dml.runtime.matrix.operators.AggregateUnaryOperator; import com.ibm.bi.dml.runtime.matrix.operators.BinaryOperator; import com.ibm.bi.dml.runtime.matrix.operators.CMOperator; import com.ibm.bi.dml.runtime.matrix.operators.COVOperator; import com.ibm.bi.dml.runtime.matrix.operators.Operator; import com.ibm.bi.dml.runtime.matrix.operators.QuaternaryOperator; import com.ibm.bi.dml.runtime.matrix.operators.ReorgOperator; import com.ibm.bi.dml.runtime.matrix.operators.ScalarOperator; import com.ibm.bi.dml.runtime.matrix.operators.UnaryOperator; import com.ibm.bi.dml.runtime.util.DataConverter; import com.ibm.bi.dml.runtime.util.FastBufferedDataInputStream; import com.ibm.bi.dml.runtime.util.FastBufferedDataOutputStream; import com.ibm.bi.dml.runtime.util.IndexRange; import com.ibm.bi.dml.runtime.util.UtilFunctions; public class MatrixBlock extends MatrixValue implements Externalizable { private static final long serialVersionUID = 7319972089143154056L; //sparsity nnz threshold, based on practical experiments on space consumption and performance public static final double SPARSITY_TURN_POINT = 0.4; //sparsity threshold for ultra-sparse matrix operations (40nnz in a 1kx1k block) public static final double ULTRA_SPARSITY_TURN_POINT = 0.00004; //basic header (int rlen, int clen, byte type) public static final int HEADER_SIZE = 9; public enum BlockType { EMPTY_BLOCK, ULTRA_SPARSE_BLOCK, //ultra sparse representation, in-mem same as sparse SPARSE_BLOCK, //sparse representation, see sparseRows DENSE_BLOCK, //dense representation, see denseBlock } //matrix meta data protected int rlen = -1; protected int clen = -1; protected boolean sparse = true; protected long nonZeros = 0; //matrix data (sparse or dense) protected double[] denseBlock = null; protected SparseRow[] sparseRows = null; //sparse-block-specific attributes (allocation only) protected int estimatedNNzsPerRow = -1; //ctable-specific attributes protected int maxrow = -1; protected int maxcolumn = -1; //grpaggregate-specific attributes (optional) protected int numGroups = -1; //diag-specific attributes (optional) protected boolean diag = false; //////// // Matrix Constructors // public MatrixBlock() { rlen = 0; clen = 0; sparse = true; nonZeros = 0; maxrow = 0; maxcolumn = 0; } public MatrixBlock(int rl, int cl, boolean sp) { rlen = rl; clen = cl; sparse = sp; nonZeros = 0; maxrow = 0; maxcolumn = 0; } public MatrixBlock(int rl, int cl, long estnnzs) { this(rl, cl, false, estnnzs); // Adjust sparsity based on estimated non-zeros double denseSize = estimateSizeDenseInMemory(rl, cl); double sparseSize = estimateSizeSparseInMemory(rl, cl, (double) estnnzs / (rl * cl)); this.sparse = (denseSize < sparseSize ? false : true); } public MatrixBlock(int rl, int cl, boolean sp, long estnnzs) { this(rl, cl, sp); estimatedNNzsPerRow = (int) Math.ceil((double) estnnzs / (double) rl); } public MatrixBlock(MatrixBlock that) { this.copy(that); } //////// // Initialization methods // (reset, init, allocate, etc) public void reset() { reset(-rlen); } public void reset(long estnnzs) { estimatedNNzsPerRow = (int) Math.ceil((double) estnnzs / (double) rlen); if (sparse) { resetSparse(); } else { if (denseBlock != null) { if (denseBlock.length < rlen * clen) denseBlock = null; else Arrays.fill(denseBlock, 0, rlen * clen, 0); } } nonZeros = 0; //operation-specific attributes maxrow = rlen; maxcolumn = clen; numGroups = -1; } public void reset(int rl, int cl) { rlen = rl; clen = cl; nonZeros = 0; reset(); } public void reset(int rl, int cl, long estnnzs) { rlen = rl; clen = cl; nonZeros = 0; reset(estnnzs); } public void reset(int rl, int cl, boolean sp) { sparse = sp; reset(rl, cl); } public void reset(int rl, int cl, boolean sp, long estnnzs) { sparse = sp; reset(rl, cl, estnnzs); } public void resetSparse() { if (sparseRows != null) { for (int i = 0; i < Math.min(rlen, sparseRows.length); i++) if (sparseRows[i] != null) sparseRows[i].reset(estimatedNNzsPerRow, clen); } } public void resetDenseWithValue(int rl, int cl, double v) throws DMLRuntimeException { estimatedNNzsPerRow = -1; rlen = rl; clen = cl; sparse = false; if (v == 0) { reset(); return; } //allocate dense block allocateDenseBlock(); //init with constant value (non-zero, see above) int limit = rlen * clen; Arrays.fill(denseBlock, 0, limit, v); nonZeros = limit; } /** * NOTE: This method is designed only for dense representation. * * @param arr * @param r * @param c * @throws DMLRuntimeException */ public void init(double[][] arr, int r, int c) throws DMLRuntimeException { //input checks if (sparse) throw new DMLRuntimeException( "MatrixBlockDSM.init() can be invoked only on matrices with dense representation."); if (r * c > rlen * clen) throw new DMLRuntimeException("MatrixBlockDSM.init() invoked with too large dimensions (" + r + "," + c + ") vs (" + rlen + "," + clen + ")"); //allocate or resize dense block allocateDenseBlock(); //copy and compute nnz for (int i = 0, ix = 0; i < r; i++, ix += clen) System.arraycopy(arr[i], 0, denseBlock, ix, arr[i].length); recomputeNonZeros(); maxrow = r; maxcolumn = c; } /** * NOTE: This method is designed only for dense representation. * * @param arr * @param r * @param c * @throws DMLRuntimeException */ public void init(double[] arr, int r, int c) throws DMLRuntimeException { //input checks if (sparse) throw new DMLRuntimeException( "MatrixBlockDSM.init() can be invoked only on matrices with dense representation."); if (r * c > rlen * clen) throw new DMLRuntimeException("MatrixBlockDSM.init() invoked with too large dimensions (" + r + "," + c + ") vs (" + rlen + "," + clen + ")"); //allocate or resize dense block allocateDenseBlock(); //copy and compute nnz System.arraycopy(arr, 0, denseBlock, 0, arr.length); recomputeNonZeros(); maxrow = r; maxcolumn = c; } /** * * @param val * @param r * @param c * @throws DMLRuntimeException */ public void init(double val, int r, int c) throws DMLRuntimeException { //input checks if (sparse) throw new DMLRuntimeException( "MatrixBlockDSM.init() can be invoked only on matrices with dense representation."); if (r * c > rlen * clen) throw new DMLRuntimeException("MatrixBlockDSM.init() invoked with too large dimensions (" + r + "," + c + ") vs (" + rlen + "," + clen + ")"); if (val != 0) { //allocate or resize dense block allocateDenseBlock(); if (r * c == rlen * clen) { //FULL MATRIX INIT //memset value Arrays.fill(denseBlock, val); } else { //PARTIAL MATRIX INIT //rowwise memset value for (int i = 0, ix = 0; i < r; i++, ix += clen) Arrays.fill(denseBlock, ix, ix + c, val); } //set non zeros to input dims nonZeros = r * c; } maxrow = r; maxcolumn = c; } /** * * @return */ public boolean isAllocated() { if (sparse) return (sparseRows != null); else return (denseBlock != null); } /** * @throws DMLRuntimeException * */ public void allocateDenseBlock() throws RuntimeException { allocateDenseBlock(true); } /** * */ public void allocateDenseOrSparseBlock() { if (sparse) allocateSparseRowsBlock(); else allocateDenseBlock(); } /** * * @param clearNNZ * @throws DMLRuntimeException */ public void allocateDenseBlock(boolean clearNNZ) throws RuntimeException { long limit = (long) rlen * clen; //check max size constraint (16GB dense), since java arrays are limited to 2^(32-1) elements) if (limit > Integer.MAX_VALUE) { throw new RuntimeException("Dense in-memory matrix block (" + rlen + "x" + clen + ") exceeds supported size of " + Integer.MAX_VALUE + " elements (16GB). " + "Please, reduce the JVM heapsize to execute this in MR."); } //allocate block if non-existing or too small (guaranteed to be 0-initialized), if (denseBlock == null || denseBlock.length < limit) { denseBlock = new double[(int) limit]; } //clear nnz if necessary if (clearNNZ) { nonZeros = 0; } } /** * */ public void allocateSparseRowsBlock() { allocateSparseRowsBlock(true); } /** * * @param clearNNZ */ public void allocateSparseRowsBlock(boolean clearNNZ) { //allocate block if non-existing or too small (guaranteed to be 0-initialized), if (sparseRows == null) { sparseRows = new SparseRow[rlen]; } else if (sparseRows.length < rlen) { SparseRow[] oldSparseRows = sparseRows; sparseRows = new SparseRow[rlen]; for (int i = 0; i < Math.min(oldSparseRows.length, rlen); i++) { sparseRows[i] = oldSparseRows[i]; } } //clear nnz if necessary if (clearNNZ) { nonZeros = 0; } } /** * This should be called only in the read and write functions for CP * This function should be called before calling any setValueDenseUnsafe() * * @param rl * @param cl * @throws DMLRuntimeException */ public void allocateDenseBlockUnsafe(int rl, int cl) throws DMLRuntimeException { sparse = false; rlen = rl; clen = cl; //allocate dense block allocateDenseBlock(); } /** * Allows to cleanup all previously allocated sparserows or denseblocks. * This is for example required in reading a matrix with many empty blocks * via distributed cache into in-memory list of blocks - not cleaning blocks * from non-empty blocks would significantly increase the total memory consumption. * */ public void cleanupBlock(boolean dense, boolean sparse) { if (dense) denseBlock = null; if (sparse) sparseRows = null; } //////// // Metadata information public int getNumRows() { return rlen; } /** * NOTE: setNumRows() and setNumColumns() are used only in tertiaryInstruction (for contingency tables) * and pmm for meta corrections. * * @param _r */ public void setNumRows(int r) { rlen = r; } public int getNumColumns() { return clen; } public void setNumColumns(int c) { clen = c; } public long getNonZeros() { return nonZeros; } public void setNonZeros(long nnz) { nonZeros = nnz; } public boolean isVector() { return (rlen == 1 || clen == 1); } /** * Return the maximum row encountered WITHIN the current block * */ public int getMaxRow() { if (!sparse) return getNumRows(); else return maxrow; } public void setMaxRow(int r) { maxrow = r; } /** * Return the maximum column encountered WITHIN the current block * */ public int getMaxColumn() { if (!sparse) return getNumColumns(); else return maxcolumn; } public void setMaxColumn(int c) { maxcolumn = c; } @Override public boolean isEmpty() { return isEmptyBlock(false); } public boolean isEmptyBlock() { return isEmptyBlock(true); } public boolean isEmptyBlock(boolean safe) { boolean ret = false; if (sparse && sparseRows == null) ret = true; else if (!sparse && denseBlock == null) ret = true; if (nonZeros == 0) { //prevent under-estimation if (safe) recomputeNonZeros(); ret = (nonZeros == 0); } return ret; } public void setDiag() { diag = true; } public boolean isDiag() { return diag; } //////// // Data handling public double[] getDenseArray() { if (sparse) return null; return denseBlock; } public SparseRow[] getSparseRows() { if (!sparse) return null; return sparseRows; } public SparseRowsIterator getSparseRowsIterator() { //check for valid format, should have been checked from outside if (!sparse) throw new RuntimeException("getSparseCellInterator should not be called for dense format"); return new SparseRowsIterator(rlen, sparseRows); } public SparseRowsIterator getSparseRowsIterator(int rl, int ru) { //check for valid format, should have been checked from outside if (!sparse) throw new RuntimeException("getSparseCellInterator should not be called for dense format"); return new SparseRowsIterator(rl, ru, sparseRows); } @Override public void getCellValues(Collection<Double> ret) { int limit = rlen * clen; if (sparse) { if (sparseRows == null) { for (int i = 0; i < limit; i++) ret.add(0.0); } else { for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; double[] container = sparseRows[r].getValueContainer(); for (int j = 0; j < sparseRows[r].size(); j++) ret.add(container[j]); } int zeros = limit - ret.size(); for (int i = 0; i < zeros; i++) ret.add(0.0); } } else { if (denseBlock == null) { for (int i = 0; i < limit; i++) ret.add(0.0); } else { for (int i = 0; i < limit; i++) ret.add(denseBlock[i]); } } } @Override public void getCellValues(Map<Double, Integer> ret) { int limit = rlen * clen; if (sparse) { if (sparseRows == null) { ret.put(0.0, limit); } else { for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; double[] container = sparseRows[r].getValueContainer(); for (int j = 0; j < sparseRows[r].size(); j++) { Double v = container[j]; Integer old = ret.get(v); if (old != null) ret.put(v, old + 1); else ret.put(v, 1); } } int zeros = limit - ret.size(); Integer old = ret.get(0.0); if (old != null) ret.put(0.0, old + zeros); else ret.put(0.0, zeros); } } else { if (denseBlock == null) { ret.put(0.0, limit); } else { for (int i = 0; i < limit; i++) { double v = denseBlock[i]; Integer old = ret.get(v); if (old != null) ret.put(v, old + 1); else ret.put(v, 1); } } } } @Override public double getValue(int r, int c) { if (r > rlen || c > clen) throw new RuntimeException("indexes (" + r + "," + c + ") out of range (" + rlen + "," + clen + ")"); if (sparse) { if (sparseRows == null || sparseRows.length <= r || sparseRows[r] == null) return 0; return sparseRows[r].get(c); } else { if (denseBlock == null) return 0; return denseBlock[r * clen + c]; } } @Override public void setValue(int r, int c, double v) { if (r > rlen || c > clen) throw new RuntimeException("indexes (" + r + "," + c + ") out of range (" + rlen + "," + clen + ")"); if (sparse) { if ((sparseRows == null || sparseRows.length <= r || sparseRows[r] == null) && v == 0.0) return; //allocation on demand allocateSparseRowsBlock(false); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(estimatedNNzsPerRow, clen); if (sparseRows[r].set(c, v)) nonZeros++; } else { if (denseBlock == null && v == 0.0) return; //allocate and init dense block (w/o overwriting nnz) allocateDenseBlock(false); int index = r * clen + c; if (denseBlock[index] == 0) nonZeros++; denseBlock[index] = v; if (v == 0) nonZeros--; } } @Override public void setValue(CellIndex index, double v) { setValue(index.row, index.column, v); } @Override /** * If (r,c) \in Block, add v to existing cell * If not, add a new cell with index (r,c). * * This function intentionally avoids the maintenance of NNZ for efficiency. * */ public void addValue(int r, int c, double v) { if (sparse) { //allocation on demand allocateSparseRowsBlock(false); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(estimatedNNzsPerRow, clen); double curV = sparseRows[r].get(c); curV += v; sparseRows[r].set(c, curV); } else { //allocate and init dense block (w/o overwriting nnz) allocateDenseBlock(false); int index = r * clen + c; denseBlock[index] += v; } } /** * * @param r * @param c * @return */ public double quickGetValue(int r, int c) { if (sparse) { if (sparseRows == null || sparseRows.length <= r || sparseRows[r] == null) return 0; return sparseRows[r].get(c); } else { if (denseBlock == null) return 0; return denseBlock[r * clen + c]; } } /** * * @param r * @param c * @param v */ public void quickSetValue(int r, int c, double v) { if (sparse) { //early abort if ((sparseRows == null || sparseRows.length <= r || sparseRows[r] == null) && v == 0) return; //allocation on demand allocateSparseRowsBlock(false); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(estimatedNNzsPerRow, clen); //set value and maintain nnz if (sparseRows[r].set(c, v)) nonZeros += (v != 0) ? 1 : -1; } else { //early abort if (denseBlock == null && v == 0) return; //allocate and init dense block (w/o overwriting nnz) allocateDenseBlock(false); //set value and maintain nnz int index = r * clen + c; if (denseBlock[index] == 0) nonZeros++; denseBlock[index] = v; if (v == 0) nonZeros--; } } public double getValueDenseUnsafe(int r, int c) { if (denseBlock == null) return 0; return denseBlock[r * clen + c]; } /** * This can be only called when you know you have properly allocated spaces for a dense representation * and r and c are in the the range of the dimension * Note: this function won't keep track of the nozeros */ public void setValueDenseUnsafe(int r, int c, double v) { denseBlock[r * clen + c] = v; } public double getValueSparseUnsafe(int r, int c) { if (sparseRows == null || sparseRows.length <= r || sparseRows[r] == null) return 0; return sparseRows[r].get(c); } /** * Append value is only used when values are appended at the end of each row for the sparse representation * This can only be called, when the caller knows the access pattern of the block * * @param r * @param c * @param v */ public void appendValue(int r, int c, double v) { //early abort (append guarantees no overwrite) if (v == 0) return; if (!sparse) //DENSE { //allocate on demand (w/o overwriting nnz) allocateDenseBlock(false); //set value and maintain nnz denseBlock[r * clen + c] = v; nonZeros++; } else //SPARSE { //allocation on demand (w/o overwriting nnz) allocateSparseRowsBlock(false); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(estimatedNNzsPerRow, clen); //set value and maintain nnz sparseRows[r].append(c, v); nonZeros++; } } public void appendRow(int r, SparseRow values) { if (values == null) return; if (sparse) { //allocation on demand allocateSparseRowsBlock(false); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(values); else sparseRows[r].copy(values); nonZeros += values.size(); } else { int[] cols = values.getIndexContainer(); double[] vals = values.getValueContainer(); for (int i = 0; i < values.size(); i++) quickSetValue(r, cols[i], vals[i]); } } /** * * @param that * @param rowoffset * @param coloffset */ public void appendToSparse(MatrixBlock that, int rowoffset, int coloffset) { if (that == null || that.isEmptyBlock(false)) return; //nothing to append //init sparse rows if necessary allocateSparseRowsBlock(false); if (that.sparse) //SPARSE <- SPARSE { for (int i = 0; i < that.rlen; i++) { SparseRow brow = that.sparseRows[i]; if (brow != null && !brow.isEmpty()) { int aix = rowoffset + i; int len = brow.size(); int[] ix = brow.getIndexContainer(); double[] val = brow.getValueContainer(); if (sparseRows[aix] == null) sparseRows[aix] = new SparseRow(estimatedNNzsPerRow, clen); for (int j = 0; j < len; j++) sparseRows[aix].append(coloffset + ix[j], val[j]); } } } else //SPARSE <- DENSE { for (int i = 0; i < that.rlen; i++) { int aix = rowoffset + i; for (int j = 0, bix = i * that.clen; j < that.clen; j++) { double val = that.denseBlock[bix + j]; if (val != 0) { if (sparseRows[aix] == null)//create sparserow only if required sparseRows[aix] = new SparseRow(estimatedNNzsPerRow, clen); sparseRows[aix].append(coloffset + j, val); } } } } } /** * */ public void sortSparseRows() { if (!sparse || sparseRows == null) return; for (SparseRow arow : sparseRows) if (arow != null && arow.size() > 1) arow.sort(); } /** * Utility function for computing the min non-zero value. * * @return * @throws DMLRuntimeException */ public double minNonZero() throws DMLRuntimeException { //check for empty block and return immediately if (isEmptyBlock()) return -1; //NOTE: usually this method is only applied on dense vectors and hence not really tuned yet. double min = Double.MAX_VALUE; for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double val = quickGetValue(i, j); if (val != 0) min = Math.min(min, val); } return min; } /** * Wrapper method for reduceall-min of a matrix. * * @return * @throws DMLRuntimeException */ public double min() throws DMLRuntimeException { //construct operator AggregateOperator aop = new AggregateOperator(Double.MAX_VALUE, Builtin.getBuiltinFnObject("min")); AggregateUnaryOperator auop = new AggregateUnaryOperator(aop, ReduceAll.getReduceAllFnObject()); //execute operation MatrixBlock out = new MatrixBlock(1, 1, false); LibMatrixAgg.aggregateUnaryMatrix(this, out, auop); return out.quickGetValue(0, 0); } /** * Wrapper method for reduceall-max of a matrix. * * @return * @throws DMLRuntimeException */ public double max() throws DMLRuntimeException { //construct operator AggregateOperator aop = new AggregateOperator(-Double.MAX_VALUE, Builtin.getBuiltinFnObject("max")); AggregateUnaryOperator auop = new AggregateUnaryOperator(aop, ReduceAll.getReduceAllFnObject()); //execute operation MatrixBlock out = new MatrixBlock(1, 1, false); LibMatrixAgg.aggregateUnaryMatrix(this, out, auop); return out.quickGetValue(0, 0); } /** * Wrapper method for reduceall-sum of a matrix. * * @return Sum of the values in the matrix. * @throws DMLRuntimeException */ public double sum() throws DMLRuntimeException { KahanPlus kplus = KahanPlus.getKahanPlusFnObject(); return sumWithFn(kplus); } /** * Wrapper method for reduceall-sumSq of a matrix. * * @return Sum of the squared values in the matrix. * @throws DMLRuntimeException */ public double sumSq() throws DMLRuntimeException { KahanPlusSq kplusSq = KahanPlusSq.getKahanPlusSqFnObject(); return sumWithFn(kplusSq); } /** * Wrapper method for reduceall-sum of a matrix using the given * Kahan function for summation. * * @param kfunc A Kahan function object to use for summation. * @return Sum of the values in the matrix with the given * function applied. * @throws DMLRuntimeException */ private double sumWithFn(KahanFunction kfunc) throws DMLRuntimeException { //construct operator CorrectionLocationType corrLoc = CorrectionLocationType.LASTCOLUMN; ReduceAll reduceAllObj = ReduceAll.getReduceAllFnObject(); AggregateOperator aop = new AggregateOperator(0, kfunc, true, corrLoc); AggregateUnaryOperator auop = new AggregateUnaryOperator(aop, reduceAllObj); //execute operation MatrixBlock out = new MatrixBlock(1, 2, false); LibMatrixAgg.aggregateUnaryMatrix(this, out, auop); return out.quickGetValue(0, 0); } //////// // sparsity handling functions /** * Returns the current representation (true for sparse). * */ public boolean isInSparseFormat() { return sparse; } /** * * @return */ public boolean isUltraSparse() { double sp = ((double) nonZeros / rlen) / clen; //check for sparse representation in order to account for vectors in dense return sparse && sp < ULTRA_SPARSITY_TURN_POINT && nonZeros < 40; } /** * Evaluates if this matrix block should be in sparse format in * memory. Note that this call does not change the representation - * for this please call examSparsity. * * @return */ public boolean evalSparseFormatInMemory() { long lrlen = (long) rlen; long lclen = (long) clen; long lnonZeros = (long) nonZeros; //ensure exact size estimates for write if (lnonZeros <= 0) { recomputeNonZeros(); lnonZeros = (long) nonZeros; } //decide on in-memory representation return evalSparseFormatInMemory(lrlen, lclen, lnonZeros); } @SuppressWarnings("unused") private boolean evalSparseFormatInMemory(boolean transpose) { int lrlen = (transpose) ? clen : rlen; int lclen = (transpose) ? rlen : clen; long lnonZeros = (long) nonZeros; //ensure exact size estimates for write if (lnonZeros <= 0) { recomputeNonZeros(); lnonZeros = (long) nonZeros; } //decide on in-memory representation return evalSparseFormatInMemory(lrlen, lclen, lnonZeros); } /** * Evaluates if this matrix block should be in sparse format on * disk. This applies to any serialized matrix representation, i.e., * when writing to in-memory buffer pool pages or writing to local fs * or hdfs. * * @return */ public boolean evalSparseFormatOnDisk() { long lrlen = (long) rlen; long lclen = (long) clen; //ensure exact size estimates for write if (nonZeros <= 0) { recomputeNonZeros(); } //decide on in-memory representation return evalSparseFormatOnDisk(lrlen, lclen, nonZeros); } /** * Evaluates if this matrix block should be in sparse format in * memory. Depending on the current representation, the state of the * matrix block is changed to the right representation if necessary. * Note that this consumes for the time of execution memory for both * representations. * * @throws DMLRuntimeException */ public void examSparsity() throws DMLRuntimeException { //determine target representation boolean sparseDst = evalSparseFormatInMemory(); //check for empty blocks (e.g., sparse-sparse) if (isEmptyBlock(false)) cleanupBlock(true, true); //change representation if required (also done for //empty blocks in order to set representation flags) if (sparse && !sparseDst) sparseToDense(); else if (!sparse && sparseDst) denseToSparse(); } /** * Evaluates if a matrix block with the given characteristics should be in sparse format * in memory. * * @param brlen * @param bclen * @param nnz * @return */ public static boolean evalSparseFormatInMemory(final long nrows, final long ncols, final long nnz) { //evaluate sparsity threshold double lsparsity = (double) nnz / nrows / ncols; boolean lsparse = (lsparsity < SPARSITY_TURN_POINT); //compare size of sparse and dense representation in order to prevent //that the sparse size exceed the dense size since we use the dense size //as worst-case estimate if unknown (and it requires less io from //main memory). double sizeSparse = estimateSizeSparseInMemory(nrows, ncols, lsparsity); double sizeDense = estimateSizeDenseInMemory(nrows, ncols); return lsparse && (sizeSparse < sizeDense); } /** * Evaluates if a matrix block with the given characteristics should be in sparse format * on disk (or in any other serialized representation). * * @param brlen * @param bclen * @param nnz * @return */ public static boolean evalSparseFormatOnDisk(final long nrows, final long ncols, final long nnz) { //evaluate sparsity threshold double lsparsity = ((double) nnz / nrows) / ncols; boolean lsparse = (lsparsity < SPARSITY_TURN_POINT); double sizeUltraSparse = estimateSizeUltraSparseOnDisk(nrows, ncols, nnz); double sizeSparse = estimateSizeSparseOnDisk(nrows, ncols, nnz); double sizeDense = estimateSizeDenseOnDisk(nrows, ncols); return lsparse && (sizeSparse < sizeDense || sizeUltraSparse < sizeDense); } //////// // basic block handling functions /** * */ private void denseToSparse() { //set target representation sparse = true; //early abort on empty blocks if (denseBlock == null) return; //allocate sparse target block (reset required to maintain nnz again) allocateSparseRowsBlock(); reset(); //copy dense to sparse double[] a = denseBlock; SparseRow[] c = sparseRows; for (int i = 0, aix = 0; i < rlen; i++) for (int j = 0; j < clen; j++, aix++) if (a[aix] != 0) { if (c[i] == null) //create sparse row only if required c[i] = new SparseRow(estimatedNNzsPerRow, clen); c[i].append(j, a[aix]); nonZeros++; } //cleanup dense block denseBlock = null; } /** * * @throws DMLRuntimeException */ private void sparseToDense() throws DMLRuntimeException { //set target representation sparse = false; //early abort on empty blocks if (sparseRows == null) return; int limit = rlen * clen; if (limit < 0) { throw new DMLRuntimeException( "Unexpected error in sparseToDense().. limit < 0: " + rlen + ", " + clen + ", " + limit); } //allocate dense target block, but keep nnz (no need to maintain) allocateDenseBlock(false); Arrays.fill(denseBlock, 0, limit, 0); //copy sparse to dense SparseRow[] a = sparseRows; double[] c = denseBlock; for (int i = 0, cix = 0; i < rlen; i++, cix += clen) if (a[i] != null && !a[i].isEmpty()) { int alen = a[i].size(); int[] aix = a[i].getIndexContainer(); double[] avals = a[i].getValueContainer(); for (int j = 0; j < alen; j++) if (avals[j] != 0) c[cix + aix[j]] = avals[j]; } //cleanup sparse rows sparseRows = null; } public void recomputeNonZeros() { nonZeros = 0; if (sparse && sparseRows != null) { int limit = Math.min(rlen, sparseRows.length); for (int i = 0; i < limit; i++) if (sparseRows[i] != null) nonZeros += sparseRows[i].size(); } else if (!sparse && denseBlock != null) { int limit = rlen * clen; for (int i = 0; i < limit; i++) { //HotSpot JVM bug causes crash in presence of NaNs //nonZeros += (denseBlock[i]!=0) ? 1 : 0; if (denseBlock[i] != 0) nonZeros++; } } } protected long recomputeNonZeros(int rl, int ru, int cl, int cu) { long nnz = 0; if (sparse) { if (sparseRows != null) { int rlimit = Math.min(ru + 1, Math.min(rlen, sparseRows.length)); if (cl == 0 && cu == clen - 1) //specific case: all cols { for (int i = rl; i < rlimit; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) nnz += sparseRows[i].size(); } else if (cl == cu) //specific case: one column { for (int i = rl; i < rlimit; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) nnz += (sparseRows[i].get(cl) != 0) ? 1 : 0; } else //general case { int astart, aend; for (int i = rl; i < rlimit; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) { SparseRow arow = sparseRows[i]; astart = arow.searchIndexesFirstGTE(cl); aend = arow.searchIndexesFirstGTE(cu); nnz += (astart != -1) ? (aend - astart + 1) : 0; } } } } else { if (denseBlock != null) { for (int i = rl, ix = rl * clen; i <= ru; i++, ix += clen) for (int j = cl; j <= cu; j++) { //HotSpot JVM bug causes crash in presence of NaNs //nnz += (denseBlock[ix+j]!=0) ? 1 : 0; if (denseBlock[ix + j] != 0) nnz++; } } } return nnz; } /** * Basic debugging primitive to check correctness of nnz, this method is not intended for * production use. * * @throws DMLRuntimeException */ public void checkNonZeros() throws DMLRuntimeException { //take non-zeros before and after recompute nnz long nnzBefore = getNonZeros(); recomputeNonZeros(); long nnzAfter = getNonZeros(); //raise exception if non-zeros don't match up if (nnzBefore != nnzAfter) throw new DMLRuntimeException("Number of non zeros incorrect: " + nnzBefore + " vs " + nnzAfter); } @Override public void copy(MatrixValue thatValue) { MatrixBlock that = checkType(thatValue); //copy into automatically determined representation copy(that, that.evalSparseFormatInMemory()); } @Override public void copy(MatrixValue thatValue, boolean sp) { MatrixBlock that = checkType(thatValue); if (this == that) //prevent data loss (e.g., on sparse-dense conversion) throw new RuntimeException("Copy must not overwrite itself!"); this.rlen = that.rlen; this.clen = that.clen; this.sparse = sp; estimatedNNzsPerRow = (int) Math.ceil((double) thatValue.getNonZeros() / (double) rlen); if (this.sparse && that.sparse) copySparseToSparse(that); else if (this.sparse && !that.sparse) copyDenseToSparse(that); else if (!this.sparse && that.sparse) copySparseToDense(that); else copyDenseToDense(that); } private void copySparseToSparse(MatrixBlock that) { this.nonZeros = that.nonZeros; if (that.isEmptyBlock(false)) { resetSparse(); return; } allocateSparseRowsBlock(false); for (int i = 0; i < Math.min(that.sparseRows.length, rlen); i++) { if (that.sparseRows[i] != null) { if (sparseRows[i] == null) sparseRows[i] = new SparseRow(that.sparseRows[i]); else sparseRows[i].copy(that.sparseRows[i]); } else if (this.sparseRows[i] != null) this.sparseRows[i].reset(estimatedNNzsPerRow, clen); } } private void copyDenseToDense(MatrixBlock that) { nonZeros = that.nonZeros; int limit = rlen * clen; //plain reset to 0 for empty input if (that.isEmptyBlock(false)) { if (denseBlock != null) Arrays.fill(denseBlock, 0, limit, 0); return; } //allocate and init dense block (w/o overwriting nnz) allocateDenseBlock(false); //actual copy System.arraycopy(that.denseBlock, 0, denseBlock, 0, limit); } private void copySparseToDense(MatrixBlock that) { this.nonZeros = that.nonZeros; if (that.isEmptyBlock(false)) { if (denseBlock != null) Arrays.fill(denseBlock, 0); return; } //allocate and init dense block (w/o overwriting nnz) allocateDenseBlock(false); int start = 0; for (int r = 0; r < Math.min(that.sparseRows.length, rlen); r++, start += clen) { if (that.sparseRows[r] == null) continue; double[] values = that.sparseRows[r].getValueContainer(); int[] cols = that.sparseRows[r].getIndexContainer(); for (int i = 0; i < that.sparseRows[r].size(); i++) { denseBlock[start + cols[i]] = values[i]; } } } private void copyDenseToSparse(MatrixBlock that) { nonZeros = that.nonZeros; if (that.isEmptyBlock(false)) { resetSparse(); return; } allocateSparseRowsBlock(false); for (int i = 0, ix = 0; i < rlen; i++) { if (sparseRows[i] != null) sparseRows[i].reset(estimatedNNzsPerRow, clen); for (int j = 0; j < clen; j++) { double val = that.denseBlock[ix++]; if (val != 0) { if (sparseRows[i] == null) //create sparse row only if required sparseRows[i] = new SparseRow(estimatedNNzsPerRow, clen); sparseRows[i].append(j, val); } } } } /** * In-place copy of matrix src into the index range of the existing current matrix. * Note that removal of existing nnz in the index range and nnz maintenance is * only done if 'awareDestNZ=true', * * @param rl * @param ru * @param cl * @param cu * @param src * @param awareDestNZ * true, forces (1) to remove existing non-zeros in the index range of the * destination if not present in src and (2) to internally maintain nnz * false, assume empty index range in destination and do not maintain nnz * (the invoker is responsible to recompute nnz after all copies are done) * @throws DMLRuntimeException */ public void copy(int rl, int ru, int cl, int cu, MatrixBlock src, boolean awareDestNZ) throws DMLRuntimeException { if (sparse && src.sparse) copySparseToSparse(rl, ru, cl, cu, src, awareDestNZ); else if (sparse && !src.sparse) copyDenseToSparse(rl, ru, cl, cu, src, awareDestNZ); else if (!sparse && src.sparse) copySparseToDense(rl, ru, cl, cu, src, awareDestNZ); else copyDenseToDense(rl, ru, cl, cu, src, awareDestNZ); } private void copySparseToSparse(int rl, int ru, int cl, int cu, MatrixBlock src, boolean awareDestNZ) { //handle empty src and dest if (src.isEmptyBlock(false)) { if (awareDestNZ && sparseRows != null) copyEmptyToSparse(rl, ru, cl, cu, true); return; } if (sparseRows == null) sparseRows = new SparseRow[rlen]; else if (awareDestNZ) { copyEmptyToSparse(rl, ru, cl, cu, true); //explicit clear if awareDestNZ because more efficient since //src will have multiple columns and only few overwriting values } //copy values int alen; int[] aix; double[] avals; for (int i = 0; i < src.rlen; i++) { SparseRow arow = src.sparseRows[i]; if (arow != null && !arow.isEmpty()) { alen = arow.size(); aix = arow.getIndexContainer(); avals = arow.getValueContainer(); if (sparseRows[rl + i] == null || sparseRows[rl + i].isEmpty()) { sparseRows[rl + i] = new SparseRow(estimatedNNzsPerRow, clen); SparseRow brow = sparseRows[rl + i]; for (int j = 0; j < alen; j++) brow.append(cl + aix[j], avals[j]); if (awareDestNZ) nonZeros += brow.size(); } else if (awareDestNZ) //general case (w/ awareness NNZ) { SparseRow brow = sparseRows[rl + i]; int lnnz = brow.size(); if (cl == cu && cl == aix[0]) { if (avals[0] == 0) brow.delete(cl); else brow.set(cl, avals[0]); } else { brow.deleteIndexRange(cl, cu); for (int j = 0; j < alen; j++) brow.set(cl + aix[j], avals[j]); } nonZeros += (brow.size() - lnnz); } else //general case (w/o awareness NNZ) { SparseRow brow = sparseRows[rl + i]; //brow.set(cl, arow); for (int j = 0; j < alen; j++) brow.set(cl + aix[j], avals[j]); } } } } private void copySparseToDense(int rl, int ru, int cl, int cu, MatrixBlock src, boolean awareDestNZ) throws DMLRuntimeException { //handle empty src and dest if (src.isEmptyBlock(false)) { if (awareDestNZ && denseBlock != null) { nonZeros -= recomputeNonZeros(rl, ru, cl, cu); copyEmptyToDense(rl, ru, cl, cu); } return; } if (denseBlock == null) allocateDenseBlock(); else if (awareDestNZ) { nonZeros -= recomputeNonZeros(rl, ru, cl, cu); copyEmptyToDense(rl, ru, cl, cu); } //copy values int alen; int[] aix; double[] avals; for (int i = 0, ix = rl * clen; i < src.rlen; i++, ix += clen) { SparseRow arow = src.sparseRows[i]; if (arow != null && !arow.isEmpty()) { alen = arow.size(); aix = arow.getIndexContainer(); avals = arow.getValueContainer(); for (int j = 0; j < alen; j++) denseBlock[ix + cl + aix[j]] = avals[j]; if (awareDestNZ) nonZeros += alen; } } } private void copyDenseToSparse(int rl, int ru, int cl, int cu, MatrixBlock src, boolean awareDestNZ) { //handle empty src and dest if (src.isEmptyBlock(false)) { if (awareDestNZ && sparseRows != null) copyEmptyToSparse(rl, ru, cl, cu, true); return; } if (sparseRows == null) sparseRows = new SparseRow[rlen]; //no need to clear for awareDestNZ since overwritten //copy values double val; for (int i = 0, ix = 0; i < src.rlen; i++, ix += src.clen) { int rix = rl + i; if (sparseRows[rix] == null || sparseRows[rix].isEmpty()) { for (int j = 0; j < src.clen; j++) if ((val = src.denseBlock[ix + j]) != 0) { if (sparseRows[rix] == null) sparseRows[rix] = new SparseRow(estimatedNNzsPerRow, clen); sparseRows[rix].append(cl + j, val); } if (awareDestNZ && sparseRows[rix] != null) nonZeros += sparseRows[rix].size(); } else if (awareDestNZ) //general case (w/ awareness NNZ) { SparseRow brow = sparseRows[rix]; int lnnz = brow.size(); if (cl == cu) { if ((val = src.denseBlock[ix]) == 0) brow.delete(cl); else brow.set(cl, val); } else { brow.setIndexRange(cl, cu, src.denseBlock, ix, src.clen); } nonZeros += (brow.size() - lnnz); } else //general case (w/o awareness NNZ) { SparseRow brow = sparseRows[rix]; for (int j = 0; j < src.clen; j++) if ((val = src.denseBlock[ix + j]) != 0) brow.set(cl + j, val); } } } private void copyDenseToDense(int rl, int ru, int cl, int cu, MatrixBlock src, boolean awareDestNZ) throws DMLRuntimeException { //handle empty src and dest if (src.isEmptyBlock(false)) { if (awareDestNZ && denseBlock != null) { nonZeros -= recomputeNonZeros(rl, ru, cl, cu); copyEmptyToDense(rl, ru, cl, cu); } return; } if (denseBlock == null) allocateDenseBlock(); //no need to clear for awareDestNZ since overwritten if (awareDestNZ) nonZeros = nonZeros - recomputeNonZeros(rl, ru, cl, cu) + src.nonZeros; //copy values int rowLen = cu - cl + 1; if (clen == src.clen) //optimization for equal width System.arraycopy(src.denseBlock, 0, denseBlock, rl * clen + cl, src.rlen * src.clen); else for (int i = 0, ix1 = 0, ix2 = rl * clen + cl; i < src.rlen; i++, ix1 += src.clen, ix2 += clen) { System.arraycopy(src.denseBlock, ix1, denseBlock, ix2, rowLen); } } private void copyEmptyToSparse(int rl, int ru, int cl, int cu, boolean updateNNZ) { if (cl == cu) //specific case: column vector { if (updateNNZ) { for (int i = rl; i <= ru; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) { int lnnz = sparseRows[i].size(); sparseRows[i].delete(cl); nonZeros += (sparseRows[i].size() - lnnz); } } else { for (int i = rl; i <= ru; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) sparseRows[i].delete(cl); } } else { if (updateNNZ) { for (int i = rl; i <= ru; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) { int lnnz = sparseRows[i].size(); sparseRows[i].deleteIndexRange(cl, cu); nonZeros += (sparseRows[i].size() - lnnz); } } else { for (int i = rl; i <= ru; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) sparseRows[i].deleteIndexRange(cl, cu); } } } private void copyEmptyToDense(int rl, int ru, int cl, int cu) { int rowLen = cu - cl + 1; if (clen == rowLen) //optimization for equal width Arrays.fill(denseBlock, rl * clen + cl, ru * clen + cu + 1, 0); else for (int i = rl, ix2 = rl * clen + cl; i <= ru; i++, ix2 += clen) Arrays.fill(denseBlock, ix2, ix2 + rowLen, 0); } /** * Merge disjoint: merges all non-zero values of the given input into the current * matrix block. Note that this method does NOT check for overlapping entries; * it's the callers reponsibility of ensuring disjoint matrix blocks. * * The appendOnly parameter is only relevant for sparse target blocks; if true, * we only append values and do not sort sparse rows for each call; this is useful * whenever we merge iterators of matrix blocks into one target block. * * @param that * @param appendOnly * @throws DMLRuntimeException */ public void merge(MatrixBlock that, boolean appendOnly) throws DMLRuntimeException { //check for empty input source (nothing to merge) if (that == null || that.isEmptyBlock(false)) return; //check dimensions (before potentially copy to prevent implicit dimension change) //this also does a best effort check for disjoint input blocks via the number of non-zeros if (rlen != that.rlen || clen != that.clen) throw new DMLRuntimeException("Dimension mismatch on merge disjoint (target=" + rlen + "x" + clen + ", source=" + that.rlen + "x" + that.clen + ")"); if ((long) nonZeros + that.nonZeros > (long) rlen * clen) throw new DMLRuntimeException("Number of non-zeros mismatch on merge disjoint (target=" + rlen + "x" + clen + ", nnz target=" + nonZeros + ", nnz source=" + that.nonZeros + ")"); //check for empty target (copy in full) if (isEmptyBlock(false)) { copy(that); return; } //core matrix block merge (guaranteed non-empty source/target, nnz maintenance not required) long nnz = nonZeros + that.nonZeros; if (sparse) mergeIntoSparse(that, appendOnly); else mergeIntoDense(that); //maintain number of nonzeros nonZeros = nnz; } /** * * @param that */ private void mergeIntoDense(MatrixBlock that) { if (that.sparse) //DENSE <- SPARSE { double[] a = denseBlock; SparseRow[] b = that.sparseRows; int m = rlen; int n = clen; for (int i = 0, aix = 0; i < m; i++, aix += n) if (b[i] != null && !b[i].isEmpty()) { SparseRow brow = b[i]; int blen = brow.size(); int[] bix = brow.getIndexContainer(); double[] bval = brow.getValueContainer(); for (int j = 0; j < blen; j++) if (bval[j] != 0) a[aix + bix[j]] = bval[j]; } } else //DENSE <- DENSE { double[] a = denseBlock; double[] b = that.denseBlock; int len = rlen * clen; for (int i = 0; i < len; i++) a[i] = (b[i] != 0) ? b[i] : a[i]; } } /** * * @param that * @param appendOnly */ private void mergeIntoSparse(MatrixBlock that, boolean appendOnly) { if (that.sparse) //SPARSE <- SPARSE { SparseRow[] a = sparseRows; SparseRow[] b = that.sparseRows; int m = rlen; for (int i = 0; i < m; i++) { if (b[i] != null && !b[i].isEmpty()) { if (a[i] == null || a[i].isEmpty()) { //copy entire sparse row (no sort required) a[i] = new SparseRow(b[i]); } else { boolean appended = false; SparseRow arow = a[i]; SparseRow brow = b[i]; int blen = brow.size(); int[] bix = brow.getIndexContainer(); double[] bval = brow.getValueContainer(); for (int j = 0; j < blen; j++) { if (bval[j] != 0) { arow.append(bix[j], bval[j]); appended = true; } } //only sort if value appended if (!appendOnly && appended) arow.sort(); } } } } else //SPARSE <- DENSE { SparseRow[] a = sparseRows; double[] b = that.denseBlock; int m = rlen; int n = clen; for (int i = 0, bix = 0; i < m; i++, bix += n) { boolean appended = false; for (int j = 0; j < n; j++) { if (b[bix + j] != 0) { appendValue(i, j, b[bix + j]); appended = true; } } //only sort if value appended if (!appendOnly && appended) a[i].sort(); } } } //////// // Input/Output functions @Override public void readFields(DataInput in) throws IOException { //read basic header (int rlen, int clen, byte type) rlen = in.readInt(); clen = in.readInt(); byte bformat = in.readByte(); //check type information if (bformat < 0 || bformat >= BlockType.values().length) throw new IOException( "invalid format: '" + bformat + "' (need to be 0-" + BlockType.values().length + ")."); BlockType format = BlockType.values()[bformat]; try { switch (format) { case ULTRA_SPARSE_BLOCK: nonZeros = readNnzInfo(in, true); sparse = evalSparseFormatInMemory(rlen, clen, nonZeros); cleanupBlock(true, true); //clean all if (sparse) readUltraSparseBlock(in); else readUltraSparseToDense(in); break; case SPARSE_BLOCK: nonZeros = readNnzInfo(in, false); sparse = evalSparseFormatInMemory(rlen, clen, nonZeros); cleanupBlock(sparse, !sparse); if (sparse) readSparseBlock(in); else readSparseToDense(in); break; case DENSE_BLOCK: sparse = false; cleanupBlock(false, true); //reuse dense readDenseBlock(in); //always dense in-mem if dense on disk break; case EMPTY_BLOCK: sparse = true; cleanupBlock(true, true); //clean all nonZeros = 0; break; } } catch (DMLRuntimeException ex) { throw new IOException("Error reading block of type '" + format.toString() + "'.", ex); } } /** * * @param in * @throws IOException * @throws DMLRuntimeException */ private void readDenseBlock(DataInput in) throws IOException, DMLRuntimeException { allocateDenseBlock(true); //allocate block, clear nnz int limit = rlen * clen; if (in instanceof MatrixBlockDataInput) //fast deserialize { MatrixBlockDataInput mbin = (MatrixBlockDataInput) in; nonZeros = mbin.readDoubleArray(limit, denseBlock); } else if (in instanceof DataInputBuffer && MRJobConfiguration.USE_BINARYBLOCK_SERIALIZATION) { //workaround because sequencefile.reader.next(key, value) does not yet support serialization framework DataInputBuffer din = (DataInputBuffer) in; MatrixBlockDataInput mbin = new FastBufferedDataInputStream(din); nonZeros = mbin.readDoubleArray(limit, denseBlock); ((FastBufferedDataInputStream) mbin).close(); } else //default deserialize { for (int i = 0; i < limit; i++) { denseBlock[i] = in.readDouble(); if (denseBlock[i] != 0) nonZeros++; } } } /** * * @param in * @throws IOException */ private void readSparseBlock(DataInput in) throws IOException { allocateSparseRowsBlock(false); resetSparse(); //reset all sparse rows if (in instanceof MatrixBlockDataInput) //fast deserialize { MatrixBlockDataInput mbin = (MatrixBlockDataInput) in; nonZeros = mbin.readSparseRows(rlen, sparseRows); } else if (in instanceof DataInputBuffer && MRJobConfiguration.USE_BINARYBLOCK_SERIALIZATION) { //workaround because sequencefile.reader.next(key, value) does not yet support serialization framework DataInputBuffer din = (DataInputBuffer) in; MatrixBlockDataInput mbin = new FastBufferedDataInputStream(din); nonZeros = mbin.readSparseRows(rlen, sparseRows); ((FastBufferedDataInputStream) mbin).close(); } else //default deserialize { for (int r = 0; r < rlen; r++) { int nr = in.readInt(); if (nr == 0) { if (sparseRows[r] != null) sparseRows[r].reset(estimatedNNzsPerRow, clen); continue; } if (sparseRows[r] == null) sparseRows[r] = new SparseRow(nr); else sparseRows[r].reset(nr, clen); for (int j = 0; j < nr; j++) sparseRows[r].append(in.readInt(), in.readDouble()); } } } /** * * @param in * @throws IOException * @throws DMLRuntimeException */ private void readSparseToDense(DataInput in) throws IOException, DMLRuntimeException { allocateDenseBlock(false); //allocate block Arrays.fill(denseBlock, 0); for (int r = 0; r < rlen; r++) { int nr = in.readInt(); for (int j = 0; j < nr; j++) { int c = in.readInt(); double val = in.readDouble(); denseBlock[r * clen + c] = val; } } } /** * * @param in * @throws IOException */ private void readUltraSparseBlock(DataInput in) throws IOException { allocateSparseRowsBlock(false); //adjust to size resetSparse(); //reset all sparse rows if (clen > 1) //ULTRA-SPARSE BLOCK { //block: read ijv-triples for (long i = 0; i < nonZeros; i++) { int r = in.readInt(); int c = in.readInt(); double val = in.readDouble(); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(1, clen); sparseRows[r].append(c, val); } } else //ULTRA-SPARSE COL { //col: read iv-pairs (should never happen since always dense) for (long i = 0; i < nonZeros; i++) { int r = in.readInt(); double val = in.readDouble(); if (sparseRows[r] == null) sparseRows[r] = new SparseRow(1, 1); sparseRows[r].append(0, val); } } } /** * * @param in * @throws IOException * @throws DMLRuntimeException */ private void readUltraSparseToDense(DataInput in) throws IOException, DMLRuntimeException { allocateDenseBlock(false); //allocate block Arrays.fill(denseBlock, 0); if (clen > 1) //ULTRA-SPARSE BLOCK { //block: read ijv-triples for (long i = 0; i < nonZeros; i++) { int r = in.readInt(); int c = in.readInt(); double val = in.readDouble(); denseBlock[r * clen + c] = val; } } else //ULTRA-SPARSE COL { //col: read iv-pairs for (long i = 0; i < nonZeros; i++) { int r = in.readInt(); double val = in.readDouble(); denseBlock[r] = val; } } } @Override public void write(DataOutput out) throws IOException { //determine format boolean sparseSrc = sparse; boolean sparseDst = evalSparseFormatOnDisk(); //write first part of header out.writeInt(rlen); out.writeInt(clen); if (sparseSrc) { //write sparse to * if (sparseRows == null || nonZeros == 0) writeEmptyBlock(out); else if (nonZeros < rlen && sparseDst) writeSparseToUltraSparse(out); else if (sparseDst) writeSparseBlock(out); else writeSparseToDense(out); } else { //write dense to * if (denseBlock == null || nonZeros == 0) writeEmptyBlock(out); else if (nonZeros < rlen && sparseDst) writeDenseToUltraSparse(out); else if (sparseDst) writeDenseToSparse(out); else writeDenseBlock(out); } } /** * * @param out * @throws IOException */ private void writeEmptyBlock(DataOutput out) throws IOException { //empty blocks do not need to materialize row information out.writeByte(BlockType.EMPTY_BLOCK.ordinal()); } /** * * @param out * @throws IOException */ private void writeDenseBlock(DataOutput out) throws IOException { out.writeByte(BlockType.DENSE_BLOCK.ordinal()); int limit = rlen * clen; if (out instanceof MatrixBlockDataOutput) //fast serialize ((MatrixBlockDataOutput) out).writeDoubleArray(limit, denseBlock); else //general case (if fast serialize not supported) for (int i = 0; i < limit; i++) out.writeDouble(denseBlock[i]); } /** * * @param out * @throws IOException */ private void writeSparseBlock(DataOutput out) throws IOException { out.writeByte(BlockType.SPARSE_BLOCK.ordinal()); writeNnzInfo(out, false); if (out instanceof MatrixBlockDataOutput) //fast serialize ((MatrixBlockDataOutput) out).writeSparseRows(rlen, sparseRows); else //general case (if fast serialize not supported) { int r = 0; for (; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) out.writeInt(0); else { int nr = sparseRows[r].size(); out.writeInt(nr); int[] cols = sparseRows[r].getIndexContainer(); double[] values = sparseRows[r].getValueContainer(); for (int j = 0; j < nr; j++) { out.writeInt(cols[j]); out.writeDouble(values[j]); } } } for (; r < rlen; r++) out.writeInt(0); } } /** * * @param out * @throws IOException */ private void writeSparseToUltraSparse(DataOutput out) throws IOException { out.writeByte(BlockType.ULTRA_SPARSE_BLOCK.ordinal()); writeNnzInfo(out, true); long wnnz = 0; if (clen > 1) //ULTRA-SPARSE BLOCK { //block: write ijv-triples for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) if (sparseRows[r] != null && !sparseRows[r].isEmpty()) { int alen = sparseRows[r].size(); int[] aix = sparseRows[r].getIndexContainer(); double[] avals = sparseRows[r].getValueContainer(); for (int j = 0; j < alen; j++) { //ultra-sparse block: write ijv-triples out.writeInt(r); out.writeInt(aix[j]); out.writeDouble(avals[j]); wnnz++; } } } else //ULTRA-SPARSE COL { //block: write iv-pairs (should never happen since always dense) for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) if (sparseRows[r] != null && !sparseRows[r].isEmpty()) { out.writeInt(r); out.writeDouble(sparseRows[r].getValueContainer()[0]); wnnz++; } } //validity check (nnz must exactly match written nnz) if (nonZeros != wnnz) { throw new IOException( "Invalid number of serialized non-zeros: " + wnnz + " (expected: " + nonZeros + ")"); } } /** * * @param out * @throws IOException */ private void writeSparseToDense(DataOutput out) throws IOException { //write block type 'dense' out.writeByte(BlockType.DENSE_BLOCK.ordinal()); //write data (from sparse to dense) if (sparseRows == null) //empty block for (int i = 0; i < rlen * clen; i++) out.writeDouble(0); else //existing sparse block { for (int i = 0; i < rlen; i++) { if (i < sparseRows.length && sparseRows[i] != null && !sparseRows[i].isEmpty()) { SparseRow arow = sparseRows[i]; int alen = arow.size(); int[] aix = arow.getIndexContainer(); double[] avals = arow.getValueContainer(); //foreach non-zero value, fill with 0s if required for (int j = 0, j2 = 0; j2 < alen; j++, j2++) { for (; j < aix[j2]; j++) out.writeDouble(0); out.writeDouble(avals[j2]); } //remaining 0 values in row for (int j = aix[alen - 1] + 1; j < clen; j++) out.writeDouble(0); } else //empty row for (int j = 0; j < clen; j++) out.writeDouble(0); } } } /** * * @param out * @throws IOException */ private void writeDenseToUltraSparse(DataOutput out) throws IOException { out.writeByte(BlockType.ULTRA_SPARSE_BLOCK.ordinal()); writeNnzInfo(out, true); long wnnz = 0; if (clen > 1) //ULTRA-SPARSE BLOCK { //block: write ijv-triples for (int r = 0, ix = 0; r < rlen; r++) for (int c = 0; c < clen; c++, ix++) if (denseBlock[ix] != 0) { out.writeInt(r); out.writeInt(c); out.writeDouble(denseBlock[ix]); wnnz++; } } else //ULTRA-SPARSE COL { //col: write iv-pairs for (int r = 0; r < rlen; r++) if (denseBlock[r] != 0) { out.writeInt(r); out.writeDouble(denseBlock[r]); wnnz++; } } //validity check (nnz must exactly match written nnz) if (nonZeros != wnnz) { throw new IOException( "Invalid number of serialized non-zeros: " + wnnz + " (expected: " + nonZeros + ")"); } } /** * * @param out * @throws IOException */ private void writeDenseToSparse(DataOutput out) throws IOException { out.writeByte(BlockType.SPARSE_BLOCK.ordinal()); //block type writeNnzInfo(out, false); int start = 0; for (int r = 0; r < rlen; r++) { //count nonzeros int nr = 0; for (int i = start; i < start + clen; i++) if (denseBlock[i] != 0.0) nr++; out.writeInt(nr); for (int c = 0; c < clen; c++) { if (denseBlock[start] != 0.0) { out.writeInt(c); out.writeDouble(denseBlock[start]); } start++; } } } /** * * @param in * @throws IOException */ private long readNnzInfo(DataInput in, boolean ultrasparse) throws IOException { //note: if ultrasparse, int always sufficient because nnz<rlen // where rlen is limited to integer long lrlen = (long) rlen; long lclen = (long) clen; //read long if required, otherwise int (see writeNnzInfo, consistency required) if (lrlen * lclen > Integer.MAX_VALUE && !ultrasparse) { nonZeros = in.readLong(); } else { nonZeros = in.readInt(); } return nonZeros; } /** * * @param out * @throws IOException */ private void writeNnzInfo(DataOutput out, boolean ultrasparse) throws IOException { //note: if ultrasparse, int always sufficient because nnz<rlen // where rlen is limited to integer long lrlen = (long) rlen; long lclen = (long) clen; //write long if required, otherwise int if (lrlen * lclen > Integer.MAX_VALUE && !ultrasparse) { out.writeLong(nonZeros); } else { out.writeInt((int) nonZeros); } } /** * Redirects the default java serialization via externalizable to our default * hadoop writable serialization for efficient broadcast/rdd deserialization. * * @param is * @throws IOException */ public void readExternal(ObjectInput is) throws IOException { if (is instanceof ObjectInputStream) { //fast deserialize of dense/sparse blocks ObjectInputStream ois = (ObjectInputStream) is; FastBufferedDataInputStream fis = new FastBufferedDataInputStream(ois); readFields(fis); } else { //default deserialize (general case) readFields(is); } } /** * Redirects the default java serialization via externalizable to our default * hadoop writable serialization for efficient broadcast/rdd serialization. * * @param is * @throws IOException */ public void writeExternal(ObjectOutput os) throws IOException { if (os instanceof ObjectOutputStream) { //fast serialize of dense/sparse blocks ObjectOutputStream oos = (ObjectOutputStream) os; FastBufferedDataOutputStream fos = new FastBufferedDataOutputStream(oos); write(fos); fos.flush(); } else { //default serialize (general case) write(os); } } /** * NOTE: The used estimates must be kept consistent with the respective write functions. * * @return */ public long getExactSizeOnDisk() { //determine format boolean sparseSrc = sparse; boolean sparseDst = evalSparseFormatOnDisk(); long lrlen = (long) rlen; long lclen = (long) clen; long lnonZeros = (long) nonZeros; //ensure exact size estimates for write if (lnonZeros <= 0) { recomputeNonZeros(); lnonZeros = (long) nonZeros; } //get exact size estimate (see write for the corresponding meaning) if (sparseSrc) { //write sparse to * if (sparseRows == null || lnonZeros == 0) return HEADER_SIZE; //empty block else if (lnonZeros < lrlen && sparseDst) return estimateSizeUltraSparseOnDisk(lrlen, lclen, lnonZeros); //ultra sparse block else if (sparseDst) return estimateSizeSparseOnDisk(lrlen, lclen, lnonZeros); //sparse block else return estimateSizeDenseOnDisk(lrlen, lclen); //dense block } else { //write dense to * if (denseBlock == null || lnonZeros == 0) return HEADER_SIZE; //empty block else if (lnonZeros < lrlen && sparseDst) return estimateSizeUltraSparseOnDisk(lrlen, lclen, lnonZeros); //ultra sparse block else if (sparseDst) return estimateSizeSparseOnDisk(lrlen, lclen, lnonZeros); //sparse block else return estimateSizeDenseOnDisk(lrlen, lclen); //dense block } } //////// // Estimates size and sparsity /** * Estimate size based on current sparse/dense representation. * * @return */ public long getSizeInMemory() { double sp = OptimizerUtils.getSparsity(rlen, clen, nonZeros); if (sparse) return estimateSizeSparseInMemory(rlen, clen, sp); else return estimateSizeDenseInMemory(rlen, clen); } /** * * @return */ public long estimateSizeInMemory() { double sp = OptimizerUtils.getSparsity(rlen, clen, nonZeros); return estimateSizeInMemory(rlen, clen, sp); } /** * * @param nrows * @param ncols * @param sparsity * @return */ public static long estimateSizeInMemory(long nrows, long ncols, double sparsity) { //determine sparse/dense representation boolean sparse = evalSparseFormatInMemory(nrows, ncols, (long) (sparsity * nrows * ncols)); //estimate memory consumption for sparse/dense if (sparse) return estimateSizeSparseInMemory(nrows, ncols, sparsity); else return estimateSizeDenseInMemory(nrows, ncols); } /** * * @param nrows * @param ncols * @return */ public static long estimateSizeDenseInMemory(long nrows, long ncols) { // basic variables and references sizes double size = 44; // core dense matrix block (double array) size += 8d * nrows * ncols; // robustness for long overflows return (long) Math.min(size, Long.MAX_VALUE); } /** * * @param nrows * @param ncols * @param sparsity * @return */ public static long estimateSizeSparseInMemory(long nrows, long ncols, double sparsity) { // basic variables and references sizes double size = 44; //NOTES: // * Each sparse row has a fixed overhead of 8B (reference) + 32B (object) + // 12B (3 int members), 32B (overhead int array), 32B (overhead double array), // * Each non-zero value requires 12B for the column-index/value pair. // * Overheads for arrays, objects, and references refer to 64bit JVMs // * If nnz < than rows we have only also empty rows. // account for sparsity and initial capacity double cnnz = Math.max(SparseRow.initialCapacity, Math.ceil(sparsity * ncols)); double rlen = Math.min(nrows, Math.ceil(sparsity * nrows * ncols)); size += rlen * (116 + 12 * cnnz); //sparse row size += nrows * 8d; //empty rows // robustness for long overflows return (long) Math.min(size, Long.MAX_VALUE); } /** * * @return */ public long estimateSizeOnDisk() { return estimateSizeOnDisk(rlen, clen, nonZeros); } /** * * @param nrows * @param ncols * @param nnz * @return */ public static long estimateSizeOnDisk(long nrows, long ncols, long nnz) { //determine sparse/dense representation boolean sparse = evalSparseFormatOnDisk(nrows, ncols, nnz); //estimate memory consumption for sparse/dense if (sparse && nnz < nrows) return estimateSizeUltraSparseOnDisk(nrows, ncols, nnz); else if (sparse) return estimateSizeSparseOnDisk(nrows, ncols, nnz); else return estimateSizeDenseOnDisk(nrows, ncols); } /** * * @param nrows * @param ncols * @return */ private static long estimateSizeDenseOnDisk(long nrows, long ncols) { //basic header (int rlen, int clen, byte type) long size = HEADER_SIZE; //data (all cells double) size += nrows * ncols * 8; return size; } /** * * @param nrows * @param ncols * @param nnz * @return */ private static long estimateSizeSparseOnDisk(long nrows, long ncols, long nnz) { //basic header: (int rlen, int clen, byte type) long size = HEADER_SIZE; //extended header (long nnz) size += (nrows * ncols > Integer.MAX_VALUE) ? 8 : 4; //data: (int num per row, int-double pair per non-zero value) size += nrows * 4 + nnz * 12; return size; } /** * * @param nrows * @param ncols * @param nnz * @return */ private static long estimateSizeUltraSparseOnDisk(long nrows, long ncols, long nnz) { //basic header (int rlen, int clen, byte type) long size = HEADER_SIZE; //extended header (int nnz, guaranteed by rlen<nnz) size += 4; //data (int-int-double triples per non-zero value) if (ncols > 1) //block: ijv-triples size += nnz * 16; else //column: iv-pairs size += nnz * 12; return size; } /** * * @param m1 * @param m2 * @param op * @return */ public static SparsityEstimate estimateSparsityOnAggBinary(MatrixBlock m1, MatrixBlock m2, AggregateBinaryOperator op) { //Since MatrixMultLib always uses a dense output (except for ultra-sparse mm) //with subsequent check for sparsity, we should always return a dense estimate. //Once, we support more aggregate binary operations, we need to change this. //WARNING: KEEP CONSISTENT WITH LIBMATRIXMULT //Note that it is crucial to report the right output representation because //in case of block reuse (e.g., mmcj) the output 'reset' refers to either //dense or sparse representation and hence would produce incorrect results //if we report the wrong representation (i.e., missing reset on ultrasparse mm). boolean ultrasparse = (m1.isUltraSparse() || m2.isUltraSparse()); return new SparsityEstimate(ultrasparse, m1.getNumRows() * m2.getNumRows()); } /** * * @param m1 * @param m2 * @param op * @return */ private static SparsityEstimate estimateSparsityOnBinary(MatrixBlock m1, MatrixBlock m2, BinaryOperator op) { SparsityEstimate est = new SparsityEstimate(); //estimate dense output for all sparse-unsafe operations, except DIV (because it commonly behaves like //sparse-safe but is not due to 0/0->NaN, this is consistent with the current hop sparsity estimate) //see also, special sparse-safe case for DIV in LibMatrixBincell if (!op.sparseSafe && !(op.fn instanceof Divide)) { est.sparse = false; return est; } BinaryAccessType atype = LibMatrixBincell.getBinaryAccessType(m1, m2); boolean outer = (atype == BinaryAccessType.OUTER_VECTOR_VECTOR); long m = m1.getNumRows(); long n = outer ? m2.getNumColumns() : m1.getNumColumns(); long nz1 = m1.getNonZeros(); long nz2 = m2.getNonZeros(); //account for matrix vector and vector/vector long estnnz = 0; if (atype == BinaryAccessType.OUTER_VECTOR_VECTOR) { //for outer vector operations the sparsity estimate is exactly known estnnz = nz1 * nz2; } else //DEFAULT CASE { if (atype == BinaryAccessType.MATRIX_COL_VECTOR) nz2 = nz2 * n; else if (atype == BinaryAccessType.MATRIX_ROW_VECTOR) nz2 = nz2 * m; //compute output sparsity consistent w/ the hop compiler OpOp2 bop = op.getBinaryOperatorOpOp2(); double sp1 = OptimizerUtils.getSparsity(m, n, nz1); double sp2 = OptimizerUtils.getSparsity(m, n, nz2); double spout = OptimizerUtils.getBinaryOpSparsity(sp1, sp2, bop, true); estnnz = UtilFunctions.toLong(spout * m * n); } est.sparse = evalSparseFormatInMemory(m, n, estnnz); est.estimatedNonZeros = estnnz; return est; } private boolean estimateSparsityOnSlice(int selectRlen, int selectClen, int finalRlen, int finalClen) { long ennz = (long) ((double) nonZeros / rlen / clen * selectRlen * selectClen); return evalSparseFormatInMemory(finalRlen, finalClen, ennz); } private boolean estimateSparsityOnLeftIndexing(long rlenm1, long clenm1, long nnzm1, long rlenm2, long clenm2, long nnzm2) { //min bound: nnzm1 - rlenm2*clenm2 + nnzm2 //max bound: min(rlenm1*rlenm2, nnzm1+nnzm2) long ennz = Math.min(rlenm1 * clenm1, nnzm1 + nnzm2); return evalSparseFormatInMemory(rlenm1, clenm1, ennz); } private boolean estimateSparsityOnGroupedAgg(long rlen, long groups) { long ennz = Math.min(groups, rlen); return evalSparseFormatInMemory(groups, 1, ennz); } //////// // Core block operations (called from instructions) /** * */ @Override public MatrixValue scalarOperations(ScalarOperator op, MatrixValue result) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock ret = checkType(result); // estimate the sparsity structure of result matrix boolean sp = this.sparse; // by default, we guess result.sparsity=input.sparsity if (!op.sparseSafe) sp = false; // if the operation is not sparse safe, then result will be in dense format //allocate the output matrix block if (ret == null) ret = new MatrixBlock(rlen, clen, sp, this.nonZeros); else ret.reset(rlen, clen, sp, this.nonZeros); //core scalar operations LibMatrixBincell.bincellOp(this, ret, op); return ret; } /** * */ @Override public MatrixValue unaryOperations(UnaryOperator op, MatrixValue result) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock ret = checkType(result); // estimate the sparsity structure of result matrix boolean sp = this.sparse; // by default, we guess result.sparsity=input.sparsity if (!op.sparseSafe) sp = false; // if the operation is not sparse safe, then result will be in dense format //allocate output if (ret == null) ret = new MatrixBlock(rlen, clen, sp, this.nonZeros); else ret.reset(rlen, clen, sp); //core execute if (LibMatrixAgg.isSupportedUnaryOperator(op)) { //e.g., cumsum/cumprod/cummin/cumax LibMatrixAgg.aggregateUnaryMatrix(this, ret, op); } else { //default execute unary operations if (op.sparseSafe) sparseUnaryOperations(op, ret); else denseUnaryOperations(op, ret); } //ensure empty results sparse representation //(no additional memory requirements) if (ret.isEmptyBlock(false)) ret.examSparsity(); return ret; } /** * * @param op * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ private void sparseUnaryOperations(UnaryOperator op, MatrixBlock ret) throws DMLUnsupportedOperationException, DMLRuntimeException { //early abort possible since sparse-safe if (isEmptyBlock(false)) return; final int m = rlen; final int n = clen; if (sparse) //SPARSE <- SPARSE { SparseRow[] a = sparseRows; for (int i = 0; i < m; i++) { if (a[i] != null && !a[i].isEmpty()) { int alen = a[i].size(); int[] aix = a[i].getIndexContainer(); double[] avals = a[i].getValueContainer(); for (int j = 0; j < alen; j++) { double val = op.fn.execute(avals[j]); ret.appendValue(i, aix[j], val); } } } } else //DENSE <- DENSE { //allocate dense output block ret.allocateDenseBlock(); double[] a = denseBlock; double[] c = ret.denseBlock; //unary op, incl nnz maintenance for (int i = 0, ix = 0; i < m; i++) { for (int j = 0; j < n; j++, ix++) { c[ix] = op.fn.execute(a[ix]); if (c[ix] != 0) ret.nonZeros++; } } } } /** * * @param op * @param ret * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ private void denseUnaryOperations(UnaryOperator op, MatrixBlock ret) throws DMLUnsupportedOperationException, DMLRuntimeException { //prepare 0-value init (determine if unnecessarily sparse-unsafe) double val0 = op.fn.execute(0); final int m = rlen; final int n = clen; //early abort possible if unnecessarily sparse unsafe //(otherwise full init with val0, no need for computation) if (isEmptyBlock(false)) { if (val0 != 0) ret.init(val0, m, n); return; } //redirection to sparse safe operation w/ init by val0 if (sparse && val0 != 0) ret.init(val0, m, n); sparseUnaryOperations(op, ret); } /** * * @param op * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ @Override public void unaryOperationsInPlace(UnaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { if (op.sparseSafe) sparseUnaryOperationsInPlace(op); else denseUnaryOperationsInPlace(op); } /** * only apply to non zero cells * * @param op * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ private void sparseUnaryOperationsInPlace(UnaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { //early abort possible since sparse-safe if (isEmptyBlock(false)) return; if (sparse) { nonZeros = 0; for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; double[] values = sparseRows[r].getValueContainer(); int[] cols = sparseRows[r].getIndexContainer(); int pos = 0; for (int i = 0; i < sparseRows[r].size(); i++) { double v = op.fn.execute(values[i]); if (v != 0) { values[pos] = v; cols[pos] = cols[i]; pos++; nonZeros++; } } sparseRows[r].truncate(pos); } } else { int limit = rlen * clen; nonZeros = 0; for (int i = 0; i < limit; i++) { denseBlock[i] = op.fn.execute(denseBlock[i]); if (denseBlock[i] != 0) nonZeros++; } } } private void denseUnaryOperationsInPlace(UnaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { if (sparse) //SPARSE MATRIX { double v; for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { v = op.fn.execute(quickGetValue(r, c)); quickSetValue(r, c, v); } } else//DENSE MATRIX { //early abort not possible because not sparsesafe if (denseBlock == null) allocateDenseBlock(); //compute values in-place and update nnz final int limit = rlen * clen; int lnnz = 0; for (int i = 0; i < limit; i++) { denseBlock[i] = op.fn.execute(denseBlock[i]); if (denseBlock[i] != 0) lnnz++; } nonZeros = lnnz; //IBM JVM bug (JDK6) causes crash for certain inputs (w/ infinities) //nonZeros = 0; //for(int i=0; i<limit; i++) //{ // denseBlock[i]=op.fn.execute(denseBlock[i]); // if(denseBlock[i]!=0) // nonZeros++; //} } } /** * */ public MatrixValue binaryOperations(BinaryOperator op, MatrixValue thatValue, MatrixValue result) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock that = checkType(thatValue); MatrixBlock ret = checkType(result); if (!LibMatrixBincell.isValidDimensionsBinary(this, that)) { throw new RuntimeException("Block sizes are not matched for binary " + "cell operations: " + this.rlen + "x" + this.clen + " vs " + that.rlen + "x" + that.clen); } //compute output dimensions boolean outer = (LibMatrixBincell.getBinaryAccessType(this, that) == BinaryAccessType.OUTER_VECTOR_VECTOR); int rows = rlen; int cols = outer ? that.clen : clen; //estimate output sparsity SparsityEstimate resultSparse = estimateSparsityOnBinary(this, that, op); if (ret == null) ret = new MatrixBlock(rows, cols, resultSparse.sparse, resultSparse.estimatedNonZeros); else ret.reset(rows, cols, resultSparse.sparse, resultSparse.estimatedNonZeros); //core binary cell operation LibMatrixBincell.bincellOp(this, that, ret, op); return ret; } /** * */ public void binaryOperationsInPlace(BinaryOperator op, MatrixValue thatValue) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock that = checkType(thatValue); if (!LibMatrixBincell.isValidDimensionsBinary(this, that)) { throw new RuntimeException("block sizes are not matched for binary " + "cell operations: " + this.rlen + "*" + this.clen + " vs " + that.rlen + "*" + that.clen); } //estimate output sparsity SparsityEstimate resultSparse = estimateSparsityOnBinary(this, that, op); if (resultSparse.sparse && !this.sparse) denseToSparse(); else if (!resultSparse.sparse && this.sparse) sparseToDense(); //core binary cell operation LibMatrixBincell.bincellOpInPlace(this, that, op); } public void incrementalAggregate(AggregateOperator aggOp, MatrixValue correction, MatrixValue newWithCorrection) throws DMLUnsupportedOperationException, DMLRuntimeException { //assert(aggOp.correctionExists); MatrixBlock cor = checkType(correction); MatrixBlock newWithCor = checkType(newWithCorrection); KahanObject buffer = new KahanObject(0, 0); if (aggOp.correctionLocation == CorrectionLocationType.LASTROW) { for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); buffer._correction = cor.quickGetValue(0, c); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.quickGetValue(r, c), newWithCor.quickGetValue(r + 1, c)); quickSetValue(r, c, buffer._sum); cor.quickSetValue(0, c, buffer._correction); } } else if (aggOp.correctionLocation == CorrectionLocationType.LASTCOLUMN) { if (aggOp.increOp.fn instanceof Builtin && (((Builtin) (aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MAXINDEX || ((Builtin) (aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MININDEX)) { // *** HACK ALERT *** HACK ALERT *** HACK ALERT *** // rowIndexMax() and its siblings don't fit very well into the standard // aggregate framework. We (ab)use the "correction factor" argument to // hold the maximum value in each row/column. // The execute() method for this aggregate takes as its argument // two candidates for the highest value. Bookkeeping about // indexes (return column/row index with highest value, breaking // ties in favor of higher indexes) is handled in this function. // Note that both versions of incrementalAggregate() contain // very similar blocks of special-case code. If one block is // modified, the other needs to be changed to match. for (int r = 0; r < rlen; r++) { double currMaxValue = cor.quickGetValue(r, 0); long newMaxIndex = (long) newWithCor.quickGetValue(r, 0); double newMaxValue = newWithCor.quickGetValue(r, 1); double update = aggOp.increOp.fn.execute(newMaxValue, currMaxValue); if (2.0 == update) { // Return value of 2 ==> both values the same, break ties // in favor of higher index. long curMaxIndex = (long) quickGetValue(r, 0); quickSetValue(r, 0, Math.max(curMaxIndex, newMaxIndex)); } else if (1.0 == update) { // Return value of 1 ==> new value is better; use its index quickSetValue(r, 0, newMaxIndex); cor.quickSetValue(r, 0, newMaxValue); } else { // Other return value ==> current answer is best } } // *** END HACK *** } else { for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); buffer._correction = cor.quickGetValue(r, 0); ; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.quickGetValue(r, c), newWithCor.quickGetValue(r, c + 1)); quickSetValue(r, c, buffer._sum); cor.quickSetValue(r, 0, buffer._correction); } } } else if (aggOp.correctionLocation == CorrectionLocationType.NONE) { //e.g., ak+ kahan plus as used in sum, mapmult, mmcj and tsmm if (aggOp.increOp.fn instanceof KahanPlus) { LibMatrixAgg.aggregateBinaryMatrix(newWithCor, this, cor); } else { if (newWithCor.isInSparseFormat() && aggOp.sparseSafe) //SPARSE { SparseRow[] bRows = newWithCor.getSparseRows(); if (bRows == null) //early abort on empty block return; for (int r = 0; r < Math.min(rlen, bRows.length); r++) { SparseRow brow = bRows[r]; if (brow != null && !brow.isEmpty()) { int blen = brow.size(); int[] bix = brow.getIndexContainer(); double[] bvals = brow.getValueContainer(); for (int j = 0; j < blen; j++) { int c = bix[j]; buffer._sum = this.quickGetValue(r, c); buffer._correction = cor.quickGetValue(r, c); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, bvals[j]); quickSetValue(r, c, buffer._sum); cor.quickSetValue(r, c, buffer._correction); } } } } else //DENSE or SPARSE (!sparsesafe) { for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); buffer._correction = cor.quickGetValue(r, c); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.quickGetValue(r, c)); quickSetValue(r, c, buffer._sum); cor.quickSetValue(r, c, buffer._correction); } } //change representation if required //(note since ak+ on blocks is currently only applied in MR, hence no need to account for this in mem estimates) examSparsity(); } } else if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOROWS) { double n, n2, mu2; for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); n = cor.quickGetValue(0, c); buffer._correction = cor.quickGetValue(1, c); mu2 = newWithCor.quickGetValue(r, c); n2 = newWithCor.quickGetValue(r + 1, c); n = n + n2; double toadd = (mu2 - buffer._sum) * n2 / n; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, toadd); quickSetValue(r, c, buffer._sum); cor.quickSetValue(0, c, n); cor.quickSetValue(1, c, buffer._correction); } } else if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOCOLUMNS) { double n, n2, mu2; for (int r = 0; r < rlen; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); n = cor.quickGetValue(r, 0); buffer._correction = cor.quickGetValue(r, 1); mu2 = newWithCor.quickGetValue(r, c); n2 = newWithCor.quickGetValue(r, c + 1); n = n + n2; double toadd = (mu2 - buffer._sum) * n2 / n; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, toadd); quickSetValue(r, c, buffer._sum); cor.quickSetValue(r, 0, n); cor.quickSetValue(r, 1, buffer._correction); } } else throw new DMLRuntimeException("unrecognized correctionLocation: " + aggOp.correctionLocation); } public void incrementalAggregate(AggregateOperator aggOp, MatrixValue newWithCorrection) throws DMLUnsupportedOperationException, DMLRuntimeException { //assert(aggOp.correctionExists); MatrixBlock newWithCor = checkType(newWithCorrection); KahanObject buffer = new KahanObject(0, 0); if (aggOp.correctionLocation == CorrectionLocationType.LASTROW) { if (aggOp.increOp.fn instanceof KahanPlus) { LibMatrixAgg.aggregateBinaryMatrix(newWithCor, this, aggOp); } else { for (int r = 0; r < rlen - 1; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); buffer._correction = this.quickGetValue(r + 1, c); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.quickGetValue(r, c), newWithCor.quickGetValue(r + 1, c)); quickSetValue(r, c, buffer._sum); quickSetValue(r + 1, c, buffer._correction); } } } else if (aggOp.correctionLocation == CorrectionLocationType.LASTCOLUMN) { if (aggOp.increOp.fn instanceof Builtin && (((Builtin) (aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MAXINDEX || ((Builtin) (aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MININDEX)) { // *** HACK ALERT *** HACK ALERT *** HACK ALERT *** // rowIndexMax() and its siblings don't fit very well into the standard // aggregate framework. We (ab)use the "correction factor" argument to // hold the maximum value in each row/column. // The execute() method for this aggregate takes as its argument // two candidates for the highest value. Bookkeeping about // indexes (return column/row index with highest value, breaking // ties in favor of higher indexes) is handled in this function. // Note that both versions of incrementalAggregate() contain // very similar blocks of special-case code. If one block is // modified, the other needs to be changed to match. for (int r = 0; r < rlen; r++) { double currMaxValue = quickGetValue(r, 1); long newMaxIndex = (long) newWithCor.quickGetValue(r, 0); double newMaxValue = newWithCor.quickGetValue(r, 1); double update = aggOp.increOp.fn.execute(newMaxValue, currMaxValue); if (2.0 == update) { // Return value of 2 ==> both values the same, break ties // in favor of higher index. long curMaxIndex = (long) quickGetValue(r, 0); quickSetValue(r, 0, Math.max(curMaxIndex, newMaxIndex)); } else if (1.0 == update) { // Return value of 1 ==> new value is better; use its index quickSetValue(r, 0, newMaxIndex); quickSetValue(r, 1, newMaxValue); } else { // Other return value ==> current answer is best } } // *** END HACK *** } else { if (aggOp.increOp.fn instanceof KahanPlus) { LibMatrixAgg.aggregateBinaryMatrix(newWithCor, this, aggOp); } else { for (int r = 0; r < rlen; r++) for (int c = 0; c < clen - 1; c++) { buffer._sum = this.quickGetValue(r, c); buffer._correction = this.quickGetValue(r, c + 1); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.quickGetValue(r, c), newWithCor.quickGetValue(r, c + 1)); quickSetValue(r, c, buffer._sum); quickSetValue(r, c + 1, buffer._correction); } } } } /*else if(aggOp.correctionLocation==0) { for(int r=0; r<rlen; r++) for(int c=0; c<clen; c++) { //buffer._sum=this.getValue(r, c); //buffer._correction=0; //buffer=(KahanObject) aggOp.increOp.fn.execute(buffer, newWithCor.getValue(r, c)); setValue(r, c, this.getValue(r, c)+newWithCor.getValue(r, c)); } }*/else if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOROWS) { double n, n2, mu2; for (int r = 0; r < rlen - 2; r++) for (int c = 0; c < clen; c++) { buffer._sum = this.quickGetValue(r, c); n = this.quickGetValue(r + 1, c); buffer._correction = this.quickGetValue(r + 2, c); mu2 = newWithCor.quickGetValue(r, c); n2 = newWithCor.quickGetValue(r + 1, c); n = n + n2; double toadd = (mu2 - buffer._sum) * n2 / n; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, toadd); quickSetValue(r, c, buffer._sum); quickSetValue(r + 1, c, n); quickSetValue(r + 2, c, buffer._correction); } } else if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOCOLUMNS) { double n, n2, mu2; for (int r = 0; r < rlen; r++) for (int c = 0; c < clen - 2; c++) { buffer._sum = this.quickGetValue(r, c); n = this.quickGetValue(r, c + 1); buffer._correction = this.quickGetValue(r, c + 2); mu2 = newWithCor.quickGetValue(r, c); n2 = newWithCor.quickGetValue(r, c + 1); n = n + n2; double toadd = (mu2 - buffer._sum) * n2 / n; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, toadd); quickSetValue(r, c, buffer._sum); quickSetValue(r, c + 1, n); quickSetValue(r, c + 2, buffer._correction); } } else throw new DMLRuntimeException("unrecognized correctionLocation: " + aggOp.correctionLocation); } @Override public MatrixValue reorgOperations(ReorgOperator op, MatrixValue ret, int startRow, int startColumn, int length) throws DMLRuntimeException { if (!(op.fn instanceof SwapIndex || op.fn instanceof DiagIndex || op.fn instanceof SortIndex)) throw new DMLRuntimeException("the current reorgOperations cannot support: " + op.fn.getClass() + "."); MatrixBlock result = checkType(ret); //compute output dimensions and sparsity flag CellIndex tempCellIndex = new CellIndex(-1, -1); op.fn.computeDimension(rlen, clen, tempCellIndex); boolean sps = evalSparseFormatInMemory(tempCellIndex.row, tempCellIndex.column, nonZeros); //prepare output matrix block w/ right meta data if (result == null) result = new MatrixBlock(tempCellIndex.row, tempCellIndex.column, sps, nonZeros); else result.reset(tempCellIndex.row, tempCellIndex.column, sps, nonZeros); if (LibMatrixReorg.isSupportedReorgOperator(op)) { //SPECIAL case (operators with special performance requirements, //or size-dependent special behavior) //currently supported opcodes: r', rdiag, rsort LibMatrixReorg.reorg(this, result, op); } else { //GENERIC case (any reorg operator) CellIndex temp = new CellIndex(0, 0); if (sparse) { if (sparseRows != null) { for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; int[] cols = sparseRows[r].getIndexContainer(); double[] values = sparseRows[r].getValueContainer(); for (int i = 0; i < sparseRows[r].size(); i++) { tempCellIndex.set(r, cols[i]); op.fn.execute(tempCellIndex, temp); result.appendValue(temp.row, temp.column, values[i]); } } } } else { if (denseBlock != null) { if (result.isInSparseFormat()) //SPARSE<-DENSE { double[] a = denseBlock; for (int i = 0, aix = 0; i < rlen; i++) for (int j = 0; j < clen; j++, aix++) { temp.set(i, j); op.fn.execute(temp, temp); result.appendValue(temp.row, temp.column, a[aix]); } } else //DENSE<-DENSE { result.allocateDenseBlock(); Arrays.fill(result.denseBlock, 0); double[] a = denseBlock; double[] c = result.denseBlock; int n = result.clen; for (int i = 0, aix = 0; i < rlen; i++) for (int j = 0; j < clen; j++, aix++) { temp.set(i, j); op.fn.execute(temp, temp); c[temp.row * n + temp.column] = a[aix]; } result.nonZeros = nonZeros; } } } } return result; } /** * * @param that * @param ret * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public MatrixBlock appendOperations(MatrixBlock that, MatrixBlock ret) throws DMLUnsupportedOperationException, DMLRuntimeException { //default append-cbind return appendOperations(that, ret, true); } /** * * @param that * @param ret * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public MatrixBlock appendOperations(MatrixBlock that, MatrixBlock ret, boolean cbind) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock result = checkType(ret); final int m = cbind ? rlen : rlen + that.rlen; final int n = cbind ? clen + that.clen : clen; final long nnz = nonZeros + that.nonZeros; boolean sp = evalSparseFormatInMemory(m, n, nnz); //init result matrix if (result == null) result = new MatrixBlock(m, n, sp, nnz); else result.reset(m, n, sp, nnz); //core append operation //copy left and right input into output if (!result.sparse) //DENSE { if (cbind) { result.copy(0, m - 1, 0, clen - 1, this, false); result.copy(0, m - 1, clen, n - 1, that, false); } else { //rbind result.copy(0, rlen - 1, 0, n - 1, this, false); result.copy(rlen, m - 1, 0, n - 1, that, false); } } else //SPARSE { //adjust sparse rows if required if (!this.isEmptyBlock(false) || !that.isEmptyBlock(false)) result.allocateSparseRowsBlock(); result.appendToSparse(this, 0, 0); if (cbind) result.appendToSparse(that, 0, clen); else //rbind result.appendToSparse(that, rlen, 0); } //update meta data result.nonZeros = nnz; return result; } /** * * @param out * @param tstype * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock transposeSelfMatrixMultOperations(MatrixBlock out, MMTSJType tstype) throws DMLRuntimeException, DMLUnsupportedOperationException { return transposeSelfMatrixMultOperations(out, tstype, 1); } /** * * @param out * @param tstype * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock transposeSelfMatrixMultOperations(MatrixBlock out, MMTSJType tstype, int k) throws DMLRuntimeException, DMLUnsupportedOperationException { //check for transpose type if (!(tstype == MMTSJType.LEFT || tstype == MMTSJType.RIGHT)) throw new DMLRuntimeException("Invalid MMTSJ type '" + tstype.toString() + "'."); //setup meta data boolean leftTranspose = (tstype == MMTSJType.LEFT); int dim = leftTranspose ? clen : rlen; //create output matrix block if (out == null) out = new MatrixBlock(dim, dim, false); else out.reset(dim, dim, false); //compute matrix mult if (k > 1) LibMatrixMult.matrixMultTransposeSelf(this, out, leftTranspose, k); else LibMatrixMult.matrixMultTransposeSelf(this, out, leftTranspose); return out; } /** * * @param v * @param w * @param out * @param ctype * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock chainMatrixMultOperations(MatrixBlock v, MatrixBlock w, MatrixBlock out, ChainType ctype) throws DMLRuntimeException, DMLUnsupportedOperationException { return chainMatrixMultOperations(v, w, out, ctype, 1); } /** * * @param v * @param w * @param out * @param ctype * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock chainMatrixMultOperations(MatrixBlock v, MatrixBlock w, MatrixBlock out, ChainType ctype, int k) throws DMLRuntimeException, DMLUnsupportedOperationException { //check for transpose type if (!(ctype == ChainType.XtXv || ctype == ChainType.XtwXv)) throw new DMLRuntimeException("Invalid mmchain type '" + ctype.toString() + "'."); //check for matching dimensions if (this.getNumColumns() != v.getNumRows()) throw new DMLRuntimeException("Dimensions mismatch on mmchain operation (" + this.getNumColumns() + " != " + v.getNumRows() + ")"); if (v != null && v.getNumColumns() != 1) throw new DMLRuntimeException( "Invalid input vector (column vector expected, but ncol=" + v.getNumColumns() + ")"); if (w != null && w.getNumColumns() != 1) throw new DMLRuntimeException( "Invalid weight vector (column vector expected, but ncol=" + w.getNumColumns() + ")"); //prepare result if (out != null) out.reset(clen, 1, false); else out = new MatrixBlock(clen, 1, false); //compute matrix mult if (k > 1) LibMatrixMult.matrixMultChain(this, v, w, out, ctype, k); else LibMatrixMult.matrixMultChain(this, v, w, out, ctype); return out; } /** * * @param m1Val * @param m2Val * @param out1Val * @param out2Val * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public void permutationMatrixMultOperations(MatrixValue m2Val, MatrixValue out1Val, MatrixValue out2Val) throws DMLRuntimeException, DMLUnsupportedOperationException { permutationMatrixMultOperations(m2Val, out1Val, out2Val, 1); } /** * * @param m1Val * @param m2Val * @param out1Val * @param out2Val * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public void permutationMatrixMultOperations(MatrixValue m2Val, MatrixValue out1Val, MatrixValue out2Val, int k) throws DMLRuntimeException, DMLUnsupportedOperationException { //check input types and dimensions MatrixBlock m2 = checkType(m2Val); MatrixBlock ret1 = checkType(out1Val); MatrixBlock ret2 = checkType(out2Val); if (this.rlen != m2.rlen) throw new RuntimeException("Dimensions do not match for permutation matrix multiplication (" + this.rlen + "!=" + m2.rlen + ")."); //compute permutation matrix multiplication if (k > 1) LibMatrixMult.matrixMultPermute(this, m2, ret1, ret2, k); else LibMatrixMult.matrixMultPermute(this, m2, ret1, ret2); } /** * Method to perform leftIndexing operation for a given lower and upper bounds in row and column dimensions. * Updated matrix is returned as the output. * * Operations to be performed: * 1) result=this; * 2) result[rowLower:rowUpper, colLower:colUpper] = rhsMatrix; * * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock leftIndexingOperations(MatrixBlock rhsMatrix, int rl, int ru, int cl, int cu, MatrixBlock ret, boolean inplace) throws DMLRuntimeException, DMLUnsupportedOperationException { // Check the validity of bounds if (rl < 0 || rl >= getNumRows() || ru < rl || ru >= getNumRows() || cl < 0 || cu >= getNumColumns() || cu < cl || cu >= getNumColumns()) { throw new DMLRuntimeException("Invalid values for matrix indexing: [" + (rl + 1) + ":" + (ru + 1) + "," + (cl + 1) + ":" + (cu + 1) + "] " + "must be within matrix dimensions [" + getNumRows() + "," + getNumColumns() + "]."); } if ((ru - rl + 1) < rhsMatrix.getNumRows() || (cu - cl + 1) < rhsMatrix.getNumColumns()) { throw new DMLRuntimeException("Invalid values for matrix indexing: " + "dimensions of the source matrix [" + rhsMatrix.getNumRows() + "x" + rhsMatrix.getNumColumns() + "] " + "do not match the shape of the matrix specified by indices [" + (rl + 1) + ":" + (ru + 1) + ", " + (cl + 1) + ":" + (cu + 1) + "]."); } MatrixBlock result = ret; boolean sp = estimateSparsityOnLeftIndexing(rlen, clen, nonZeros, rhsMatrix.getNumRows(), rhsMatrix.getNumColumns(), rhsMatrix.getNonZeros()); if (!inplace) //general case { if (result == null) result = new MatrixBlock(rlen, clen, sp); else result.reset(rlen, clen, sp); result.copy(this, sp); } else //update in-place { //use current block as in-place result result = this; //ensure that the current block adheres to the sparsity estimate //and thus implicitly the memory budget used by the compiler if (result.sparse && !sp) result.sparseToDense(); else if (!result.sparse && sp) result.denseToSparse(); } //NOTE conceptually we could directly use a zeroout and copy(..., false) but // since this was factors slower, we still use a full copy and subsequently // copy(..., true) - however, this can be changed in the future once we // improved the performance of zeroout. //result = (MatrixBlockDSM) zeroOutOperations(result, new IndexRange(rowLower,rowUpper, colLower, colUpper ), false); MatrixBlock src = (MatrixBlock) rhsMatrix; if (rl == ru && cl == cu) //specific case: cell update { //copy single value and update nnz result.quickSetValue(rl, cl, src.quickGetValue(0, 0)); } else //general case { //copy submatrix into result result.copy(rl, ru, cl, cu, src, true); } return result; } /** * Explicitly allow left indexing for scalars. Note: This operation is now 0-based. * * * Operations to be performed: * 1) result=this; * 2) result[row,column] = scalar.getDoubleValue(); * * @param scalar * @param row * @param col * @param ret * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock leftIndexingOperations(ScalarObject scalar, int rl, int cl, MatrixBlock ret, boolean inplace) throws DMLRuntimeException, DMLUnsupportedOperationException { double inVal = scalar.getDoubleValue(); boolean sp = estimateSparsityOnLeftIndexing(rlen, clen, nonZeros, 1, 1, (inVal != 0) ? 1 : 0); if (!inplace) //general case { if (ret == null) ret = new MatrixBlock(rlen, clen, sp); else ret.reset(rlen, clen, sp); ret.copy(this, sp); } else //update in-place ret = this; ret.quickSetValue(rl, cl, inVal); return ret; } /** * Method to perform rangeReIndex operation for a given lower and upper bounds in row and column dimensions. * Extracted submatrix is returned as "result". Note: This operation is now 0-based. * * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock sliceOperations(int rl, int ru, int cl, int cu, MatrixBlock ret) throws DMLRuntimeException { // check the validity of bounds if (rl < 0 || rl >= getNumRows() || ru < rl || ru >= getNumRows() || cl < 0 || cu >= getNumColumns() || cu < cl || cu >= getNumColumns()) { throw new DMLRuntimeException("Invalid values for matrix indexing: [" + (rl + 1) + ":" + (ru + 1) + "," + (cl + 1) + ":" + (cu + 1) + "] " + "must be within matrix dimensions [" + getNumRows() + "," + getNumColumns() + "]"); } // Output matrix will have the same sparsity as that of the input matrix. // (assuming a uniform distribution of non-zeros in the input) MatrixBlock result = checkType(ret); long estnnz = (long) ((double) this.nonZeros / rlen / clen * (ru - rl + 1) * (cu - cl + 1)); boolean result_sparsity = this.sparse && MatrixBlock.evalSparseFormatInMemory(ru - rl + 1, cu - cl + 1, estnnz); if (result == null) result = new MatrixBlock(ru - rl + 1, cu - cl + 1, result_sparsity, estnnz); else result.reset(ru - rl + 1, cu - cl + 1, result_sparsity, estnnz); // actual slice operation if (rl == 0 && ru == rlen - 1 && cl == 0 && cu == clen - 1) { // copy if entire matrix required result.copy(this); } else //general case { //core slicing operation (nnz maintained internally) if (sparse) sliceSparse(rl, ru, cl, cu, result); else sliceDense(rl, ru, cl, cu, result); } return result; } /** * * @param rl * @param ru * @param cl * @param cu * @param dest * @throws DMLRuntimeException */ private void sliceSparse(int rl, int ru, int cl, int cu, MatrixBlock dest) throws DMLRuntimeException { //check for early abort if (isEmptyBlock(false)) return; if (cl == cu) //COLUMN VECTOR { //note: always dense dest dest.allocateDenseBlock(); for (int i = rl; i <= ru; i++) { SparseRow arow = sparseRows[i]; if (arow != null && !arow.isEmpty()) { double val = arow.get(cl); if (val != 0) { dest.denseBlock[i - rl] = val; dest.nonZeros++; } } } } else if (rl == ru && cl == 0 && cu == clen - 1) //ROW VECTOR { //note: always sparse dest, but also works for dense dest.appendRow(0, sparseRows[rl]); } else //general case (sparse/dense dest) { for (int i = rl; i <= ru; i++) if (sparseRows[i] != null && !sparseRows[i].isEmpty()) { SparseRow arow = sparseRows[i]; int alen = arow.size(); int[] aix = arow.getIndexContainer(); double[] avals = arow.getValueContainer(); int astart = (cl > 0) ? arow.searchIndexesFirstGTE(cl) : 0; if (astart != -1) for (int j = astart; j < alen && aix[j] <= cu; j++) dest.appendValue(i - rl, aix[j] - cl, avals[j]); } } } /** * * @param rl * @param ru * @param cl * @param cu * @param dest * @throws DMLRuntimeException */ private void sliceDense(int rl, int ru, int cl, int cu, MatrixBlock dest) throws DMLRuntimeException { //ensure allocated input/output blocks if (denseBlock == null) return; dest.allocateDenseBlock(); //indexing operation if (cl == cu) //COLUMN INDEXING { if (clen == 1) //vector -> vector { System.arraycopy(denseBlock, rl, dest.denseBlock, 0, ru - rl + 1); } else //matrix -> vector { //IBM JVM bug (JDK7) causes crash for certain cl/cu values (e.g., divide by zero for 4) //for( int i=rl*clen+cl, ix=0; i<=ru*clen+cu; i+=clen, ix++ ) // dest.denseBlock[ix] = denseBlock[i]; int len = clen; for (int i = rl * len + cl, ix = 0; i <= ru * len + cu; i += len, ix++) dest.denseBlock[ix] = denseBlock[i]; } } else // GENERAL RANGE INDEXING { //IBM JVM bug (JDK7) causes crash for certain cl/cu values (e.g., divide by zero for 4) //for(int i = rl, ix1 = rl*clen+cl, ix2=0; i <= ru; i++, ix1+=clen, ix2+=dest.clen) // System.arraycopy(denseBlock, ix1, dest.denseBlock, ix2, dest.clen); int len1 = clen; int len2 = dest.clen; for (int i = rl, ix1 = rl * len1 + cl, ix2 = 0; i <= ru; i++, ix1 += len1, ix2 += len2) System.arraycopy(denseBlock, ix1, dest.denseBlock, ix2, len2); } //compute nnz of output (not maintained due to native calls) dest.recomputeNonZeros(); } public void sliceOperations(ArrayList<IndexedMatrixValue> outlist, IndexRange range, int rowCut, int colCut, int normalBlockRowFactor, int normalBlockColFactor, int boundaryRlen, int boundaryClen) { MatrixBlock topleft = null, topright = null, bottomleft = null, bottomright = null; Iterator<IndexedMatrixValue> p = outlist.iterator(); int blockRowFactor = normalBlockRowFactor, blockColFactor = normalBlockColFactor; if (rowCut > range.rowEnd) blockRowFactor = boundaryRlen; if (colCut > range.colEnd) blockColFactor = boundaryClen; int minrowcut = (int) Math.min(rowCut, range.rowEnd); int mincolcut = (int) Math.min(colCut, range.colEnd); int maxrowcut = (int) Math.max(rowCut, range.rowStart); int maxcolcut = (int) Math.max(colCut, range.colStart); if (range.rowStart < rowCut && range.colStart < colCut) { topleft = (MatrixBlock) p.next().getValue(); //topleft.reset(blockRowFactor, blockColFactor, // checkSparcityOnSlide(rowCut-(int)range.rowStart, colCut-(int)range.colStart, blockRowFactor, blockColFactor)); topleft.reset(blockRowFactor, blockColFactor, estimateSparsityOnSlice(minrowcut - (int) range.rowStart, mincolcut - (int) range.colStart, blockRowFactor, blockColFactor)); } if (range.rowStart < rowCut && range.colEnd >= colCut) { topright = (MatrixBlock) p.next().getValue(); topright.reset(blockRowFactor, boundaryClen, estimateSparsityOnSlice(minrowcut - (int) range.rowStart, (int) range.colEnd - maxcolcut + 1, blockRowFactor, boundaryClen)); } if (range.rowEnd >= rowCut && range.colStart < colCut) { bottomleft = (MatrixBlock) p.next().getValue(); bottomleft.reset(boundaryRlen, blockColFactor, estimateSparsityOnSlice((int) range.rowEnd - maxrowcut + 1, mincolcut - (int) range.colStart, boundaryRlen, blockColFactor)); } if (range.rowEnd >= rowCut && range.colEnd >= colCut) { bottomright = (MatrixBlock) p.next().getValue(); bottomright.reset(boundaryRlen, boundaryClen, estimateSparsityOnSlice((int) range.rowEnd - maxrowcut + 1, (int) range.colEnd - maxcolcut + 1, boundaryRlen, boundaryClen)); } if (sparse) { if (sparseRows != null) { int r = (int) range.rowStart; for (; r < Math.min(Math.min(rowCut, sparseRows.length), range.rowEnd + 1); r++) sliceHelp(r, range, colCut, topleft, topright, normalBlockRowFactor - rowCut, normalBlockRowFactor, normalBlockColFactor); for (; r <= Math.min(range.rowEnd, sparseRows.length - 1); r++) sliceHelp(r, range, colCut, bottomleft, bottomright, -rowCut, normalBlockRowFactor, normalBlockColFactor); //System.out.println("in: \n"+this); //System.out.println("outlist: \n"+outlist); } } else { if (denseBlock != null) { int i = ((int) range.rowStart) * clen; int r = (int) range.rowStart; for (; r < Math.min(rowCut, range.rowEnd + 1); r++) { int c = (int) range.colStart; for (; c < Math.min(colCut, range.colEnd + 1); c++) topleft.appendValue(r + normalBlockRowFactor - rowCut, c + normalBlockColFactor - colCut, denseBlock[i + c]); for (; c <= range.colEnd; c++) topright.appendValue(r + normalBlockRowFactor - rowCut, c - colCut, denseBlock[i + c]); i += clen; } for (; r <= range.rowEnd; r++) { int c = (int) range.colStart; for (; c < Math.min(colCut, range.colEnd + 1); c++) bottomleft.appendValue(r - rowCut, c + normalBlockColFactor - colCut, denseBlock[i + c]); for (; c <= range.colEnd; c++) bottomright.appendValue(r - rowCut, c - colCut, denseBlock[i + c]); i += clen; } } } } private void sliceHelp(int r, IndexRange range, int colCut, MatrixBlock left, MatrixBlock right, int rowOffset, int normalBlockRowFactor, int normalBlockColFactor) { if (sparseRows[r] == null) return; int[] cols = sparseRows[r].getIndexContainer(); double[] values = sparseRows[r].getValueContainer(); int start = sparseRows[r].searchIndexesFirstGTE((int) range.colStart); if (start < 0) return; int end = sparseRows[r].searchIndexesFirstLTE((int) range.colEnd); if (end < 0 || start > end) return; //actual slice operation for (int i = start; i <= end; i++) { if (cols[i] < colCut) left.appendValue(r + rowOffset, cols[i] + normalBlockColFactor - colCut, values[i]); else right.appendValue(r + rowOffset, cols[i] - colCut, values[i]); } } @Override //This the append operations for MR side //nextNCol is the number columns for the block right of block v2 public void appendOperations(MatrixValue v2, ArrayList<IndexedMatrixValue> outlist, int blockRowFactor, int blockColFactor, boolean cbind, boolean m2IsLast, int nextNCol) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock m2 = (MatrixBlock) v2; //case 1: copy lhs and rhs to output if (cbind && clen == blockColFactor || !cbind && rlen == blockRowFactor) { ((MatrixBlock) outlist.get(0).getValue()).copy(this); ((MatrixBlock) outlist.get(1).getValue()).copy(m2); } //case 2: append part of rhs to lhs, append to 2nd output if necessary else { //single output block (via plain append operation) if (cbind && clen + m2.clen < blockColFactor || !cbind && rlen + m2.rlen < blockRowFactor) { appendOperations(m2, (MatrixBlock) outlist.get(0).getValue(), cbind); } //two output blocks (via slice and append) else { //prepare output block 1 MatrixBlock ret1 = (MatrixBlock) outlist.get(0).getValue(); int lrlen1 = cbind ? rlen - 1 : blockRowFactor - rlen - 1; int lclen1 = cbind ? blockColFactor - clen - 1 : clen - 1; MatrixBlock tmp1 = m2.sliceOperations(0, lrlen1, 0, lclen1, new MatrixBlock()); appendOperations(tmp1, ret1, cbind); //prepare output block 2 MatrixBlock ret2 = (MatrixBlock) outlist.get(1).getValue(); if (cbind) m2.sliceOperations(0, rlen - 1, lclen1 + 1, m2.clen - 1, ret2); else m2.sliceOperations(lrlen1 + 1, m2.rlen - 1, 0, clen - 1, ret2); } } } /** * */ public MatrixValue zeroOutOperations(MatrixValue result, IndexRange range, boolean complementary) throws DMLUnsupportedOperationException, DMLRuntimeException { checkType(result); double currentSparsity = (double) nonZeros / (double) rlen / (double) clen; double estimatedSps = currentSparsity * (double) (range.rowEnd - range.rowStart + 1) * (double) (range.colEnd - range.colStart + 1) / (double) rlen / (double) clen; if (!complementary) estimatedSps = currentSparsity - estimatedSps; boolean lsparse = evalSparseFormatInMemory(rlen, clen, (long) (estimatedSps * rlen * clen)); if (result == null) result = new MatrixBlock(rlen, clen, lsparse, (int) (estimatedSps * rlen * clen)); else result.reset(rlen, clen, lsparse, (int) (estimatedSps * rlen * clen)); if (sparse) { if (sparseRows != null) { if (!complementary)//if zero out { for (int r = 0; r < Math.min((int) range.rowStart, sparseRows.length); r++) ((MatrixBlock) result).appendRow(r, sparseRows[r]); for (int r = Math.min((int) range.rowEnd + 1, sparseRows.length); r < Math.min(rlen, sparseRows.length); r++) ((MatrixBlock) result).appendRow(r, sparseRows[r]); } for (int r = (int) range.rowStart; r <= Math.min(range.rowEnd, sparseRows.length - 1); r++) { if (sparseRows[r] == null) continue; int[] cols = sparseRows[r].getIndexContainer(); double[] values = sparseRows[r].getValueContainer(); if (complementary)//if selection { int start = sparseRows[r].searchIndexesFirstGTE((int) range.colStart); if (start < 0) continue; int end = sparseRows[r].searchIndexesFirstGT((int) range.colEnd); if (end < 0 || start > end) continue; for (int i = start; i < end; i++) { ((MatrixBlock) result).appendValue(r, cols[i], values[i]); } } else { int start = sparseRows[r].searchIndexesFirstGTE((int) range.colStart); if (start < 0) start = sparseRows[r].size(); int end = sparseRows[r].searchIndexesFirstGT((int) range.colEnd); if (end < 0) end = sparseRows[r].size(); for (int i = 0; i < start; i++) { ((MatrixBlock) result).appendValue(r, cols[i], values[i]); } for (int i = end; i < sparseRows[r].size(); i++) { ((MatrixBlock) result).appendValue(r, cols[i], values[i]); } } } } } else { if (denseBlock != null) { if (complementary)//if selection { int offset = ((int) range.rowStart) * clen; for (int r = (int) range.rowStart; r <= range.rowEnd; r++) { for (int c = (int) range.colStart; c <= range.colEnd; c++) ((MatrixBlock) result).appendValue(r, c, denseBlock[offset + c]); offset += clen; } } else { int offset = 0; int r = 0; for (; r < (int) range.rowStart; r++) for (int c = 0; c < clen; c++, offset++) ((MatrixBlock) result).appendValue(r, c, denseBlock[offset]); for (; r <= (int) range.rowEnd; r++) { for (int c = 0; c < (int) range.colStart; c++) ((MatrixBlock) result).appendValue(r, c, denseBlock[offset + c]); for (int c = (int) range.colEnd + 1; c < clen; c++) ((MatrixBlock) result).appendValue(r, c, denseBlock[offset + c]); offset += clen; } for (; r < rlen; r++) for (int c = 0; c < clen; c++, offset++) ((MatrixBlock) result).appendValue(r, c, denseBlock[offset]); } } } return result; } public MatrixValue aggregateUnaryOperations(AggregateUnaryOperator op, MatrixValue result, int blockingFactorRow, int blockingFactorCol, MatrixIndexes indexesIn) throws DMLUnsupportedOperationException, DMLRuntimeException { return aggregateUnaryOperations(op, result, blockingFactorRow, blockingFactorCol, indexesIn, false); } public MatrixValue aggregateUnaryOperations(AggregateUnaryOperator op, MatrixValue result, int blockingFactorRow, int blockingFactorCol, MatrixIndexes indexesIn, boolean inCP) throws DMLUnsupportedOperationException, DMLRuntimeException { CellIndex tempCellIndex = new CellIndex(-1, -1); op.indexFn.computeDimension(rlen, clen, tempCellIndex); if (op.aggOp.correctionExists) { switch (op.aggOp.correctionLocation) { case LASTROW: tempCellIndex.row++; break; case LASTCOLUMN: tempCellIndex.column++; break; case LASTTWOROWS: tempCellIndex.row += 2; break; case LASTTWOCOLUMNS: tempCellIndex.column += 2; break; default: throw new DMLRuntimeException("unrecognized correctionLocation: " + op.aggOp.correctionLocation); } } if (result == null) result = new MatrixBlock(tempCellIndex.row, tempCellIndex.column, false); else result.reset(tempCellIndex.row, tempCellIndex.column, false); MatrixBlock ret = (MatrixBlock) result; if (LibMatrixAgg.isSupportedUnaryAggregateOperator(op)) { if (op.getNumThreads() > 1) LibMatrixAgg.aggregateUnaryMatrix(this, ret, op, op.getNumThreads()); else LibMatrixAgg.aggregateUnaryMatrix(this, ret, op); LibMatrixAgg.recomputeIndexes(ret, op, blockingFactorRow, blockingFactorCol, indexesIn); } else if (op.sparseSafe) sparseAggregateUnaryHelp(op, ret, blockingFactorRow, blockingFactorCol, indexesIn); else denseAggregateUnaryHelp(op, ret, blockingFactorRow, blockingFactorCol, indexesIn); if (op.aggOp.correctionExists && inCP) ((MatrixBlock) result).dropLastRowsOrColums(op.aggOp.correctionLocation); return ret; } private void sparseAggregateUnaryHelp(AggregateUnaryOperator op, MatrixBlock result, int blockingFactorRow, int blockingFactorCol, MatrixIndexes indexesIn) throws DMLRuntimeException { //initialize result if (op.aggOp.initialValue != 0) result.resetDenseWithValue(result.rlen, result.clen, op.aggOp.initialValue); CellIndex tempCellIndex = new CellIndex(-1, -1); KahanObject buffer = new KahanObject(0, 0); int r = 0, c = 0; if (sparse) { if (sparseRows != null) { for (r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; int[] cols = sparseRows[r].getIndexContainer(); double[] values = sparseRows[r].getValueContainer(); for (int i = 0; i < sparseRows[r].size(); i++) { tempCellIndex.set(r, cols[i]); op.indexFn.execute(tempCellIndex, tempCellIndex); incrementalAggregateUnaryHelp(op.aggOp, result, tempCellIndex.row, tempCellIndex.column, values[i], buffer); } } } } else { if (denseBlock != null) { int limit = rlen * clen; for (int i = 0; i < limit; i++) { r = i / clen; c = i % clen; tempCellIndex.set(r, c); op.indexFn.execute(tempCellIndex, tempCellIndex); incrementalAggregateUnaryHelp(op.aggOp, result, tempCellIndex.row, tempCellIndex.column, denseBlock[i], buffer); } } } } private void denseAggregateUnaryHelp(AggregateUnaryOperator op, MatrixBlock result, int blockingFactorRow, int blockingFactorCol, MatrixIndexes indexesIn) throws DMLRuntimeException { //initialize if (op.aggOp.initialValue != 0) result.resetDenseWithValue(result.rlen, result.clen, op.aggOp.initialValue); CellIndex tempCellIndex = new CellIndex(-1, -1); KahanObject buffer = new KahanObject(0, 0); for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { tempCellIndex.set(i, j); op.indexFn.execute(tempCellIndex, tempCellIndex); if (op.aggOp.correctionExists && op.aggOp.correctionLocation == CorrectionLocationType.LASTCOLUMN && op.aggOp.increOp.fn instanceof Builtin && (((Builtin) (op.aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MAXINDEX || ((Builtin) (op.aggOp.increOp.fn)).bFunc == Builtin.BuiltinFunctionCode.MININDEX)) { double currMaxValue = result.quickGetValue(i, 1); long newMaxIndex = UtilFunctions.cellIndexCalculation(indexesIn.getColumnIndex(), blockingFactorCol, j); double newMaxValue = quickGetValue(i, j); double update = op.aggOp.increOp.fn.execute(newMaxValue, currMaxValue); //System.out.println("currV="+currMaxValue+",newV="+newMaxValue+",newIX="+newMaxIndex+",update="+update); if (update == 1) { result.quickSetValue(i, 0, newMaxIndex); result.quickSetValue(i, 1, newMaxValue); } } else incrementalAggregateUnaryHelp(op.aggOp, result, tempCellIndex.row, tempCellIndex.column, quickGetValue(i, j), buffer); } } private void incrementalAggregateUnaryHelp(AggregateOperator aggOp, MatrixBlock result, int row, int column, double newvalue, KahanObject buffer) throws DMLRuntimeException { if (aggOp.correctionExists) { if (aggOp.correctionLocation == CorrectionLocationType.LASTROW || aggOp.correctionLocation == CorrectionLocationType.LASTCOLUMN) { int corRow = row, corCol = column; if (aggOp.correctionLocation == CorrectionLocationType.LASTROW)//extra row corRow++; else if (aggOp.correctionLocation == CorrectionLocationType.LASTCOLUMN) corCol++; else throw new DMLRuntimeException("unrecognized correctionLocation: " + aggOp.correctionLocation); buffer._sum = result.quickGetValue(row, column); buffer._correction = result.quickGetValue(corRow, corCol); buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newvalue); result.quickSetValue(row, column, buffer._sum); result.quickSetValue(corRow, corCol, buffer._correction); } else if (aggOp.correctionLocation == CorrectionLocationType.NONE) { throw new DMLRuntimeException("unrecognized correctionLocation: " + aggOp.correctionLocation); } else// for mean { int corRow = row, corCol = column; int countRow = row, countCol = column; if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOROWS) { countRow++; corRow += 2; } else if (aggOp.correctionLocation == CorrectionLocationType.LASTTWOCOLUMNS) { countCol++; corCol += 2; } else throw new DMLRuntimeException("unrecognized correctionLocation: " + aggOp.correctionLocation); buffer._sum = result.quickGetValue(row, column); buffer._correction = result.quickGetValue(corRow, corCol); double count = result.quickGetValue(countRow, countCol) + 1.0; buffer = (KahanObject) aggOp.increOp.fn.execute(buffer, newvalue, count); result.quickSetValue(row, column, buffer._sum); result.quickSetValue(corRow, corCol, buffer._correction); result.quickSetValue(countRow, countCol, count); } } else { newvalue = aggOp.increOp.fn.execute(result.quickGetValue(row, column), newvalue); result.quickSetValue(row, column, newvalue); } } /** * * @param correctionLocation */ public void dropLastRowsOrColums(CorrectionLocationType correctionLocation) { //do nothing if (correctionLocation == CorrectionLocationType.NONE || correctionLocation == CorrectionLocationType.INVALID) { return; } //determine number of rows/cols to be removed int step = (correctionLocation == CorrectionLocationType.LASTTWOROWS || correctionLocation == CorrectionLocationType.LASTTWOCOLUMNS) ? 2 : 1; //e.g., colSums, colMeans, colMaxs, colMeans if (correctionLocation == CorrectionLocationType.LASTROW || correctionLocation == CorrectionLocationType.LASTTWOROWS) { if (sparse) //SPARSE { if (sparseRows != null) for (int i = 1; i <= step; i++) if (sparseRows[rlen - i] != null) this.nonZeros -= sparseRows[rlen - i].size(); } else //DENSE { if (denseBlock != null) for (int i = (rlen - step) * clen; i < rlen * clen; i++) if (denseBlock[i] != 0) this.nonZeros--; } //just need to shrink the dimension, the deleted rows won't be accessed rlen -= step; } //e.g., rowSums, rowsMeans, rowsMaxs, rowsMeans if (correctionLocation == CorrectionLocationType.LASTCOLUMN || correctionLocation == CorrectionLocationType.LASTTWOCOLUMNS) { if (sparse) //SPARSE { if (sparseRows != null) { for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) if (sparseRows[r] != null) { int newSize = sparseRows[r].searchIndexesFirstGTE(clen - step); if (newSize >= 0) { this.nonZeros -= sparseRows[r].size() - newSize; sparseRows[r].truncate(newSize); } } } } else //DENSE { if (this.denseBlock != null) { //the first row doesn't need to be copied int targetIndex = clen - step; int sourceOffset = clen; this.nonZeros = 0; for (int i = 0; i < targetIndex; i++) if (denseBlock[i] != 0) this.nonZeros++; //start from the 2nd row for (int r = 1; r < rlen; r++) { for (int c = 0; c < clen - step; c++) { if ((denseBlock[targetIndex] = denseBlock[sourceOffset + c]) != 0) this.nonZeros++; targetIndex++; } sourceOffset += clen; } } } clen -= step; } } /** * * @param op * @return * @throws DMLRuntimeException */ public CM_COV_Object cmOperations(CMOperator op) throws DMLRuntimeException { // dimension check for input column vectors if (this.getNumColumns() != 1) { throw new DMLRuntimeException("Central Moment can not be computed on [" + this.getNumRows() + "," + this.getNumColumns() + "] matrix."); } CM_COV_Object cmobj = new CM_COV_Object(); // empty block handling (important for result corretness, otherwise // we get a NaN due to 0/0 on reading out the required result) if (isEmptyBlock(false)) { op.fn.execute(cmobj, 0.0, getNumRows()); return cmobj; } int nzcount = 0; if (sparse && sparseRows != null) //SPARSE { for (int r = 0; r < Math.min(rlen, sparseRows.length); r++) { if (sparseRows[r] == null) continue; double[] values = sparseRows[r].getValueContainer(); for (int i = 0; i < sparseRows[r].size(); i++) { op.fn.execute(cmobj, values[i]); nzcount++; } } // account for zeros in the vector op.fn.execute(cmobj, 0.0, this.getNumRows() - nzcount); } else if (denseBlock != null) //DENSE { //always vector (see check above) for (int i = 0; i < rlen; i++) op.fn.execute(cmobj, denseBlock[i]); } return cmobj; } public CM_COV_Object cmOperations(CMOperator op, MatrixBlock weights) throws DMLRuntimeException { /* this._data must be a 1 dimensional vector */ if (this.getNumColumns() != 1 || weights.getNumColumns() != 1) { throw new DMLRuntimeException("Central Moment can be computed only on 1-dimensional column matrices."); } if (this.getNumRows() != weights.getNumRows() || this.getNumColumns() != weights.getNumColumns()) { throw new DMLRuntimeException("Covariance: Mismatching dimensions between input and weight matrices - " + "[" + this.getNumRows() + "," + this.getNumColumns() + "] != [" + weights.getNumRows() + "," + weights.getNumColumns() + "]"); } CM_COV_Object cmobj = new CM_COV_Object(); if (sparse && sparseRows != null) //SPARSE { for (int i = 0; i < rlen; i++) op.fn.execute(cmobj, this.quickGetValue(i, 0), weights.quickGetValue(i, 0)); } else if (denseBlock != null) //DENSE { //always vectors (see check above) if (!weights.sparse) { //both dense vectors (default case) if (weights.denseBlock != null) for (int i = 0; i < rlen; i++) op.fn.execute(cmobj, denseBlock[i], weights.denseBlock[i]); } else { for (int i = 0; i < rlen; i++) op.fn.execute(cmobj, denseBlock[i], weights.quickGetValue(i, 0)); } } return cmobj; } public CM_COV_Object covOperations(COVOperator op, MatrixBlock that) throws DMLRuntimeException { /* this._data must be a 1 dimensional vector */ if (this.getNumColumns() != 1 || that.getNumColumns() != 1) { throw new DMLRuntimeException("Covariance can be computed only on 1-dimensional column matrices."); } if (this.getNumRows() != that.getNumRows() || this.getNumColumns() != that.getNumColumns()) { throw new DMLRuntimeException("Covariance: Mismatching input matrix dimensions - " + "[" + this.getNumRows() + "," + this.getNumColumns() + "] != [" + that.getNumRows() + "," + that.getNumColumns() + "]"); } CM_COV_Object covobj = new CM_COV_Object(); if (sparse && sparseRows != null) //SPARSE { for (int i = 0; i < rlen; i++) op.fn.execute(covobj, this.quickGetValue(i, 0), that.quickGetValue(i, 0)); } else if (denseBlock != null) //DENSE { //always vectors (see check above) if (!that.sparse) { //both dense vectors (default case) if (that.denseBlock != null) for (int i = 0; i < rlen; i++) op.fn.execute(covobj, denseBlock[i], that.denseBlock[i]); } else { for (int i = 0; i < rlen; i++) op.fn.execute(covobj, denseBlock[i], that.quickGetValue(i, 0)); } } return covobj; } public CM_COV_Object covOperations(COVOperator op, MatrixBlock that, MatrixBlock weights) throws DMLRuntimeException { /* this._data must be a 1 dimensional vector */ if (this.getNumColumns() != 1 || that.getNumColumns() != 1 || weights.getNumColumns() != 1) { throw new DMLRuntimeException("Covariance can be computed only on 1-dimensional column matrices."); } if (this.getNumRows() != that.getNumRows() || this.getNumColumns() != that.getNumColumns()) { throw new DMLRuntimeException("Covariance: Mismatching input matrix dimensions - " + "[" + this.getNumRows() + "," + this.getNumColumns() + "] != [" + that.getNumRows() + "," + that.getNumColumns() + "]"); } if (this.getNumRows() != weights.getNumRows() || this.getNumColumns() != weights.getNumColumns()) { throw new DMLRuntimeException("Covariance: Mismatching dimensions between input and weight matrices - " + "[" + this.getNumRows() + "," + this.getNumColumns() + "] != [" + weights.getNumRows() + "," + weights.getNumColumns() + "]"); } CM_COV_Object covobj = new CM_COV_Object(); if (sparse && sparseRows != null) //SPARSE { for (int i = 0; i < rlen; i++) op.fn.execute(covobj, this.quickGetValue(i, 0), that.quickGetValue(i, 0), weights.quickGetValue(i, 0)); } else if (denseBlock != null) //DENSE { //always vectors (see check above) if (!that.sparse && !weights.sparse) { //all dense vectors (default case) if (that.denseBlock != null) for (int i = 0; i < rlen; i++) op.fn.execute(covobj, denseBlock[i], that.denseBlock[i], weights.denseBlock[i]); } else { for (int i = 0; i < rlen; i++) op.fn.execute(covobj, denseBlock[i], that.quickGetValue(i, 0), weights.quickGetValue(i, 0)); } } return covobj; } public MatrixValue sortOperations(MatrixValue weights, MatrixValue result) throws DMLRuntimeException, DMLUnsupportedOperationException { boolean wtflag = (weights != null); MatrixBlock wts = (weights == null ? null : checkType(weights)); checkType(result); if (getNumColumns() != 1) { throw new DMLRuntimeException( "Invalid input dimensions (" + getNumRows() + "x" + getNumColumns() + ") to sort operation."); } if (wts != null && wts.getNumColumns() != 1) { throw new DMLRuntimeException("Invalid weight dimensions (" + wts.getNumRows() + "x" + wts.getNumColumns() + ") to sort operation."); } // prepare result, currently always dense // #rows in temp matrix = 1 + #nnz in the input ( 1 is for the "zero" value) int dim1 = (int) (1 + this.getNonZeros()); if (result == null) result = new MatrixBlock(dim1, 2, false); else result.reset(dim1, 2, false); // Copy the input elements into a temporary array for sorting // First column is data and second column is weights // (since the inputs are vectors, they are likely dense - hence quickget is sufficient) MatrixBlock tdw = new MatrixBlock(dim1, 2, false); double d, w, zero_wt = 0; int ind = 1; if (wtflag) // w/ weights { for (int i = 0; i < rlen; i++) { d = quickGetValue(i, 0); w = wts.quickGetValue(i, 0); if (d != 0) { tdw.quickSetValue(ind, 0, d); tdw.quickSetValue(ind, 1, w); ind++; } else zero_wt += w; } } else //w/o weights { zero_wt = getNumRows() - getNonZeros(); for (int i = 0; i < rlen; i++) { d = quickGetValue(i, 0); if (d != 0) { tdw.quickSetValue(ind, 0, d); tdw.quickSetValue(ind, 1, 1); ind++; } } } tdw.quickSetValue(0, 0, 0.0); tdw.quickSetValue(0, 1, zero_wt); //num zeros in input // Sort td and tw based on values inside td (ascending sort), incl copy into result SortIndex sfn = SortIndex.getSortIndexFnObject(1, false, false); ReorgOperator rop = new ReorgOperator(sfn); LibMatrixReorg.reorg(tdw, (MatrixBlock) result, rop); return result; } public double interQuartileMean() throws DMLRuntimeException { double sum_wt = sumWeightForQuantile(); double q25d = 0.25 * sum_wt; double q75d = 0.75 * sum_wt; int q25i = (int) Math.ceil(q25d); int q75i = (int) Math.ceil(q75d); // skip until (but excluding) q25 int t = 0, i = -1; while (i < getNumRows() && t < q25i) { i++; //System.out.println(" " + i + ": " + quickGetValue(i,0) + "," + quickGetValue(i,1)); t += quickGetValue(i, 1); } // compute the portion of q25 double runningSum = (t - q25d) * quickGetValue(i, 0); // add until (including) q75 while (i < getNumRows() && t < q75i) { i++; t += quickGetValue(i, 1); runningSum += quickGetValue(i, 0) * quickGetValue(i, 1); } // subtract additional portion of q75 runningSum -= (t - q75d) * quickGetValue(i, 0); return runningSum / (sum_wt * 0.5); } /** * Computes the weighted interQuartileMean. * The matrix block ("this" pointer) has two columns, in which the first column * refers to the data and second column denotes corresponding weights. * * @return InterQuartileMean * @throws DMLRuntimeException */ public double interQuartileMeanOLD() throws DMLRuntimeException { double sum_wt = sumWeightForQuantile(); int fromPos = (int) Math.ceil(0.25 * sum_wt); int toPos = (int) Math.ceil(0.75 * sum_wt); int selectRange = toPos - fromPos; // range: (fromPos,toPos] if (selectRange == 0) return 0.0; int index, count = 0; // The first row (0^th row) has value 0. // If it has a non-zero weight i.e., input data has zero values // then "index" must start from 0, otherwise we skip the first row // and start with the next value in the data, which is in the 1st row. if (quickGetValue(0, 1) > 0) index = 0; else index = 1; // keep scanning the weights, until we hit the required position <code>fromPos</code> while (count < fromPos) { count += quickGetValue(index, 1); ++index; } double runningSum; double val; int wt, selectedCount; runningSum = (count - fromPos) * quickGetValue(index - 1, 0); selectedCount = (count - fromPos); while (count <= toPos) { val = quickGetValue(index, 0); wt = (int) quickGetValue(index, 1); runningSum += (val * Math.min(wt, selectRange - selectedCount)); selectedCount += Math.min(wt, selectRange - selectedCount); count += wt; ++index; } //System.out.println(fromPos + ", " + toPos + ": " + count + ", "+ runningSum + ", " + selectedCount); return runningSum / selectedCount; } public MatrixValue pickValues(MatrixValue quantiles, MatrixValue ret) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock qs = checkType(quantiles); if (qs.clen != 1) { throw new DMLRuntimeException("Multiple quantiles can only be computed on a 1D matrix"); } MatrixBlock output = checkType(ret); if (output == null) output = new MatrixBlock(qs.rlen, qs.clen, false); // resulting matrix is mostly likely be dense else output.reset(qs.rlen, qs.clen, false); for (int i = 0; i < qs.rlen; i++) { output.quickSetValue(i, 0, this.pickValue(qs.quickGetValue(i, 0))); } return output; } public double median() throws DMLRuntimeException { double sum_wt = sumWeightForQuantile(); return pickValue(0.5, sum_wt % 2 == 0); } public double pickValue(double quantile) throws DMLRuntimeException { return pickValue(quantile, false); } public double pickValue(double quantile, boolean average) throws DMLRuntimeException { double sum_wt = sumWeightForQuantile(); // do averaging only if it is asked for; and sum_wt is even average = average && (sum_wt % 2 == 0); int pos = (int) Math.ceil(quantile * sum_wt); int t = 0, i = -1; do { i++; t += quickGetValue(i, 1); } while (t < pos && i < getNumRows()); //System.out.println("values: " + quickGetValue(i,0) + "," + quickGetValue(i,1) + " -- " + quickGetValue(i+1,0) + "," + quickGetValue(i+1,1)); if (quickGetValue(i, 1) != 0) { // i^th value is present in the data set, simply return it if (average) { if (pos < t) { return quickGetValue(i, 0); } if (quickGetValue(i + 1, 1) != 0) return (quickGetValue(i, 0) + quickGetValue(i + 1, 0)) / 2; else // (i+1)^th value is 0. So, fetch (i+2)^th value return (quickGetValue(i, 0) + quickGetValue(i + 2, 0)) / 2; } else return quickGetValue(i, 0); } else { // i^th value is not present in the data set. // It can only happen in the case where i^th value is 0.0; and 0.0 is not present in the data set (but introduced by sort). if (i + 1 < getNumRows()) // when 0.0 is not the last element in the sorted order return quickGetValue(i + 1, 0); else // when 0.0 is the last element in the sorted order (input data is all negative) return quickGetValue(i - 1, 0); } } /** * In a given two column matrix, the second column denotes weights. * This function computes the total weight * * @return * @throws DMLRuntimeException */ private double sumWeightForQuantile() throws DMLRuntimeException { double sum_wt = 0; for (int i = 0; i < getNumRows(); i++) sum_wt += quickGetValue(i, 1); if (Math.floor(sum_wt) < sum_wt) { throw new DMLRuntimeException("Unexpected error while computing quantile -- weights must be integers."); } return sum_wt; } /** * * @param m1Index * @param m1Value * @param m2Index * @param m2Value * @param result * @param op * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public MatrixValue aggregateBinaryOperations(MatrixIndexes m1Index, MatrixValue m1Value, MatrixIndexes m2Index, MatrixValue m2Value, MatrixValue result, AggregateBinaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { return aggregateBinaryOperations(m1Value, m2Value, result, op); } /** * */ public MatrixValue aggregateBinaryOperations(MatrixValue m1Value, MatrixValue m2Value, MatrixValue result, AggregateBinaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { //check input types, dimensions, configuration MatrixBlock m1 = checkType(m1Value); MatrixBlock m2 = checkType(m2Value); MatrixBlock ret = checkType(result); if (m1.clen != m2.rlen) { throw new RuntimeException( "Dimensions do not match for matrix multiplication (" + m1.clen + "!=" + m2.rlen + ")."); } if (!(op.binaryFn instanceof Multiply && op.aggOp.increOp.fn instanceof Plus)) { throw new DMLRuntimeException( "Unsupported binary aggregate operation: (" + op.binaryFn + ", " + op.aggOp + ")."); } //setup meta data (dimensions, sparsity) int rl = m1.rlen; int cl = m2.clen; SparsityEstimate sp = estimateSparsityOnAggBinary(m1, m2, op); //create output matrix block if (ret == null) ret = new MatrixBlock(rl, cl, sp.sparse, sp.estimatedNonZeros); else ret.reset(rl, cl, sp.sparse, sp.estimatedNonZeros); //compute matrix multiplication (only supported binary aggregate operation) if (op.getNumThreads() > 1) LibMatrixMult.matrixMult(m1, m2, ret, op.getNumThreads()); else LibMatrixMult.matrixMult(m1, m2, ret); return ret; } /** * * @param m1 * @param m2 * @param m3 * @param op * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public ScalarObject aggregateTernaryOperations(MatrixBlock m1, MatrixBlock m2, MatrixBlock m3, AggregateBinaryOperator op) throws DMLUnsupportedOperationException, DMLRuntimeException { //check input dimensions and operators if (m1.rlen != m2.rlen || m1.clen != m2.clen || (m3 != null && (m2.rlen != m3.rlen || m2.clen != m3.clen))) throw new DMLRuntimeException("Invalid dimensions for aggregate tertiary (" + m1.rlen + "x" + m1.clen + ", " + m2.rlen + "x" + m2.clen + ", " + m3.rlen + "x" + m3.clen + ")."); if (!(op.aggOp.increOp.fn instanceof KahanPlus && op.binaryFn instanceof Multiply)) throw new DMLRuntimeException("Unsupported operator for aggregate tertiary operations."); //execute ternary aggregate function double val = -1; if (op.getNumThreads() > 1) val = LibMatrixAgg.aggregateTernary(m1, m2, m3, op.getNumThreads()); else val = LibMatrixAgg.aggregateTernary(m1, m2, m3); //create output return new DoubleObject(val); } /** * * @param mbLeft * @param mbRight * @param mbOut * @param bOp * @oaram uaggOp * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public MatrixBlock uaggouterchainOperations(MatrixBlock mbLeft, MatrixBlock mbRight, MatrixBlock mbOut, BinaryOperator bOp, AggregateUnaryOperator uaggOp) throws DMLRuntimeException, DMLUnsupportedOperationException { double bv[] = DataConverter.convertToDoubleVector(mbRight); int bvi[] = null; //process instruction if (LibMatrixOuterAgg.isSupportedUaggOp(uaggOp, bOp)) { if ((LibMatrixOuterAgg.isRowIndexMax(uaggOp)) || (LibMatrixOuterAgg.isRowIndexMin(uaggOp))) { bvi = LibMatrixOuterAgg.prepareRowIndices(bv.length, bv, bOp, uaggOp); } else { Arrays.sort(bv); } int iRows = (uaggOp.indexFn instanceof ReduceCol ? mbLeft.getNumRows() : 2); int iCols = (uaggOp.indexFn instanceof ReduceRow ? mbLeft.getNumColumns() : 2); if (mbOut == null) mbOut = new MatrixBlock(iRows, iCols, false); // Output matrix will be dense matrix most of the time. else mbOut.reset(iRows, iCols, false); LibMatrixOuterAgg.aggregateMatrix(mbLeft, mbOut, bv, bvi, bOp, uaggOp); } else throw new DMLRuntimeException("Unsupported operator for unary aggregate operations."); return mbOut; } /** * Invocation from CP instructions. The aggregate is computed on the groups object * against target and weights. * * Notes: * * The computed number of groups is reused for multiple invocations with different target. * * This implementation supports that the target is passed as column or row vector, * in case of row vectors we also use sparse-safe implementations for sparse safe * aggregation operators. * * @param tgt * @param wghts * @param ret * @param op * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixValue groupedAggOperations(MatrixValue tgt, MatrixValue wghts, MatrixValue ret, int ngroups, Operator op) throws DMLRuntimeException, DMLUnsupportedOperationException { //setup input matrices // this <- groups MatrixBlock target = checkType(tgt); MatrixBlock weights = checkType(wghts); //check valid dimensions if (this.getNumColumns() != 1 || (weights != null && weights.getNumColumns() != 1)) throw new DMLRuntimeException( "groupedAggregate can only operate on 1-dimensional column matrices for groups and weights."); if (target.getNumColumns() != 1 && op instanceof CMOperator) throw new DMLRuntimeException( "groupedAggregate can only operate on 1-dimensional column matrices for target (for this aggregation function)."); if (target.getNumColumns() != 1 && target.getNumRows() != 1) throw new DMLRuntimeException( "groupedAggregate can only operate on 1-dimensional column or row matrix for target."); if (this.getNumRows() != Math.max(target.getNumRows(), target.getNumColumns()) || (weights != null && this.getNumRows() != weights.getNumRows())) throw new DMLRuntimeException("groupedAggregate can only operate on matrices with equal dimensions."); // obtain numGroups from instruction, if provided if (ngroups > 0) numGroups = ngroups; // Determine the number of groups if (numGroups <= 0) //reuse if available { double min = this.min(); double max = this.max(); if (min <= 0) throw new DMLRuntimeException( "Invalid value (" + min + ") encountered in 'groups' while computing groupedAggregate"); if (max <= 0) throw new DMLRuntimeException( "Invalid value (" + max + ") encountered in 'groups' while computing groupedAggregate."); numGroups = (int) max; } // Allocate result matrix MatrixBlock result = checkType(ret); boolean result_sparsity = estimateSparsityOnGroupedAgg(rlen, numGroups); if (result == null) result = new MatrixBlock(numGroups, 1, result_sparsity); else result.reset(numGroups, 1, result_sparsity); // Compute the result double w = 1; // default weight //CM operator for count, mean, variance //note: current support only for column vectors if (op instanceof CMOperator) { // initialize required objects for storing the result of CM operations CM cmFn = CM.getCMFnObject(((CMOperator) op).getAggOpType()); CM_COV_Object[] cmValues = new CM_COV_Object[numGroups]; for (int i = 0; i < numGroups; i++) cmValues[i] = new CM_COV_Object(); for (int i = 0; i < this.getNumRows(); i++) { int g = (int) this.quickGetValue(i, 0); if (g > numGroups) continue; double d = target.quickGetValue(i, 0); if (weights != null) w = weights.quickGetValue(i, 0); // cmValues is 0-indexed, whereas range of values for g = [1,numGroups] cmFn.execute(cmValues[g - 1], d, w); } // extract the required value from each CM_COV_Object for (int i = 0; i < numGroups; i++) // result is 0-indexed, so is cmValues result.quickSetValue(i, 0, cmValues[i].getRequiredResult(op)); } //Aggregate operator for sum (via kahan sum) //note: support for row/column vectors and dense/sparse else if (op instanceof AggregateOperator) { //the only aggregate operator that is supported here is sum, //furthermore, we always use KahanPlus and hence aggop.correctionExists is true AggregateOperator aggop = (AggregateOperator) op; //default case for aggregate(sum) groupedAggregateKahanPlus(target, weights, result, aggop); } else throw new DMLRuntimeException( "Invalid operator (" + op + ") encountered while processing groupedAggregate."); return result; } /** * This is a specific implementation for aggregate(fn="sum"), where we use KahanPlus for numerical * stability. In contrast to other functions of aggregate, this implementation supports row and column * vectors for target and exploits sparse representations since KahanPlus is sparse-safe. * * @param target * @param weights * @param op * @throws DMLRuntimeException */ private void groupedAggregateKahanPlus(MatrixBlock target, MatrixBlock weights, MatrixBlock result, AggregateOperator aggop) throws DMLRuntimeException { boolean rowVector = target.getNumColumns() > 1; double w = 1; //default weight //skip empty blocks (sparse-safe operation) if (target.isEmptyBlock(false)) return; //init group buffers KahanObject[] buffer = new KahanObject[numGroups]; for (int i = 0; i < numGroups; i++) buffer[i] = new KahanObject(aggop.initialValue, 0); if (rowVector) //target is rowvector { if (target.sparse) //SPARSE target { if (target.sparseRows[0] != null) { int len = target.sparseRows[0].size(); int[] aix = target.sparseRows[0].getIndexContainer(); double[] avals = target.sparseRows[0].getValueContainer(); for (int j = 0; j < len; j++) //for each nnz { int g = (int) this.quickGetValue(aix[j], 0); if (g > numGroups) continue; if (weights != null) w = weights.quickGetValue(aix[j], 0); aggop.increOp.fn.execute(buffer[g - 1], avals[j] * w); } } } else //DENSE target { for (int i = 0; i < target.getNumColumns(); i++) { double d = target.denseBlock[i]; if (d != 0) //sparse-safe { int g = (int) this.quickGetValue(i, 0); if (g > numGroups) continue; if (weights != null) w = weights.quickGetValue(i, 0); // buffer is 0-indexed, whereas range of values for g = [1,numGroups] aggop.increOp.fn.execute(buffer[g - 1], d * w); } } } } else //column vector (always dense, but works for sparse as well) { for (int i = 0; i < this.getNumRows(); i++) { double d = target.quickGetValue(i, 0); if (d != 0) //sparse-safe { int g = (int) this.quickGetValue(i, 0); if (g > numGroups) continue; if (weights != null) w = weights.quickGetValue(i, 0); // buffer is 0-indexed, whereas range of values for g = [1,numGroups] aggop.increOp.fn.execute(buffer[g - 1], d * w); } } } // extract the results from group buffers for (int i = 0; i < numGroups; i++) result.quickSetValue(i, 0, buffer[i]._sum); } /** * * @param ret * @param rows * @param select * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock removeEmptyOperations(MatrixBlock ret, boolean rows, MatrixBlock select) throws DMLRuntimeException, DMLUnsupportedOperationException { MatrixBlock result = checkType(ret); return LibMatrixReorg.rmempty(this, result, rows, select); } /** * * @param ret * @param rows * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock removeEmptyOperations(MatrixBlock ret, boolean rows) throws DMLRuntimeException, DMLUnsupportedOperationException { return removeEmptyOperations(ret, rows, null); } /** * * @param ret * @param max * @param rows * @param cast * @param ignore * @return * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ public MatrixBlock rexpandOperations(MatrixBlock ret, double max, boolean rows, boolean cast, boolean ignore) throws DMLRuntimeException, DMLUnsupportedOperationException { MatrixBlock result = checkType(ret); return LibMatrixReorg.rexpand(this, result, max, rows, cast, ignore); } @Override public MatrixValue replaceOperations(MatrixValue result, double pattern, double replacement) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock ret = checkType(result); examSparsity(); //ensure its in the right format ret.reset(rlen, clen, sparse); if (nonZeros == 0 && pattern != 0) return ret; //early abort boolean NaNpattern = Double.isNaN(pattern); if (sparse) //SPARSE { if (pattern != 0d) //SPARSE <- SPARSE (sparse-safe) { ret.allocateSparseRowsBlock(); SparseRow[] a = sparseRows; SparseRow[] c = ret.sparseRows; for (int i = 0; i < rlen; i++) { SparseRow arow = a[i]; if (arow != null && !arow.isEmpty()) { SparseRow crow = new SparseRow(arow.size()); int alen = arow.size(); int[] aix = arow.getIndexContainer(); double[] avals = arow.getValueContainer(); for (int j = 0; j < alen; j++) { double val = avals[j]; if (val == pattern || (NaNpattern && Double.isNaN(val))) crow.append(aix[j], replacement); else crow.append(aix[j], val); } c[i] = crow; } } } else //DENSE <- SPARSE { ret.sparse = false; ret.allocateDenseBlock(); SparseRow[] a = sparseRows; double[] c = ret.denseBlock; //initialize with replacement (since all 0 values, see SPARSITY_TURN_POINT) Arrays.fill(c, replacement); //overwrite with existing values (via scatter) if (a != null) //check for empty matrix for (int i = 0, cix = 0; i < rlen; i++, cix += clen) { SparseRow arow = a[i]; if (arow != null && !arow.isEmpty()) { int alen = arow.size(); int[] aix = arow.getIndexContainer(); double[] avals = arow.getValueContainer(); for (int j = 0; j < alen; j++) if (avals[j] != 0) c[cix + aix[j]] = avals[j]; } } } } else //DENSE <- DENSE { int mn = ret.rlen * ret.clen; ret.allocateDenseBlock(); double[] a = denseBlock; double[] c = ret.denseBlock; for (int i = 0; i < mn; i++) { double val = a[i]; if (val == pattern || (NaNpattern && Double.isNaN(val))) c[i] = replacement; else c[i] = a[i]; } } ret.recomputeNonZeros(); ret.examSparsity(); return ret; } /** * D = ctable(A,v2,W) * this <- A; scalarThat <- v2; that2 <- W; result <- D * * (i1,j1,v1) from input1 (this) * (v2) from sclar_input2 (scalarThat) * (i3,j3,w) from input3 (that2) * @throws DMLRuntimeException * @throws DMLUnsupportedOperationException */ @Override public void ternaryOperations(Operator op, double scalarThat, MatrixValue that2Val, CTableMap resultMap, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock that2 = checkType(that2Val); CTable ctable = CTable.getCTableFnObject(); double v2 = scalarThat; //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) if (resultBlock == null) { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); double w = that2.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultMap); } } else { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); double w = that2.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultBlock); } resultBlock.recomputeNonZeros(); } } /** * D = ctable(A,v2,w) * this <- A; scalar_that <- v2; scalar_that2 <- w; result <- D * * (i1,j1,v1) from input1 (this) * (v2) from sclar_input2 (scalarThat) * (w) from scalar_input3 (scalarThat2) */ @Override public void ternaryOperations(Operator op, double scalarThat, double scalarThat2, CTableMap resultMap, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { CTable ctable = CTable.getCTableFnObject(); double v2 = scalarThat; double w = scalarThat2; //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) if (resultBlock == null) { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultMap); } } else { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultBlock); } resultBlock.recomputeNonZeros(); } } /** * Specific ctable case of ctable(seq(...),X), where X is the only * matrix input. The 'left' input parameter specifies if the seq appeared * on the left, otherwise it appeared on the right. * */ @Override public void ternaryOperations(Operator op, MatrixIndexes ix1, double scalarThat, boolean left, int brlen, CTableMap resultMap, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { CTable ctable = CTable.getCTableFnObject(); double w = scalarThat; int offset = (int) ((ix1.getRowIndex() - 1) * brlen); //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) if (resultBlock == null) { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); if (left) ctable.execute(offset + i + 1, v1, w, false, resultMap); else ctable.execute(v1, offset + i + 1, w, false, resultMap); } } else { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); if (left) ctable.execute(offset + i + 1, v1, w, false, resultBlock); else ctable.execute(v1, offset + i + 1, w, false, resultBlock); } resultBlock.recomputeNonZeros(); } } /** * D = ctable(A,B,w) * this <- A; that <- B; scalar_that2 <- w; result <- D * * (i1,j1,v1) from input1 (this) * (i1,j1,v2) from input2 (that) * (w) from scalar_input3 (scalarThat2) * * NOTE: This method supports both vectors and matrices. In case of matrices and ignoreZeros=true * we can also use a sparse-safe implementation */ @Override public void ternaryOperations(Operator op, MatrixValue thatVal, double scalarThat2, boolean ignoreZeros, CTableMap resultMap, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { //setup ctable computation MatrixBlock that = checkType(thatVal); CTable ctable = CTable.getCTableFnObject(); double w = scalarThat2; if (ignoreZeros //SPARSE-SAFE & SPARSE INPUTS && this.sparse && that.sparse) { //note: only used if both inputs have aligned zeros, which //allows us to infer that the nnz both inputs are equivalent //early abort on empty blocks possible if (this.isEmptyBlock(false) && that.isEmptyBlock(false)) return; SparseRow[] a = this.sparseRows; SparseRow[] b = that.sparseRows; for (int i = 0; i < rlen; i++) { SparseRow arow = a[i]; SparseRow brow = b[i]; if (arow != null && !arow.isEmpty()) { int alen = arow.size(); double[] avals = arow.getValueContainer(); double[] bvals = brow.getValueContainer(); if (resultBlock == null) { for (int j = 0; j < alen; j++) ctable.execute(avals[j], bvals[j], w, ignoreZeros, resultMap); } else { for (int j = 0; j < alen; j++) ctable.execute(avals[j], bvals[j], w, ignoreZeros, resultBlock); } } } } else //SPARSE-UNSAFE | GENERIC INPUTS { //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); double v2 = that.quickGetValue(i, j); if (resultBlock == null) ctable.execute(v1, v2, w, ignoreZeros, resultMap); else ctable.execute(v1, v2, w, ignoreZeros, resultBlock); } } //maintain nnz (if necessary) if (resultBlock != null) resultBlock.recomputeNonZeros(); } /** * D = ctable(seq,A,w) * this <- seq; thatMatrix <- A; thatScalar <- w; result <- D * * (i1,j1,v1) from input1 (this) * (i1,j1,v2) from input2 (that) * (w) from scalar_input3 (scalarThat2) */ public void ternaryOperations(Operator op, MatrixValue thatMatrix, double thatScalar, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock that = checkType(thatMatrix); CTable ctable = CTable.getCTableFnObject(); double w = thatScalar; //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) //resultBlock guaranteed to be allocated for ctableexpand //each row in resultBlock will be allocated and will contain exactly one value int maxCol = 0; for (int i = 0; i < rlen; i++) { double v2 = that.quickGetValue(i, 0); maxCol = ctable.execute(i + 1, v2, w, maxCol, resultBlock); } //update meta data (initially unknown number of columns) //note: nnz maintained in ctable (via quickset) resultBlock.clen = maxCol; } /** * D = ctable(A,B,W) * this <- A; that <- B; that2 <- W; result <- D * * (i1,j1,v1) from input1 (this) * (i1,j1,v2) from input2 (that) * (i1,j1,w) from input3 (that2) */ public void ternaryOperations(Operator op, MatrixValue thatVal, MatrixValue that2Val, CTableMap resultMap) throws DMLUnsupportedOperationException, DMLRuntimeException { ternaryOperations(op, thatVal, that2Val, resultMap, null); } @Override public void ternaryOperations(Operator op, MatrixValue thatVal, MatrixValue that2Val, CTableMap resultMap, MatrixBlock resultBlock) throws DMLUnsupportedOperationException, DMLRuntimeException { MatrixBlock that = checkType(thatVal); MatrixBlock that2 = checkType(that2Val); CTable ctable = CTable.getCTableFnObject(); //sparse-unsafe ctable execution //(because input values of 0 are invalid and have to result in errors) if (resultBlock == null) { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); double v2 = that.quickGetValue(i, j); double w = that2.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultMap); } } else { for (int i = 0; i < rlen; i++) for (int j = 0; j < clen; j++) { double v1 = this.quickGetValue(i, j); double v2 = that.quickGetValue(i, j); double w = that2.quickGetValue(i, j); ctable.execute(v1, v2, w, false, resultBlock); } resultBlock.recomputeNonZeros(); } } @Override public MatrixValue quaternaryOperations(QuaternaryOperator qop, MatrixValue um, MatrixValue vm, MatrixValue wm, MatrixValue out) throws DMLUnsupportedOperationException, DMLRuntimeException { return quaternaryOperations(qop, um, vm, wm, out, 1); } /** * * @param op * @param um * @param vm * @param wm * @param out * @param wt * @param k * @return * @throws DMLUnsupportedOperationException * @throws DMLRuntimeException */ public MatrixValue quaternaryOperations(QuaternaryOperator qop, MatrixValue um, MatrixValue vm, MatrixValue wm, MatrixValue out, int k) throws DMLUnsupportedOperationException, DMLRuntimeException { //check input dimensions if (getNumRows() != um.getNumRows()) throw new DMLRuntimeException( "Dimension mismatch rows on quaternary operation: " + getNumRows() + "!=" + um.getNumRows()); if (getNumColumns() != vm.getNumRows()) throw new DMLRuntimeException( "Dimension mismatch columns quaternary operation: " + getNumColumns() + "!=" + vm.getNumRows()); //check input data types MatrixBlock X = this; MatrixBlock U = checkType(um); MatrixBlock V = checkType(vm); MatrixBlock R = checkType(out); //prepare intermediates and output if (qop.wtype1 != null || qop.wtype4 != null) R.reset(1, 1, false); else if (qop.wtype2 != null) R.reset(rlen, clen, sparse); else if (qop.wtype3 != null) { MatrixCharacteristics mc = qop.wtype3.computeOutputCharacteristics(X.rlen, X.clen, U.clen); R.reset((int) mc.getRows(), (int) mc.getCols(), qop.wtype3.isBasic() ? X.isInSparseFormat() : false); } //core block operation if (qop.wtype1 != null) { //wsloss MatrixBlock W = qop.wtype1.hasFourInputs() ? checkType(wm) : null; if (k > 1) LibMatrixMult.matrixMultWSLoss(X, U, V, W, R, qop.wtype1, k); else LibMatrixMult.matrixMultWSLoss(X, U, V, W, R, qop.wtype1); } else if (qop.wtype2 != null) { //wsigmoid if (k > 1) LibMatrixMult.matrixMultWSigmoid(X, U, V, R, qop.wtype2, k); else LibMatrixMult.matrixMultWSigmoid(X, U, V, R, qop.wtype2); } else if (qop.wtype3 != null) { //wdivmm //note: for wdivmm-minus X and W interchanged because W always present if (k > 1) LibMatrixMult.matrixMultWDivMM(X, U, V, R, qop.wtype3, k); else LibMatrixMult.matrixMultWDivMM(X, U, V, R, qop.wtype3); } else if (qop.wtype4 != null) { //wcemm if (k > 1) LibMatrixMult.matrixMultWCeMM(X, U, V, R, qop.wtype4, k); else LibMatrixMult.matrixMultWCeMM(X, U, V, R, qop.wtype4); } return R; } //////// // Data Generation Methods // (rand, sequence) /** * Function to generate the random matrix with specified dimensions (block sizes are not specified). * * @param rows * @param cols * @param sparsity * @param min * @param max * @param pdf * @param seed * @return * @throws DMLRuntimeException */ public static MatrixBlock randOperations(int rows, int cols, double sparsity, double min, double max, String pdf, long seed) throws DMLRuntimeException { return randOperations(rows, cols, sparsity, min, max, pdf, seed, 1); } /** * Function to generate the random matrix with specified dimensions (block sizes are not specified). * * @param rows * @param cols * @param sparsity * @param min * @param max * @param pdf * @param seed * @return * @throws DMLRuntimeException */ public static MatrixBlock randOperations(int rows, int cols, double sparsity, double min, double max, String pdf, long seed, int k) throws DMLRuntimeException { DMLConfig conf = ConfigurationManager.getConfig(); int blocksize = (conf != null) ? ConfigurationManager.getConfig().getIntValue(DMLConfig.DEFAULT_BLOCK_SIZE) : DMLTranslator.DMLBlockSize; RandomMatrixGenerator rgen = new RandomMatrixGenerator(pdf, rows, cols, blocksize, blocksize, sparsity, min, max); if (k > 1) return randOperations(rgen, seed, k); else return randOperations(rgen, seed); } /** * Function to generate the random matrix with specified dimensions and block dimensions. * @param rgen * @param seed * @return * @throws DMLRuntimeException */ public static MatrixBlock randOperations(RandomMatrixGenerator rgen, long seed) throws DMLRuntimeException { return randOperations(rgen, seed, 1); } /** * Function to generate the random matrix with specified dimensions and block dimensions. * @param rows * @param cols * @param rowsInBlock * @param colsInBlock * @param sparsity * @param min * @param max * @param pdf * @param seed * @return * @throws DMLRuntimeException */ public static MatrixBlock randOperations(RandomMatrixGenerator rgen, long seed, int k) throws DMLRuntimeException { MatrixBlock out = new MatrixBlock(); Well1024a bigrand = null; long[] nnzInBlock = null; //setup seeds and nnz per block if (!LibMatrixDatagen.isShortcutRandOperation(rgen._min, rgen._max, rgen._sparsity, rgen._pdf)) { bigrand = LibMatrixDatagen.setupSeedsForRand(seed); nnzInBlock = LibMatrixDatagen.computeNNZperBlock(rgen._rows, rgen._cols, rgen._rowsPerBlock, rgen._colsPerBlock, rgen._sparsity); } //generate rand data if (k > 1) out.randOperationsInPlace(rgen, nnzInBlock, bigrand, -1, k); else out.randOperationsInPlace(rgen, nnzInBlock, bigrand, -1); return out; } /** * Function to generate a matrix of random numbers. This is invoked both * from CP as well as from MR. In case of CP, it generates an entire matrix * block-by-block. A <code>bigrand</code> is passed so that block-level * seeds are generated internally. In case of MR, it generates a single * block for given block-level seed <code>bSeed</code>. * * When pdf="uniform", cell values are drawn from uniform distribution in * range <code>[min,max]</code>. * * When pdf="normal", cell values are drawn from standard normal * distribution N(0,1). The range of generated values will always be * (-Inf,+Inf). * * @param rgen * @param nnzInBlock * @param bigrand * @param bSeed * @return * @throws DMLRuntimeException */ public MatrixBlock randOperationsInPlace(RandomMatrixGenerator rgen, long[] nnzInBlock, Well1024a bigrand, long bSeed) throws DMLRuntimeException { LibMatrixDatagen.generateRandomMatrix(this, rgen, nnzInBlock, bigrand, bSeed); return this; } /** * Function to generate a matrix of random numbers. This is invoked both * from CP as well as from MR. In case of CP, it generates an entire matrix * block-by-block. A <code>bigrand</code> is passed so that block-level * seeds are generated internally. In case of MR, it generates a single * block for given block-level seed <code>bSeed</code>. * * When pdf="uniform", cell values are drawn from uniform distribution in * range <code>[min,max]</code>. * * When pdf="normal", cell values are drawn from standard normal * distribution N(0,1). The range of generated values will always be * (-Inf,+Inf). * * @param rgen * @param nnzInBlock * @param bigrand * @param bSeed * @param k * @return * @throws DMLRuntimeException */ public MatrixBlock randOperationsInPlace(RandomMatrixGenerator rgen, long[] nnzInBlock, Well1024a bigrand, long bSeed, int k) throws DMLRuntimeException { LibMatrixDatagen.generateRandomMatrix(this, rgen, nnzInBlock, bigrand, bSeed, k); return this; } /** * Method to generate a sequence according to the given parameters. The * generated sequence is always in dense format. * * Both end points specified <code>from</code> and <code>to</code> must be * included in the generated sequence i.e., [from,to] both inclusive. Note * that, <code>to</code> is included only if (to-from) is perfectly * divisible by <code>incr</code>. * * For example, seq(0,1,0.5) generates (0.0 0.5 1.0) * whereas seq(0,1,0.6) generates (0.0 0.6) but not (0.0 0.6 1.0) * * @param from * @param to * @param incr * @return * @throws DMLRuntimeException */ public static MatrixBlock seqOperations(double from, double to, double incr) throws DMLRuntimeException { MatrixBlock out = new MatrixBlock(); LibMatrixDatagen.generateSequence(out, from, to, incr); return out; } /** * * @param from * @param to * @param incr * @return * @throws DMLRuntimeException */ public MatrixBlock seqOperationsInPlace(double from, double to, double incr) throws DMLRuntimeException { LibMatrixDatagen.generateSequence(this, from, to, incr); return this; } /** * * @param range * @param size * @param replace * @return * @throws DMLRuntimeException */ public static MatrixBlock sampleOperations(long range, int size, boolean replace, long seed) throws DMLRuntimeException { MatrixBlock out = new MatrixBlock(); LibMatrixDatagen.generateSample(out, range, size, replace, seed); return out; } //////// // Misc methods private static MatrixBlock checkType(MatrixValue block) throws RuntimeException { if (block != null && !(block instanceof MatrixBlock)) throw new RuntimeException("Unsupported matrix value: " + block.getClass().getSimpleName()); return (MatrixBlock) block; } public void print() { System.out.println("sparse = " + sparse); if (!sparse) System.out.println("nonzeros = " + nonZeros); for (int i = 0; i < rlen; i++) { for (int j = 0; j < clen; j++) { System.out.print(quickGetValue(i, j) + "\t"); } System.out.println(); } } @Override public int compareTo(Object arg0) { throw new RuntimeException("CompareTo should never be called for matrix blocks."); } @Override public boolean equals(Object arg0) { throw new RuntimeException("equals should never be called for matrix blocks."); } @Override public int hashCode() { throw new RuntimeException("HashCode should never be called for matrix blocks."); } @Override public String toString() { StringBuilder sb = new StringBuilder(); sb.append("sparse? = "); sb.append(sparse); sb.append("\n"); sb.append("nonzeros = "); sb.append(nonZeros); sb.append("\n"); sb.append("size: "); sb.append(rlen); sb.append(" X "); sb.append(clen); sb.append("\n"); if (sparse) { int len = 0; if (sparseRows != null) len = Math.min(rlen, sparseRows.length); int i = 0; for (; i < len; i++) { sb.append("row +"); sb.append(i); sb.append(": "); sb.append(sparseRows[i]); sb.append("\n"); } for (; i < rlen; i++) { sb.append("row +"); sb.append(i); sb.append(": null\n"); } } else { if (denseBlock != null) { for (int i = 0, ix = 0; i < rlen; i++, ix += clen) { for (int j = 0; j < clen; j++) { sb.append(this.denseBlock[ix + j]); sb.append("\t"); } sb.append("\n"); } } } return sb.toString(); } /////////////////////////// // Helper classes public static class SparsityEstimate { public long estimatedNonZeros = 0; public boolean sparse = false; public SparsityEstimate(boolean sps, long nnzs) { sparse = sps; estimatedNonZeros = nnzs; } public SparsityEstimate() { } } }