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

/**
9 point stencil operation.
Applies a function to a moving <tt>3 x 3</tt> window.
@param A the matrix to operate on.
@param function the function to be applied to each window.
@param maxIterations the maximum number of times the stencil shall be applied to the matrix. 
	Should be a multiple of 2 because two iterations are always done in one atomic step.
@param hasConverged Convergence condition; will return before maxIterations are done when <tt>hasConverged.apply(A)==true</tt>.
	Set this parameter to <tt>null</tt> to indicate that no convergence checks shall be made.
@param convergenceIterations the number of iterations to pass between each convergence check.
	(Since a convergence may be expensive, you may want to do it only every 2,4 or 8 iterations.)
@return the number of iterations actually executed. 
*/
public static int stencil9(DoubleMatrix2D A, cern.colt.function.Double9Function function, int maxIterations, DoubleMatrix2DProcedure hasConverged, int convergenceIterations) {
	DoubleMatrix2D B = A.copy();
	if (convergenceIterations <= 1) convergenceIterations=2;
	if (convergenceIterations%2 != 0) convergenceIterations++; // odd -> make it even

	int i=0;
	while (i<maxIterations) { // do two steps at a time for efficiency
		A.zAssign8Neighbors(B,function);
		B.zAssign8Neighbors(A,function);
		i=i+2;
		if (i%convergenceIterations == 0 && hasConverged!=null) {
			if (hasConverged.apply(A)) return i;
		}
	}
	return i;
}
 

/**
 * Linear algebrax matrix-matrix multiply.
 */
protected static Double2DProcedure funMatMultLarge() {
	return new Double2DProcedure() {
		public String toString() { return "xxxxxxx";	}
		public void setParameters(DoubleMatrix2D A, DoubleMatrix2D B) {
			// do not allocate mem for "D" --> safe some mem
			this.A = A;
			this.B = B;
			this.C = A.copy();
		}
		public void init() { C.assign(0); }
		public void apply(cern.colt.Timer timer) { A.zMult(B,C); }
		public double operations() { // Mflops
			double m = A.rows();
			double n = A.columns();
			double p = B.columns();
			return 2.0*m*n*p / 1.0E6; 
		}
	};
}
 

/**
9 point stencil operation.
Applies a function to a moving <tt>3 x 3</tt> window.
@param A the matrix to operate on.
@param function the function to be applied to each window.
@param maxIterations the maximum number of times the stencil shall be applied to the matrix. 
	Should be a multiple of 2 because two iterations are always done in one atomic step.
@param hasConverged Convergence condition; will return before maxIterations are done when <tt>hasConverged.apply(A)==true</tt>.
	Set this parameter to <tt>null</tt> to indicate that no convergence checks shall be made.
@param convergenceIterations the number of iterations to pass between each convergence check.
	(Since a convergence may be expensive, you may want to do it only every 2,4 or 8 iterations.)
@return the number of iterations actually executed. 
*/
public static int stencil9(DoubleMatrix2D A, cern.colt.function.Double9Function function, int maxIterations, DoubleMatrix2DProcedure hasConverged, int convergenceIterations) {
	DoubleMatrix2D B = A.copy();
	if (convergenceIterations <= 1) convergenceIterations=2;
	if (convergenceIterations%2 != 0) convergenceIterations++; // odd -> make it even

	int i=0;
	while (i<maxIterations) { // do two steps at a time for efficiency
		A.zAssign8Neighbors(B,function);
		B.zAssign8Neighbors(A,function);
		i=i+2;
		if (i%convergenceIterations == 0 && hasConverged!=null) {
			if (hasConverged.apply(A)) return i;
		}
	}
	return i;
}
 

/** 
Solves <tt>A*X = B</tt>; returns <tt>X</tt>.
@param  B   A Matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> so that <tt>L*L'*X = B</tt>.
@exception  IllegalArgumentException  if <tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if <tt>!isSymmetricPositiveDefinite()</tt>.
*/
private DoubleMatrix2D XXXsolveBuggy(DoubleMatrix2D B) {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	if (B.rows() != n) {
		throw new IllegalArgumentException("Matrix row dimensions must agree.");
	}
	if (!isSymmetricPositiveDefinite) {
		throw new IllegalArgumentException("Matrix is not symmetric positive definite.");
	}
	
	// Copy right hand side.
	DoubleMatrix2D X = B.copy();
	int nx = B.columns();
	
	// precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] Xrows = new DoubleMatrix1D[n];
	for (int k = 0; k < n; k++) Xrows[k] = X.viewRow(k);
	
	// Solve L*Y = B;
	for (int k = 0; k < n; k++) {
		for (int i = k+1; i < n; i++) {
			// X[i,j] -= X[k,j]*L[i,k]
			Xrows[i].assign(Xrows[k], F.minusMult(L.getQuick(i,k)));
		}
		Xrows[k].assign(F.div(L.getQuick(k,k)));
	}
	
	// Solve L'*X = Y;
	for (int k = n-1; k >= 0; k--) {
		Xrows[k].assign(F.div(L.getQuick(k,k)));
		for (int i = 0; i < k; i++) {
			// X[i,j] -= X[k,j]*L[k,i]
			Xrows[i].assign(Xrows[k], F.minusMult(L.getQuick(k,i)));
		}
	}
	return X;
}
 

