Java Code Examples for cern.colt.matrix.DoubleMatrix2D#viewPart()

/**
 */
public static void doubleTest() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
System.out.println(master);
master.assign(1); // set all cells to 1
System.out.println("\n"+master);
master.viewPart(2,1,2,3).assign(2); // set [2,1] .. [3,3] to 2
System.out.println("\n"+master);

DoubleMatrix2D copyPart = master.viewPart(2,1,2,3).copy();
copyPart.assign(3); // modify an independent copy
copyPart.set(0,0,4); 
System.out.println("\n"+copyPart); // has changed
System.out.println("\n"+master); // master has not changed

DoubleMatrix2D view1 = master.viewPart(0,3,4,2); // [0,3] .. [3,4]
DoubleMatrix2D view2 = view1.viewPart(0,0,4,1); // a view from a view 
System.out.println("\n"+view1);
System.out.println("\n"+view2);
}
 

/**
 */
public static void doubleTest7() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = Factory2D.ascending(rows,columns);
//master.assign(1); // set all cells to 1
System.out.println("\n"+master);
//master.viewPart(2,0,2,3).assign(2); // set [2,1] .. [3,3] to 2
//System.out.println("\n"+master);

int[] rowIndexes = {0,1,3,0};
int[] columnIndexes = {0,2};
DoubleMatrix2D view1 = master.viewSelection(rowIndexes,columnIndexes);
System.out.println("view1="+view1);
DoubleMatrix2D view2 = view1.viewPart(0,0,2,2);
System.out.println("view2="+view2);

view2.assign(-1);
System.out.println("master replaced"+master);
System.out.println("flip1 replaced"+view1);
System.out.println("flip2 replaced"+view2);

}
 

/**
 * Copies the rows of the indicated columns into a new sub matrix.
 * <tt>sub[0..rowTo-rowFrom,0..columnIndexes.length-1] = A[rowFrom..rowTo,columnIndexes(:)]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowFrom the index of the first row to copy (inclusive).
 * @param   rowTo the index of the last row to copy (inclusive).
 * @param   columnIndexes the indexes of the columns to copy. May be unsorted.
 * @return  a new sub matrix; with <tt>sub.rows()==rowTo-rowFrom+1; sub.columns()==columnIndexes.length</tt>.
 * @throws	IndexOutOfBoundsException if <tt>rowFrom<0 || rowTo-rowFrom+1<0 || rowTo+1>matrix.rows() || for any col=columnIndexes[i]: col < 0 || col >= matrix.columns()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int rowFrom, int rowTo, int[] columnIndexes) {
	if (rowTo-rowFrom >= A.rows()) throw new IndexOutOfBoundsException("Too many rows");
	int height = rowTo-rowFrom+1;
	int columns = A.columns();
	A = A.viewPart(rowFrom,0,height,columns);
	DoubleMatrix2D sub = A.like(height, columnIndexes.length);
	
	for (int c = columnIndexes.length; --c >= 0; ) {
		int column = columnIndexes[c];
		if (column < 0 || column >= columns)
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewColumn(c).assign(A.viewColumn(column));
	}
	return sub;
}
 

/**
 * Copies the rows of the indicated columns into a new sub matrix.
 * <tt>sub[0..rowTo-rowFrom,0..columnIndexes.length-1] = A[rowFrom..rowTo,columnIndexes(:)]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowFrom the index of the first row to copy (inclusive).
 * @param   rowTo the index of the last row to copy (inclusive).
 * @param   columnIndexes the indexes of the columns to copy. May be unsorted.
 * @return  a new sub matrix; with <tt>sub.rows()==rowTo-rowFrom+1; sub.columns()==columnIndexes.length</tt>.
 * @throws	IndexOutOfBoundsException if <tt>rowFrom<0 || rowTo-rowFrom+1<0 || rowTo+1>matrix.rows() || for any col=columnIndexes[i]: col < 0 || col >= matrix.columns()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int rowFrom, int rowTo, int[] columnIndexes) {
	if (rowTo-rowFrom >= A.rows()) throw new IndexOutOfBoundsException("Too many rows");
	int height = rowTo-rowFrom+1;
	int columns = A.columns();
	A = A.viewPart(rowFrom,0,height,columns);
	DoubleMatrix2D sub = A.like(height, columnIndexes.length);
	
	for (int c = columnIndexes.length; --c >= 0; ) {
		int column = columnIndexes[c];
		if (column < 0 || column >= columns)
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewColumn(c).assign(A.viewColumn(column));
	}
	return sub;
}
 

/**
 */