/**
 * Demonstrates sorting with precomputation of aggregates, comparing mergesort with quicksort.
 */
public static void zdemo7(int rows, int columns, boolean print) {
	// for reliable benchmarks, call this method twice: once with small dummy parameters to "warm up" the jitter, then with your real work-load
		
	System.out.println("\n\n");
	System.out.println("now initializing... ");
		
	final cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D A = cern.colt.matrix.DoubleFactory2D.dense.make(rows,columns);
	A.assign(new cern.jet.random.engine.DRand()); // initialize randomly
	DoubleMatrix2D B = A.copy();

	double[] v1 = A.viewColumn(0).toArray();
	double[] v2 = A.viewColumn(0).toArray();
	System.out.print("now quick sorting... ");
	cern.colt.Timer timer = new cern.colt.Timer().start();
	quickSort.sort(A,0);
	timer.stop().display();

	System.out.print("now merge sorting... ");
	timer.reset().start();
	mergeSort.sort(A,0);
	timer.stop().display();

	System.out.print("now quick sorting with simple aggregation... ");
	timer.reset().start();
	quickSort.sort(A,v1);
	timer.stop().display();

	System.out.print("now merge sorting with simple aggregation... ");
	timer.reset().start();
	mergeSort.sort(A,v2);
	timer.stop().display();
}
 

/** 
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);
}
 

/**
 * Linear algebraic matrix power; <tt>B = A<sup>k</sup> <==> B = A*A*...*A</tt>.
 * @param A the source matrix; must be square.
 * @param k the exponent, can be any number.
 * @return a new result matrix.
 * 
 * @throws IllegalArgumentException if <tt>!Testing.isSquare(A)</tt>.
 */
private DoubleMatrix2D xpowSlow(DoubleMatrix2D A, int k) {
	//cern.colt.Timer timer = new cern.colt.Timer().start();
	DoubleMatrix2D result = A.copy();
	for (int i=0; i<k-1; i++) {
		result = mult(result,A);
	}
	//timer.stop().display();
	return result;
}
 

/**
 * Returns the inverse or pseudo-inverse of matrix <tt>A</tt>.
 * @return a new independent matrix; inverse(matrix) if the matrix is square, pseudoinverse otherwise.
 */
public DoubleMatrix2D inverse(DoubleMatrix2D A) {
	if (property.isSquare(A) && property.isDiagonal(A)) {
		DoubleMatrix2D inv = A.copy();
		boolean isNonSingular = Diagonal.inverse(inv);
		if (!isNonSingular) throw new IllegalArgumentException("A is singular.");
		return inv;
	}
	return solve(A, DoubleFactory2D.dense.identity(A.rows()));
}
 

/**
 * Sets the matrices to operate upon.
 */
public void setParameters(DoubleMatrix2D A, DoubleMatrix2D B) {
	this.A=A;
	this.B=B;
	this.C=A.copy();
	this.D=B.copy();
}
 

/**
 * Demonstrates sorting with precomputation of aggregates, comparing mergesort with quicksort.
 */
public static void zdemo7(int rows, int columns, boolean print) {
	// for reliable benchmarks, call this method twice: once with small dummy parameters to "warm up" the jitter, then with your real work-load
		
	System.out.println("\n\n");
	System.out.println("now initializing... ");
		
	final cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D A = cern.colt.matrix.DoubleFactory2D.dense.make(rows,columns);
	A.assign(new cern.jet.random.engine.DRand()); // initialize randomly
	DoubleMatrix2D B = A.copy();

	double[] v1 = A.viewColumn(0).toArray();
	double[] v2 = A.viewColumn(0).toArray();
	System.out.print("now quick sorting... ");
	cern.colt.Timer timer = new cern.colt.Timer().start();
	quickSort.sort(A,0);
	timer.stop().display();

	System.out.print("now merge sorting... ");
	timer.reset().start();
	mergeSort.sort(A,0);
	timer.stop().display();

	System.out.print("now quick sorting with simple aggregation... ");
	timer.reset().start();
	quickSort.sort(A,v1);
	timer.stop().display();

	System.out.print("now merge sorting with simple aggregation... ");
	timer.reset().start();
	mergeSort.sort(A,v2);
	timer.stop().display();
}
 