public static void doubleTest7() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = Factory2D.ascending(rows,columns);
//master.assign(1); // set all cells to 1
System.out.println("\n"+master);
//master.viewPart(2,0,2,3).assign(2); // set [2,1] .. [3,3] to 2
//System.out.println("\n"+master);

int[] rowIndexes = {0,1,3,0};
int[] columnIndexes = {0,2};
DoubleMatrix2D view1 = master.viewSelection(rowIndexes,columnIndexes);
System.out.println("view1="+view1);
DoubleMatrix2D view2 = view1.viewPart(0,0,2,2);
System.out.println("view2="+view2);

view2.assign(-1);
System.out.println("master replaced"+master);
System.out.println("flip1 replaced"+view1);
System.out.println("flip2 replaced"+view2);

}
 

/**
 */
public static void doubleTest() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
System.out.println(master);
master.assign(1); // set all cells to 1
System.out.println("\n"+master);
master.viewPart(2,1,2,3).assign(2); // set [2,1] .. [3,3] to 2
System.out.println("\n"+master);

DoubleMatrix2D copyPart = master.viewPart(2,1,2,3).copy();
copyPart.assign(3); // modify an independent copy
copyPart.set(0,0,4); 
System.out.println("\n"+copyPart); // has changed
System.out.println("\n"+master); // master has not changed

DoubleMatrix2D view1 = master.viewPart(0,3,4,2); // [0,3] .. [3,4]
DoubleMatrix2D view2 = view1.viewPart(0,0,4,1); // a view from a view 
System.out.println("\n"+view1);
System.out.println("\n"+view2);
}
 

/**
 */
public static void doubleTest10() {
int rows = 6;
int columns = 7; // make a 4*5 matrix
//DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
DoubleMatrix2D master = Factory2D.ascending(rows,columns);
//Basic.ascending(master);
//master.assign(1); // set all cells to 1
Transform.mult(master,Math.sin(0.3));
System.out.println("\n"+master);
//master.viewPart(2,0,2,3).assign(2); // set [2,1] .. [3,3] to 2
//System.out.println("\n"+master);

int[] rowIndexes = {0,1,2,3};
int[] columnIndexes = {0,1,2,3};

int[] rowIndexes2 = {3,0,3};
int[] columnIndexes2 = {3,0,3};
DoubleMatrix2D view1 = master.viewPart(1,1,4,5).viewSelection(rowIndexes,columnIndexes);
System.out.println("\nview1="+view1);
DoubleMatrix2D view9 = view1.viewStrides(2,2).viewStrides(2,1);
System.out.println("\nview9="+view9);
view1 = view1.viewSelection(rowIndexes2,columnIndexes2);
System.out.println("\nview1="+view1);
DoubleMatrix2D view2 = view1.viewPart(1,1,2,2);
System.out.println("\nview2="+view2);
DoubleMatrix2D view3 = view2.viewRowFlip();
System.out.println("\nview3="+view3);
view3.assign(Factory2D.ascending(view3.rows(),view3.columns()));
//Basic.ascending(view3);
System.out.println("\nview3="+view3);

//view2.assign(-1);
System.out.println("\nmaster replaced"+master);
System.out.println("\nview1 replaced"+view1);
System.out.println("\nview2 replaced"+view2);
System.out.println("\nview3 replaced"+view3);

}
 

/** 
Least squares solution of <tt>A*X = B</tt>; <tt>returns X</tt>.
@param B    A matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> that minimizes the two norm of <tt>Q*R*X - B</tt>.
@exception  IllegalArgumentException  if <tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if <tt>!this.hasFullRank()</tt> (<tt>A</tt> is rank deficient).
*/
public DoubleMatrix2D solve(DoubleMatrix2D B) {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	if (B.rows() != m) {
		throw new IllegalArgumentException("Matrix row dimensions must agree.");
	}
	if (!this.hasFullRank()) {
		throw new IllegalArgumentException("Matrix is rank deficient.");
	}
	
	// Copy right hand side
	int nx = B.columns();
	DoubleMatrix2D X = B.copy();
	
	// Compute Y = transpose(Q)*B
	for (int k = 0; k < n; k++) {
		for (int j = 0; j < nx; j++) {
			double s = 0.0; 
			for (int i = k; i < m; i++) {
				s += QR.getQuick(i,k)*X.getQuick(i,j);
			}
			s = -s / QR.getQuick(k,k);
			for (int i = k; i < m; i++) {
				X.setQuick(i,j, X.getQuick(i,j) + s*QR.getQuick(i,k));
			}
		}
	}
	// Solve R*X = Y;
	for (int k = n-1; k >= 0; k--) {
		for (int j = 0; j < nx; j++) {
			X.setQuick(k,j, X.getQuick(k,j) / Rdiag.getQuick(k));
		}
		for (int i = 0; i < k; i++) {
			for (int j = 0; j < nx; j++) {
				X.setQuick(i,j, X.getQuick(i,j) - X.getQuick(k,j)*QR.getQuick(i,k));
			}
		}
	}
	return X.viewPart(0,0,n,nx);
}
 

/**
 * Copies the columns of the indicated rows into a new sub matrix.
 * <tt>sub[0..rowIndexes.length-1,0..columnTo-columnFrom] = A[rowIndexes(:),columnFrom..columnTo]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowIndexes the indexes of the rows to copy. May be unsorted.
 * @param   columnFrom the index of the first column to copy (inclusive).
 * @param   columnTo the index of the last column to copy (inclusive).
 * @return  a new sub matrix; with <tt>sub.rows()==rowIndexes.length; sub.columns()==columnTo-columnFrom+1</tt>.
 * @throws	IndexOutOfBoundsException if <tt>columnFrom<0 || columnTo-columnFrom+1<0 || columnTo+1>matrix.columns() || for any row=rowIndexes[i]: row < 0 || row >= matrix.rows()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int[] rowIndexes, int columnFrom, int columnTo) {
	int width = columnTo-columnFrom+1;
	int rows = A.rows();
	A = A.viewPart(0,columnFrom,rows,width);
	DoubleMatrix2D sub = A.like(rowIndexes.length, width);
	
	for (int r = rowIndexes.length; --r >= 0; ) {
		int row = rowIndexes[r];
		if (row < 0 || row >= rows) 
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewRow(r).assign(A.viewRow(row));
	}
	return sub;
}
 

/**
 */
public static void doubleTest10() {
int rows = 6;
int columns = 7; // make a 4*5 matrix
//DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
DoubleMatrix2D master = Factory2D.ascending(rows,columns);
//Basic.ascending(master);
//master.assign(1); // set all cells to 1
Transform.mult(master,Math.sin(0.3));
System.out.println("\n"+master);
//master.viewPart(2,0,2,3).assign(2); // set [2,1] .. [3,3] to 2
//System.out.println("\n"+master);

int[] rowIndexes = {0,1,2,3};
int[] columnIndexes = {0,1,2,3};

int[] rowIndexes2 = {3,0,3};
int[] columnIndexes2 = {3,0,3};
DoubleMatrix2D view1 = master.viewPart(1,1,4,5).viewSelection(rowIndexes,columnIndexes);
System.out.println("\nview1="+view1);
DoubleMatrix2D view9 = view1.viewStrides(2,2).viewStrides(2,1);
System.out.println("\nview9="+view9);
view1 = view1.viewSelection(rowIndexes2,columnIndexes2);
System.out.println("\nview1="+view1);
DoubleMatrix2D view2 = view1.viewPart(1,1,2,2);
System.out.println("\nview2="+view2);
DoubleMatrix2D view3 = view2.viewRowFlip();
System.out.println("\nview3="+view3);
view3.assign(Factory2D.ascending(view3.rows(),view3.columns()));
//Basic.ascending(view3);
System.out.println("\nview3="+view3);

//view2.assign(-1);
System.out.println("\nmaster replaced"+master);
System.out.println("\nview1 replaced"+view1);
System.out.println("\nview2 replaced"+view2);
System.out.println("\nview3 replaced"+view3);

}
 