/** 
Solves <tt>A*X = B</tt>; returns <tt>X</tt>.
@param  B   A Matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> so that <tt>L*L'*X = B</tt>.
@exception  IllegalArgumentException  if <tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if <tt>!isSymmetricPositiveDefinite()</tt>.
*/
private DoubleMatrix2D XXXsolveBuggy(DoubleMatrix2D B) {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	if (B.rows() != n) {
		throw new IllegalArgumentException("Matrix row dimensions must agree.");
	}
	if (!isSymmetricPositiveDefinite) {
		throw new IllegalArgumentException("Matrix is not symmetric positive definite.");
	}
	
	// Copy right hand side.
	DoubleMatrix2D X = B.copy();
	int nx = B.columns();
	
	// precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] Xrows = new DoubleMatrix1D[n];
	for (int k = 0; k < n; k++) Xrows[k] = X.viewRow(k);
	
	// Solve L*Y = B;
	for (int k = 0; k < n; k++) {
		for (int i = k+1; i < n; i++) {
			// X[i,j] -= X[k,j]*L[i,k]
			Xrows[i].assign(Xrows[k], F.minusMult(L.getQuick(i,k)));
		}
		Xrows[k].assign(F.div(L.getQuick(k,k)));
	}
	
	// Solve L'*X = Y;
	for (int k = n-1; k >= 0; k--) {
		Xrows[k].assign(F.div(L.getQuick(k,k)));
		for (int i = 0; i < k; i++) {
			// X[i,j] -= X[k,j]*L[k,i]
			Xrows[i].assign(Xrows[k], F.minusMult(L.getQuick(k,i)));
		}
	}
	return X;
}
 

/** 
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);
}
 

/**
 * Linear algebraic matrix power; <tt>B = A<sup>k</sup> <==> B = A*A*...*A</tt>.
 * @param A the source matrix; must be square.
 * @param k the exponent, can be any number.
 * @return a new result matrix.
 * 
 * @throws IllegalArgumentException if <tt>!Testing.isSquare(A)</tt>.
 */
private DoubleMatrix2D xpowSlow(DoubleMatrix2D A, int k) {
	//cern.colt.Timer timer = new cern.colt.Timer().start();
	DoubleMatrix2D result = A.copy();
	for (int i=0; i<k-1; i++) {
		result = mult(result,A);
	}
	//timer.stop().display();
	return result;
}
 

/**
 * Returns the inverse or pseudo-inverse of matrix <tt>A</tt>.
 * @return a new independent matrix; inverse(matrix) if the matrix is square, pseudoinverse otherwise.
 */
public DoubleMatrix2D inverse(DoubleMatrix2D A) {
	if (property.isSquare(A) && property.isDiagonal(A)) {
		DoubleMatrix2D inv = A.copy();
		boolean isNonSingular = Diagonal.inverse(inv);
		if (!isNonSingular) throw new IllegalArgumentException("A is singular.");
		return inv;
	}
	return solve(A, DoubleFactory2D.dense.identity(A.rows()));
}
 

/**
 * Sets the matrices to operate upon.
 */
public void setParameters(DoubleMatrix2D A, DoubleMatrix2D B) {
	this.A=A;
	this.B=B;
	this.C=A.copy();
	this.D=B.copy();
}
 

/** 
Constructs and returns a new QR decomposition object;  computed by Householder reflections;
The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
@param A    A rectangular matrix.
@return     a decomposition object to access <tt>R</tt> and the Householder vectors <tt>H</tt>, and to compute <tt>Q</tt>.
@throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
*/

public QRDecomposition (DoubleMatrix2D A) {
	Property.DEFAULT.checkRectangular(A);

	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	// Initialize.
	QR = A.copy();
	m = A.rows();
	n = A.columns();
	Rdiag = A.like1D(n);
	//Rdiag = new double[n];
	cern.colt.function.DoubleDoubleFunction hypot = Algebra.hypotFunction();
	
	// precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] QRcolumns = new DoubleMatrix1D[n];
	DoubleMatrix1D[] QRcolumnsPart = new DoubleMatrix1D[n];
	for (int k = 0; k < n; k++) {
		QRcolumns[k] = QR.viewColumn(k);
		QRcolumnsPart[k] = QR.viewColumn(k).viewPart(k,m-k);
	}
	
	// Main loop.
	for (int k = 0; k < n; k++) {
		//DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		// Compute 2-norm of k-th column without under/overflow.
		double nrm = 0;
		//if (k<m) nrm = QRcolumnsPart[k].aggregate(hypot,F.identity);
		
		for (int i = k; i < m; i++) { // fixes bug reported by [email protected]
			nrm = Algebra.hypot(nrm,QR.getQuick(i,k));
		}
		
		
		if (nrm != 0.0) {
			// Form k-th Householder vector.
			if (QR.getQuick(k,k) < 0) nrm = -nrm;
			QRcolumnsPart[k].assign(cern.jet.math.Functions.div(nrm));		
			/*
			for (int i = k; i < m; i++) {
			   QR[i][k] /= nrm;
			}
			*/
			
			QR.setQuick(k,k, QR.getQuick(k,k) + 1);
			
			// Apply transformation to remaining columns.
			for (int j = k+1; j < n; j++) {
				DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
				double s = QRcolumnsPart[k].zDotProduct(QRcolj);
				/*
				// fixes bug reported by John Chambers
				DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
				double s = QRcolumnsPart[k].zDotProduct(QRcolumns[j]);
				double s = 0.0; 
				for (int i = k; i < m; i++) {
				  s += QR[i][k]*QR[i][j];
				}
				*/
				s = -s / QR.getQuick(k,k);
				//QRcolumnsPart[j].assign(QRcolumns[k], F.plusMult(s));
				
				for (int i = k; i < m; i++) {
				  QR.setQuick(i,j, QR.getQuick(i,j) + s*QR.getQuick(i,k));
				}
				
			}
		}
		Rdiag.setQuick(k, -nrm);
	}
}
 

/**
 * Linear algebraic matrix power; <tt>B = A<sup>k</sup> <==> B = A*A*...*A</tt>.
 * <ul>
 * <li><tt>p &gt;= 1: B = A*A*...*A</tt>.</li>
 * <li><tt>p == 0: B = identity matrix</tt>.</li>
 * <li><tt>p &lt;  0: B = pow(inverse(A),-p)</tt>.</li>
 * </ul>
 * Implementation: Based on logarithms of 2, memory usage minimized.
 * @param A the source matrix; must be square; stays unaffected by this operation.
 * @param p the exponent, can be any number.
 * @return <tt>B</tt>, a newly constructed result matrix; storage-independent of <tt>A</tt>.
 * 
 * @throws IllegalArgumentException if <tt>!property().isSquare(A)</tt>.
 */
public DoubleMatrix2D pow(DoubleMatrix2D A, int p) {
	// matrix multiplication based on log2 method: A*A*....*A is slow, ((A * A)^2)^2 * ... is faster
	// allocates two auxiliary matrices as work space
/*
	Blas blas = SmpBlas.smpBlas; // for parallel matrix mult; if not initialized defaults to sequential blas
	Property.DEFAULT.checkSquare(A);
	if (p<0) {
		A = inverse(A);
		p = -p;
	}
	if (p==0) return DoubleFactory2D.dense.identity(A.rows());
	DoubleMatrix2D T = A.like(); // temporary
	if (p==1) return T.assign(A);  // safes one auxiliary matrix allocation
	if (p==2) {
		blas.dgemm(false,false,1,A,A,0,T); // mult(A,A); // safes one auxiliary matrix allocation
		return T;
	}

	int k = cern.colt.bitvector.QuickBitVector.mostSignificantBit(p); // index of highest bit in state "true"
*/	
	/*
	this is the naive version:*/
	DoubleMatrix2D B = A.copy();
	for (int i=0; i<p-1; i++) {
		B = mult(B,A);
	}
	return B;
	

	// here comes the optimized version:
	//cern.colt.Timer timer = new cern.colt.Timer().start();

/*	int i=0;
	while (i<=k && (p & (1<<i)) == 0) { // while (bit i of p == false)
		// A = mult(A,A); would allocate a lot of temporary memory
		blas.dgemm(false,false,1,A,A,0,T); // A.zMult(A,T);
		DoubleMatrix2D swap = A; A = T; T = swap; // swap A with T
		i++;
	}

	DoubleMatrix2D B = A.copy();
	i++;
	for (; i<=k; i++) {
		// A = mult(A,A); would allocate a lot of temporary memory
		blas.dgemm(false,false,1,A,A,0,T); // A.zMult(A,T);	
		DoubleMatrix2D swap = A; A = T; T = swap; // swap A with T

		if ((p & (1<<i)) != 0) { // if (bit i of p == true)
			// B = mult(B,A); would allocate a lot of temporary memory
			blas.dgemm(false,false,1,B,A,0,T); // B.zMult(A,T);		
			swap = B; B = T; T = swap; // swap B with T
		}
	}
	//timer.stop().display();
	return B;*/
}
 