/** 
Least squares solution of <tt>A*X = B</tt>; <tt>returns X</tt>.
@param B    A matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> that minimizes the two norm of <tt>Q*R*X - B</tt>.
@exception  IllegalArgumentException  if <tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if <tt>!this.hasFullRank()</tt> (<tt>A</tt> is rank deficient).
*/
public DoubleMatrix2D solve(DoubleMatrix2D B) {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	if (B.rows() != m) {
		throw new IllegalArgumentException("Matrix row dimensions must agree.");
	}
	if (!this.hasFullRank()) {
		throw new IllegalArgumentException("Matrix is rank deficient.");
	}
	
	// Copy right hand side
	int nx = B.columns();
	DoubleMatrix2D X = B.copy();
	
	// Compute Y = transpose(Q)*B
	for (int k = 0; k < n; k++) {
		for (int j = 0; j < nx; j++) {
			double s = 0.0; 
			for (int i = k; i < m; i++) {
				s += QR.getQuick(i,k)*X.getQuick(i,j);
			}
			s = -s / QR.getQuick(k,k);
			for (int i = k; i < m; i++) {
				X.setQuick(i,j, X.getQuick(i,j) + s*QR.getQuick(i,k));
			}
		}
	}
	// Solve R*X = Y;
	for (int k = n-1; k >= 0; k--) {
		for (int j = 0; j < nx; j++) {
			X.setQuick(k,j, X.getQuick(k,j) / Rdiag.getQuick(k));
		}
		for (int i = 0; i < k; i++) {
			for (int j = 0; j < nx; j++) {
				X.setQuick(i,j, X.getQuick(i,j) - X.getQuick(k,j)*QR.getQuick(i,k));
			}
		}
	}
	return X.viewPart(0,0,n,nx);
}
 

/**
 * Copies the columns of the indicated rows into a new sub matrix.
 * <tt>sub[0..rowIndexes.length-1,0..columnTo-columnFrom] = A[rowIndexes(:),columnFrom..columnTo]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowIndexes the indexes of the rows to copy. May be unsorted.
 * @param   columnFrom the index of the first column to copy (inclusive).
 * @param   columnTo the index of the last column to copy (inclusive).
 * @return  a new sub matrix; with <tt>sub.rows()==rowIndexes.length; sub.columns()==columnTo-columnFrom+1</tt>.
 * @throws	IndexOutOfBoundsException if <tt>columnFrom<0 || columnTo-columnFrom+1<0 || columnTo+1>matrix.columns() || for any row=rowIndexes[i]: row < 0 || row >= matrix.rows()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int[] rowIndexes, int columnFrom, int columnTo) {
	int width = columnTo-columnFrom+1;
	int rows = A.rows();
	A = A.viewPart(0,columnFrom,rows,width);
	DoubleMatrix2D sub = A.like(rowIndexes.length, width);
	
	for (int r = rowIndexes.length; --r >= 0; ) {
		int row = rowIndexes[r];
		if (row < 0 || row >= rows) 
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewRow(r).assign(A.viewRow(row));
	}
	return sub;
}
 

protected DoubleMatrix2D[] splitBlockedNN(DoubleMatrix2D A, int threshold, long flops) {
	/*
	determine how to split and parallelize best into blocks
	if more B.columns than tasks --> split B.columns, as follows:
	
			xx|xx|xxx B
			xx|xx|xxx
			xx|xx|xxx
	A
	xxx     xx|xx|xxx C 
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx

	if less B.columns than tasks --> split A.rows, as follows:
	
			xxxxxxx B
			xxxxxxx
			xxxxxxx
	A
	xxx     xxxxxxx C
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx

	*/
	//long flops = 2L*A.rows()*A.columns()*A.columns();
	int noOfTasks = (int) Math.min(flops / threshold, this.maxThreads); // each thread should process at least 30000 flops
	boolean splitHoriz = (A.columns() < noOfTasks);
	//boolean splitHoriz = (A.columns() >= noOfTasks);
	int p = splitHoriz ? A.rows() : A.columns();
	noOfTasks = Math.min(p,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		return null;
	}

	// set up concurrent tasks
	int span = p/noOfTasks;
	final DoubleMatrix2D[] blocks = new DoubleMatrix2D[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = p - span*i; // last span may be a bit larger

		final DoubleMatrix2D AA,BB,CC; 
		if (!splitHoriz) { 	// split B along columns into blocks
			blocks[i] = A.viewPart(0,offset, A.rows(), span);
		}
		else { // split A along rows into blocks
			blocks[i] = A.viewPart(offset,0,span,A.columns());
		}
	}
	return blocks;
}
 