/** 
Constructs and returns a new QR decomposition object;  computed by Householder reflections;
The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
@param A    A rectangular matrix.
@return     a decomposition object to access <tt>R</tt> and the Householder vectors <tt>H</tt>, and to compute <tt>Q</tt>.
@throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
*/

public QRDecomposition (DoubleMatrix2D A) {
	Property.DEFAULT.checkRectangular(A);

	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	// Initialize.
	QR = A.copy();
	m = A.rows();
	n = A.columns();
	Rdiag = A.like1D(n);
	//Rdiag = new double[n];
	cern.colt.function.DoubleDoubleFunction hypot = Algebra.hypotFunction();
	
	// precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] QRcolumns = new DoubleMatrix1D[n];
	DoubleMatrix1D[] QRcolumnsPart = new DoubleMatrix1D[n];
	for (int k = 0; k < n; k++) {
		QRcolumns[k] = QR.viewColumn(k);
		QRcolumnsPart[k] = QR.viewColumn(k).viewPart(k,m-k);
	}
	
	// Main loop.
	for (int k = 0; k < n; k++) {
		//DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		// Compute 2-norm of k-th column without under/overflow.
		double nrm = 0;
		//if (k<m) nrm = QRcolumnsPart[k].aggregate(hypot,F.identity);
		
		for (int i = k; i < m; i++) { // fixes bug reported by [email protected]
			nrm = Algebra.hypot(nrm,QR.getQuick(i,k));
		}
		
		
		if (nrm != 0.0) {
			// Form k-th Householder vector.
			if (QR.getQuick(k,k) < 0) nrm = -nrm;
			QRcolumnsPart[k].assign(cern.jet.math.Functions.div(nrm));		
			/*
			for (int i = k; i < m; i++) {
			   QR[i][k] /= nrm;
			}
			*/
			
			QR.setQuick(k,k, QR.getQuick(k,k) + 1);
			
			// Apply transformation to remaining columns.
			for (int j = k+1; j < n; j++) {
				DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
				double s = QRcolumnsPart[k].zDotProduct(QRcolj);
				/*
				// fixes bug reported by John Chambers
				DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
				double s = QRcolumnsPart[k].zDotProduct(QRcolumns[j]);
				double s = 0.0; 
				for (int i = k; i < m; i++) {
				  s += QR[i][k]*QR[i][j];
				}
				*/
				s = -s / QR.getQuick(k,k);
				//QRcolumnsPart[j].assign(QRcolumns[k], F.plusMult(s));
				
				for (int i = k; i < m; i++) {
				  QR.setQuick(i,j, QR.getQuick(i,j) + s*QR.getQuick(i,k));
				}
				
			}
		}
		Rdiag.setQuick(k, -nrm);
	}
}
 

/** 
Solves <tt>A*X = B</tt>.
@param  B   A matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> so that <tt>L*U*X = B(piv,:)</tt>.
@exception  IllegalArgumentException if </tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if A is singular, that is, if <tt>!this.isNonsingular()</tt>.
@exception  IllegalArgumentException  if <tt>A.rows() < A.columns()</tt>.
*/

public DoubleMatrix2D solve(DoubleMatrix2D B) {
	DoubleMatrix2D X = B.copy();
	quick.solve(X);
	return X;
}
 

/** 
Solves <tt>A*X = B</tt>.
@param  B   A matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> so that <tt>L*U*X = B(piv,:)</tt>.
@exception  IllegalArgumentException if </tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if A is singular, that is, if <tt>!this.isNonsingular()</tt>.
@exception  IllegalArgumentException  if <tt>A.rows() < A.columns()</tt>.
*/

public DoubleMatrix2D solve(DoubleMatrix2D B) {
	DoubleMatrix2D X = B.copy();
	quick.solve(X);
	return X;
}