public void dgemv(final boolean transposeA, final double alpha, DoubleMatrix2D A, final DoubleMatrix1D x, final double beta, DoubleMatrix1D y) {
	/*
	split A, as follows:
	
			x x
			x
			x
	A
	xxx     x y
	xxx     x
	---     -
	xxx     x
	xxx     x
	---     -
	xxx     x

	*/
	if (transposeA) {
		dgemv(false, alpha, A.viewDice(), x, beta, y);
		return;
	}
	int m = A.rows();
	int n = A.columns();
	long flops = 2L*m*n;
	int noOfTasks = (int) Math.min(flops / 30000, this.maxThreads); // each thread should process at least 30000 flops
	int width = A.rows();
	noOfTasks = Math.min(width,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		seqBlas.dgemv(transposeA, alpha, A, x, beta, y);
		return;
	}
	
	// set up concurrent tasks
	int span = width/noOfTasks;
	final RecursiveAction[] subTasks = new RecursiveAction[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = width - span*i; // last span may be a bit larger

		// split A along rows into blocks
		final DoubleMatrix2D AA = A.viewPart(offset,0,span,n);
		final DoubleMatrix1D yy = y.viewPart(offset,span);
				
		subTasks[i] = new RecursiveAction() { 
			@Override
			protected void compute() { 
				seqBlas.dgemv(transposeA,alpha,AA,x,beta,yy); 
				//System.out.println("Hello "+offset); 
			}
		};
	}
	
	// run tasks and wait for completion
		this.smp.taskGroup.invoke(
			new RecursiveAction() {
				@Override
				public void compute() {	
					invokeAll(subTasks);	
				}
			}
		);
}
 

public void dgemm(final boolean transposeA, final boolean transposeB, final double alpha, final DoubleMatrix2D A, final DoubleMatrix2D B, final double beta, final DoubleMatrix2D C) {
	/*
	determine how to split and parallelize best into blocks
	if more B.columns than tasks --> split B.columns, as follows:
	
			xx|xx|xxx B
			xx|xx|xxx
			xx|xx|xxx
	A
	xxx     xx|xx|xxx C 
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx

	if less B.columns than tasks --> split A.rows, as follows:
	
			xxxxxxx B
			xxxxxxx
			xxxxxxx
	A
	xxx     xxxxxxx C
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx

	*/
	if (transposeA) {
		dgemm(false, transposeB, alpha, A.viewDice(), B, beta, C);
		return;
	}
	if (transposeB) {
		dgemm(transposeA, false, alpha, A, B.viewDice(), beta, C);
		return;
	}
	int m = A.rows();
	int n = A.columns();
	int p = B.columns();

	if (B.rows() != n)
		throw new IllegalArgumentException("Matrix2D inner dimensions must agree:"+A.toStringShort()+", "+B.toStringShort());
	if (C.rows() != m || C.columns() != p)
		throw new IllegalArgumentException("Incompatibel result matrix: "+A.toStringShort()+", "+B.toStringShort()+", "+C.toStringShort());
	if (A == C || B == C)
		throw new IllegalArgumentException("Matrices must not be identical");

	long flops = 2L*m*n*p;
	int noOfTasks = (int) Math.min(flops / 30000, this.maxThreads); // each thread should process at least 30000 flops
	boolean splitB = (p >= noOfTasks);
	int width = splitB ? p : m;
	noOfTasks = Math.min(width,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		seqBlas.dgemm(transposeA, transposeB, alpha, A, B, beta, C);
		return;
	}
	
	// set up concurrent tasks
	int span = width/noOfTasks;
	final RecursiveAction[] subTasks = new RecursiveAction[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = width - span*i; // last span may be a bit larger

		final DoubleMatrix2D AA,BB,CC; 
		if (splitB) { 
			// split B along columns into blocks
			AA = A;
			BB = B.viewPart(0,offset,n,span);
			CC = C.viewPart(0,offset,m,span);
		}
		else { 
			// split A along rows into blocks
			AA = A.viewPart(offset,0,span,n);
			BB = B;
			CC = C.viewPart(offset,0,span,p);
		}
				
		subTasks[i] = new RecursiveAction() { 
			@Override
			protected void compute() { 
				seqBlas.dgemm(transposeA,transposeB,alpha,AA,BB,beta,CC); 
				//System.out.println("Hello "+offset); 
			}
		};
	}
	
	// run tasks and wait for completion
		this.smp.taskGroup.invoke(
			new RecursiveAction() {
 				@Override
				protected void compute() {	
					invokeAll(subTasks);	
				}
			}
		);
}
 

public void dgemm(final boolean transposeA, final boolean transposeB, final double alpha, final DoubleMatrix2D A, final DoubleMatrix2D B, final double beta, final DoubleMatrix2D C) {
	/*
	determine how to split and parallelize best into blocks
	if more B.columns than tasks --> split B.columns, as follows:
	
			xx|xx|xxx B
			xx|xx|xxx
			xx|xx|xxx
	A
	xxx     xx|xx|xxx C 
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx

	if less B.columns than tasks --> split A.rows, as follows:
	
			xxxxxxx B
			xxxxxxx
			xxxxxxx
	A
	xxx     xxxxxxx C
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx

	*/
	if (transposeA) {
		dgemm(false, transposeB, alpha, A.viewDice(), B, beta, C);
		return;
	}
	if (transposeB) {
		dgemm(transposeA, false, alpha, A, B.viewDice(), beta, C);
		return;
	}
	int m = A.rows();
	int n = A.columns();
	int p = B.columns();

	if (B.rows() != n)
		throw new IllegalArgumentException("Matrix2D inner dimensions must agree:"+A.toStringShort()+", "+B.toStringShort());
	if (C.rows() != m || C.columns() != p)
		throw new IllegalArgumentException("Incompatibel result matrix: "+A.toStringShort()+", "+B.toStringShort()+", "+C.toStringShort());
	if (A == C || B == C)
		throw new IllegalArgumentException("Matrices must not be identical");

	long flops = 2L*m*n*p;
	int noOfTasks = (int) Math.min(flops / 30000, this.maxThreads); // each thread should process at least 30000 flops
	boolean splitB = (p >= noOfTasks);
	int width = splitB ? p : m;
	noOfTasks = Math.min(width,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		seqBlas.dgemm(transposeA, transposeB, alpha, A, B, beta, C);
		return;
	}
	
	// set up concurrent tasks
	int span = width/noOfTasks;
	final FJTask[] subTasks = new FJTask[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = width - span*i; // last span may be a bit larger

		final DoubleMatrix2D AA,BB,CC; 
		if (splitB) { 
			// split B along columns into blocks
			AA = A;
			BB = B.viewPart(0,offset,n,span);
			CC = C.viewPart(0,offset,m,span);
		}
		else { 
			// split A along rows into blocks
			AA = A.viewPart(offset,0,span,n);
			BB = B;
			CC = C.viewPart(offset,0,span,p);
		}
				
		subTasks[i] = new FJTask() { 
			public void run() { 
				seqBlas.dgemm(transposeA,transposeB,alpha,AA,BB,beta,CC); 
				//System.out.println("Hello "+offset); 
			}
		};
	}
	
	// run tasks and wait for completion
	try { 
		this.smp.taskGroup.invoke(
			new FJTask() {
				public void run() {	
					coInvoke(subTasks);	
				}
			}
		);
	} catch (InterruptedException exc) {}
}
 

public void dgemv(final boolean transposeA, final double alpha, DoubleMatrix2D A, final DoubleMatrix1D x, final double beta, DoubleMatrix1D y) {
	/*
	split A, as follows:
	
			x x
			x
			x
	A
	xxx     x y
	xxx     x
	---     -
	xxx     x
	xxx     x
	---     -
	xxx     x

	*/
	if (transposeA) {
		dgemv(false, alpha, A.viewDice(), x, beta, y);
		return;
	}
	int m = A.rows();
	int n = A.columns();
	long flops = 2L*m*n;
	int noOfTasks = (int) Math.min(flops / 30000, this.maxThreads); // each thread should process at least 30000 flops
	int width = A.rows();
	noOfTasks = Math.min(width,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		seqBlas.dgemv(transposeA, alpha, A, x, beta, y);
		return;
	}
	
	// set up concurrent tasks
	int span = width/noOfTasks;
	final FJTask[] subTasks = new FJTask[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = width - span*i; // last span may be a bit larger

		// split A along rows into blocks
		final DoubleMatrix2D AA = A.viewPart(offset,0,span,n);
		final DoubleMatrix1D yy = y.viewPart(offset,span);
				
		subTasks[i] = new FJTask() { 
			public void run() { 
				seqBlas.dgemv(transposeA,alpha,AA,x,beta,yy); 
				//System.out.println("Hello "+offset); 
			}
		};
	}
	
	// run tasks and wait for completion
	try { 
		this.smp.taskGroup.invoke(
			new FJTask() {
				public void run() {	
					coInvoke(subTasks);	
				}
			}
		);
	} catch (InterruptedException exc) {}
}
 

protected DoubleMatrix2D[] splitBlockedNN(DoubleMatrix2D A, int threshold, long flops) {
	/*
	determine how to split and parallelize best into blocks
	if more B.columns than tasks --> split B.columns, as follows:
	
			xx|xx|xxx B
			xx|xx|xxx
			xx|xx|xxx
	A
	xxx     xx|xx|xxx C 
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx
	xxx		xx|xx|xxx

	if less B.columns than tasks --> split A.rows, as follows:
	
			xxxxxxx B
			xxxxxxx
			xxxxxxx
	A
	xxx     xxxxxxx C
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx
	xxx     xxxxxxx
	---     -------
	xxx     xxxxxxx

	*/
	//long flops = 2L*A.rows()*A.columns()*A.columns();
	int noOfTasks = (int) Math.min(flops / threshold, this.maxThreads); // each thread should process at least 30000 flops
	boolean splitHoriz = (A.columns() < noOfTasks);
	//boolean splitHoriz = (A.columns() >= noOfTasks);
	int p = splitHoriz ? A.rows() : A.columns();
	noOfTasks = Math.min(p,noOfTasks);
	
	if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
		return null;
	}

	// set up concurrent tasks
	int span = p/noOfTasks;
	final DoubleMatrix2D[] blocks = new DoubleMatrix2D[noOfTasks];
	for (int i=0; i<noOfTasks; i++) {
		final int offset = i*span;
		if (i==noOfTasks-1) span = p - span*i; // last span may be a bit larger

		final DoubleMatrix2D AA,BB,CC; 
		if (!splitHoriz) { 	// split B along columns into blocks
			blocks[i] = A.viewPart(0,offset, A.rows(), span);
		}
		else { // split A along rows into blocks
			blocks[i] = A.viewPart(offset,0,span,A.columns());
		}
	}
	return blocks;
}
 

/**
Constructs and returns a new <i>sub-range view</i> which is the sub matrix <tt>A[fromRow..toRow,fromColumn..toColumn]</tt>.
The returned matrix is backed by this matrix, so changes in the returned matrix are reflected in this matrix, and vice-versa.
Use idioms like <tt>result = subMatrix(...).copy()</tt> to generate an independent sub matrix.

@param A the source matrix.
@param fromRow   The index of the first row (inclusive).
@param toRow   The index of the last row (inclusive).
@param fromColumn   The index of the first column (inclusive).
@param toColumn   The index of the last column (inclusive).
@return a new sub-range view.
@throws	IndexOutOfBoundsException if <tt>fromColumn<0 || toColumn-fromColumn+1<0 || toColumn>=A.columns() || fromRow<0 || toRow-fromRow+1<0 || toRow>=A.rows()</tt>
*/
public DoubleMatrix2D subMatrix(DoubleMatrix2D A, int fromRow, int toRow, int fromColumn, int toColumn) {
	return A.viewPart(fromRow, fromColumn, toRow-fromRow+1, toColumn-fromColumn+1);
}
 

/**
Constructs and returns a new <i>sub-range view</i> which is the sub matrix <tt>A[fromRow..toRow,fromColumn..toColumn]</tt>.
The returned matrix is backed by this matrix, so changes in the returned matrix are reflected in this matrix, and vice-versa.
Use idioms like <tt>result = subMatrix(...).copy()</tt> to generate an independent sub matrix.

@param A the source matrix.
@param fromRow   The index of the first row (inclusive).
@param toRow   The index of the last row (inclusive).
@param fromColumn   The index of the first column (inclusive).
@param toColumn   The index of the last column (inclusive).
@return a new sub-range view.
@throws	IndexOutOfBoundsException if <tt>fromColumn<0 || toColumn-fromColumn+1<0 || toColumn>=A.columns() || fromRow<0 || toRow-fromRow+1<0 || toRow>=A.rows()</tt>
*/
public DoubleMatrix2D subMatrix(DoubleMatrix2D A, int fromRow, int toRow, int fromColumn, int toColumn) {
	return A.viewPart(fromRow, fromColumn, toRow-fromRow+1, toColumn-fromColumn+1);
}