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

private static DoubleMatrix2D blockDiagonalExponential(double distance, DoubleMatrix2D mat) {
    for (int i = 0; i < mat.rows(); i++) {
        if ((i + 1) < mat.rows() && mat.getQuick(i, i + 1) != 0) {
            double a = mat.getQuick(i, i);
            double b = mat.getQuick(i, i + 1);
            double expat = Math.exp(distance * a);
            double cosbt = Math.cos(distance * b);
            double sinbt = Math.sin(distance * b);
            mat.setQuick(i, i, expat * cosbt);
            mat.setQuick(i + 1, i + 1, expat * cosbt);
            mat.setQuick(i, i + 1, expat * sinbt);
            mat.setQuick(i + 1, i, -expat * sinbt);
            i++; // processed two entries in loop
        } else
            mat.setQuick(i, i, Math.exp(distance * mat.getQuick(i, i))); // 1x1 block
    }
    return mat;
}
 

public static DoubleMatrix2D getOrthogVectors(DoubleMatrix2D M){
    //Get a matrix whose columns are orthogonal to the row space of M
    //which expected to be nonsingular
    DoubleMatrix2D Maug = DoubleFactory2D.dense.make(M.columns(), M.columns());
    Maug.viewPart(0, 0, M.rows(), M.columns()).assign(M);

    SingularValueDecomposition svd = new SingularValueDecomposition(Maug);

    int numOrthVecs = M.columns() - M.rows();
    if(svd.rank() != M.rows()){
        throw new RuntimeException("ERROR: Singularity in constr jac.  Rank: "+svd.rank());
    }

    DoubleMatrix2D orthVecs = svd.getV().viewPart(0, M.rows(), M.columns(), numOrthVecs);

    //DEBUG!!!  Should be 0 and identity respecitvely
    /*DoubleMatrix2D check = Algebra.DEFAULT.mult(M, orthVecs);
    DoubleMatrix2D checkOrth = Algebra.DEFAULT.mult( Algebra.DEFAULT.transpose(orthVecs), orthVecs);
    */


    return orthVecs;
}
 

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

/**
 * A matrix <tt>A</tt> is an <i>identity</i> matrix if <tt>A[i,i] == 1</tt> and all other cells are zero.
 * Matrix may but need not be square.
 */
public boolean isIdentity(DoubleMatrix2D A) {
	double epsilon = tolerance();
	int rows = A.rows();
	int columns = A.columns();
	for (int row = rows; --row >=0; ) {
		for (int column = columns; --column >= 0; ) {
			double v = A.getQuick(row,column);
			if (row==column) {
				if (!(Math.abs(1-v) < epsilon)) return false;
			}
			else if (!(Math.abs(v) <= epsilon)) return false;
		}
	}
	return true;
}
 

/**
Modifies the matrix to be a lower trapezoidal matrix.
@return <tt>A</tt> (for convenience only).
@see #triangulateLower(DoubleMatrix2D)
*/
protected DoubleMatrix2D trapezoidalLower(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	for (int r = rows; --r >= 0; ) {
		for (int c = columns; --c >= 0; ) {
			if (r < c) A.setQuick(r,c, 0);
		}
	}
	return A;
}
 

/**
 * A matrix <tt>A</tt> is <i>strictly upper triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i &gt;= j</tt>.
 * Matrix may but need not be square.
 */
public boolean isStrictlyUpperTriangular(DoubleMatrix2D A) {
	double epsilon = tolerance();
	int rows = A.rows();
	int columns = A.columns();
	for (int column = columns; --column >= 0; ) {
		for (int row = rows; --row >= column; ) {
			//if (A.getQuick(row,column) != 0) return false;
			if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
		}
	}
	return true;
}
 

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

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

/**
 * Returns the infinity norm of matrix <tt>A</tt>, which is the maximum absolute row sum.
 */
public double normInfinity(DoubleMatrix2D A) {
	double max = 0;
	for (int row = A.rows(); --row >=0; ) {
		//max = Math.max(max, normInfinity(A.viewRow(row)));
		max = Math.max(max, norm1(A.viewRow(row)));
	}
	return max;
}
 

/**
 * A matrix <tt>A</tt> is <i>strictly lower triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i &lt;= j</tt>.
 * Matrix may but need not be square.
 */
public boolean isStrictlyLowerTriangular(DoubleMatrix2D A) {
	double epsilon = tolerance();
	int rows = A.rows();
	int columns = A.columns();
	for (int column = columns; --column >= 0; ) {
		for (int row = Math.min(rows,column+1); --row >= 0; ) {
			//if (A.getQuick(row,column) != 0) return false;
			if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
		}
	}
	return true;
}
 

/**
 * A matrix <tt>A</tt> is <i>upper bidiagonal</i> if <tt>A[i,j]==0</tt> unless <tt>i==j || i==j-1</tt>.
 * Matrix may but need not be square.
 */
public boolean isUpperBidiagonal(DoubleMatrix2D A) {
	double epsilon = tolerance();
	int rows = A.rows();
	int columns = A.columns();
	for (int row = rows; --row >=0; ) {
		for (int column = columns; --column >= 0; ) {
			if (!(row==column || row==column-1)) {
				//if (A.getQuick(row,column) != 0) return false;
				if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
			}
		}
	}
	return true;
}
 

/**
Modifies the matrix to be an upper triangular matrix.
@return <tt>A</tt> (for convenience only).
@see #triangulateLower(DoubleMatrix2D)
*/
protected DoubleMatrix2D upperTriangular(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	int min = Math.min(rows,columns);
	for (int r = min; --r >= 0; ) {
		for (int c = min; --c >= 0; ) {
			if (r > c) A.setQuick(r,c, 0);
		}
	}
	if (columns<rows) A.viewPart(min,0,rows-min,columns).assign(0);

	return A;
}
 

/**
 * Returns the infinity norm of matrix <tt>A</tt>, which is the maximum absolute row sum.
 */
public double normInfinity(DoubleMatrix2D A) {
	double max = 0;
	for (int row = A.rows(); --row >=0; ) {
		//max = Math.max(max, normInfinity(A.viewRow(row)));
		max = Math.max(max, norm1(A.viewRow(row)));
	}
	return max;
}
 

/**
Modifies the given matrix <tt>A</tt> such that it's rows are permuted as specified; Useful for pivoting.
Row <tt>A[i]</tt> will go into row <tt>A[indexes[i]]</tt>.
<p>
<b>Example:</b>
<pre>
Reordering
[A,B,C,D,E] with indexes [0,4,2,3,1] yields 
[A,E,C,D,B]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[2], A[3]<--A[3], A[4]<--A[1].

Reordering
[A,B,C,D,E] with indexes [0,4,1,2,3] yields 
[A,E,B,C,D]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[1], A[3]<--A[2], A[4]<--A[3].
</pre>

@param   A   the matrix to permute.
@param   indexes the permutation indexes, must satisfy <tt>indexes.length==A.rows() && indexes[i] >= 0 && indexes[i] < A.rows()</tt>;
@param   work the working storage, must satisfy <tt>work.length >= A.rows()</tt>; set <tt>work==null</tt> if you don't care about performance.
@return the modified <tt>A</tt> (for convenience only).
@throws	IndexOutOfBoundsException if <tt>indexes.length != A.rows()</tt>.
*/
public DoubleMatrix2D permuteRows(final DoubleMatrix2D A, int[] indexes, int[] work) {
	// check validity
	int size = A.rows();
	if (indexes.length != size) throw new IndexOutOfBoundsException("invalid permutation");

	/*
	int i=size;
	int a;
	while (--i >= 0 && (a=indexes[i])==i) if (a < 0 || a >= size) throw new IndexOutOfBoundsException("invalid permutation");
	if (i<0) return; // nothing to permute
	*/

	int columns = A.columns();
	if (columns < size/10) { // quicker
		double[] doubleWork = new double[size];
		for (int j=A.columns(); --j >= 0; ) permute(A.viewColumn(j), indexes, doubleWork);
		return A;
	}

	cern.colt.Swapper swapper = new cern.colt.Swapper() {
		public void swap(int a, int b) {
			A.viewRow(a).swap(A.viewRow(b));
		}
	};

	cern.colt.GenericPermuting.permute(indexes, swapper, work, null);
	return A;
}
 

/**
Decomposes matrix <tt>A</tt> into <tt>L</tt> and <tt>U</tt> (in-place).
Upon return <tt>A</tt> is overridden with the result <tt>LU</tt>, such that <tt>L*U = A</tt>.
Uses a "left-looking", dot-product, Crout/Doolittle algorithm.
@param  A   any matrix.
*/	
public void decompose(DoubleMatrix2D A) {
	final int CUT_OFF = 10;
	// setup
	LU = A;
	int m = A.rows();
	int n = A.columns();

	// setup pivot vector
	if (this.piv==null || this.piv.length != m) this.piv = new int[m];
	for (int i = m; --i >= 0; ) piv[i] = i;
	pivsign = 1;

	if (m*n == 0) {
		setLU(LU);
		return; // nothing to do
	}
	
	//precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] LUrows = new DoubleMatrix1D[m];
	for (int i = 0; i < m; i++) LUrows[i] = LU.viewRow(i);
	
	cern.colt.list.IntArrayList nonZeroIndexes = new cern.colt.list.IntArrayList(); // sparsity
	DoubleMatrix1D LUcolj = LU.viewColumn(0).like();  // blocked column j
	cern.jet.math.Mult multFunction = cern.jet.math.Mult.mult(0); 

	// Outer loop.
	for (int j = 0; j < n; j++) {
		// blocking (make copy of j-th column to localize references)
		LUcolj.assign(LU.viewColumn(j));
		
		// sparsity detection
		int maxCardinality = m/CUT_OFF; // == heuristic depending on speedup
		LUcolj.getNonZeros(nonZeroIndexes,null,maxCardinality);
		int cardinality = nonZeroIndexes.size(); 
		boolean sparse = (cardinality < maxCardinality);

		// Apply previous transformations.
		for (int i = 0; i < m; i++) {
			int kmax = Math.min(i,j);
			double s;
			if (sparse) {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax,nonZeroIndexes);
			}
			else {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax);
			}
			double before = LUcolj.getQuick(i);
			double after = before -s;
			LUcolj.setQuick(i, after); // LUcolj is a copy
			LU.setQuick(i,j, after);   // this is the original
			if (sparse) {
				if (before==0 && after!=0) { // nasty bug fixed!
					int pos = nonZeroIndexes.binarySearch(i);
					pos = -pos -1;
					nonZeroIndexes.beforeInsert(pos,i);
				}
				if (before!=0 && after==0) {
					nonZeroIndexes.remove(nonZeroIndexes.binarySearch(i));
				}
			}
		}
	
		// Find pivot and exchange if necessary.
		int p = j;
		if (p < m) {
			double max = Math.abs(LUcolj.getQuick(p));
			for (int i = j+1; i < m; i++) {
				double v = Math.abs(LUcolj.getQuick(i));
				if (v > max) {
					p = i;
					max = v;
				}
			}
		}
		if (p != j) {
			LUrows[p].swap(LUrows[j]);
			int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
			pivsign = -pivsign;
		}
		
		// Compute multipliers.
		double jj;
		if (j < m && (jj=LU.getQuick(j,j)) != 0.0) {
			multFunction.multiplicator = 1 / jj;
			LU.viewColumn(j).viewPart(j+1,m-(j+1)).assign(multFunction);
		}
		
	}
	setLU(LU);
}
 

/**
 * Checks whether the given matrix <tt>A</tt> is <i>rectangular</i>.
 * @throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
 */
public void checkRectangular(DoubleMatrix2D A) {
	if (A.rows() < A.columns()) {
		throw new IllegalArgumentException("Matrix must be rectangular: "+cern.colt.matrix.doublealgo.Formatter.shape(A));
	}
}
 

/**
Same as {@link cern.colt.Partitioning#partition(int[],int,int,int[],int,int,int[])}
except that it <i>synchronously</i> partitions the rows of the given matrix by the values of the given matrix column;
This is essentially the same as partitioning a list of composite objects by some instance variable;
In other words, two entire rows of the matrix are swapped, whenever two column values indicate so.
<p>
Let's say, a "row" is an "object" (tuple, d-dimensional point).
A "column" is the list of "object" values of a given variable (field, dimension).
A "matrix" is a list of "objects" (tuples, points).
<p>
Now, rows (objects, tuples) are partially sorted according to their values in one given variable (dimension).
Two entire rows of the matrix are swapped, whenever two column values indicate so.
<p>
Note that arguments are not checked for validity.
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt>8 x 3 matrix:<br>
	  23, 22, 21<br>
	  20, 19, 18<br>
	  17, 16, 15<br>
	  14, 13, 12<br>
	  11, 10, 9<br>
	  8,  7,  6<br>
	  5,  4,  3<br>
	  2,  1,  0 </tt></td>
	<td align="left" valign="top"> 
	    <tt>column = 0;<br>
		splitters = {5,10,12}<br>
		partition(matrix,column,splitters,splitIndexes);<br>
		==><br>
		splitIndexes == {0, 2, 3}</tt></p>
	  </td>
	<td valign="top">
	  The matrix IS NOT REORDERED.<br>
	  The new VIEW IS REORDERED:<br>
	  <tt>8 x 3 matrix:<br>
	  2,  1,  0<br>
	  5,  4,  3<br>
	  8,  7,  6<br>
	  11, 10, 9<br>
	  23, 22, 21<br>
	  20, 19, 18<br>
	  17, 16, 15<br>
	  14, 13, 12 </tt></td>
  </tr>
</table>
@param matrix the matrix to be partitioned.
@param column the index of the column to partition on.
@param splitters the values at which the rows shall be split into intervals.
	Must be sorted ascending and must not contain multiple identical values.
	These preconditions are not checked; be sure that they are met.
 
@param splitIndexes a list into which this method fills the indexes of rows delimiting intervals.
Therefore, must satisfy <tt>splitIndexes.length >= splitters.length</tt>.

@return a new matrix view having rows partitioned by the given column and splitters.
*/
public static DoubleMatrix2D partition(DoubleMatrix2D matrix, int column, final double[] splitters, int[] splitIndexes) {
	int rowFrom = 0;
	int rowTo = matrix.rows()-1;
	int splitFrom = 0;
	int splitTo = splitters.length-1;
	int[] rowIndexes = new int[matrix.rows()]; // row indexes to reorder instead of matrix itself
	for (int i=rowIndexes.length; --i >= 0; ) rowIndexes[i] = i;

	partition(matrix,rowIndexes,rowFrom,rowTo,column,splitters,splitFrom,splitTo,splitIndexes);

	// take all columns in the original order
	int[] columnIndexes = new int[matrix.columns()];
	for (int i=columnIndexes.length; --i >= 0; ) columnIndexes[i] = i;

	// view the matrix according to the reordered row indexes
	return matrix.viewSelection(rowIndexes,columnIndexes);
}
 

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

/**
Modifies the matrix to be a lower triangular matrix.
<p>
<b>Examples:</b> 
<table border="0">
  <tr nowrap> 
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9 </td>
	<td align="center">triang.Upper<br>
	  ==></td>
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  0, 9, 9, 9, 9<br>
	  0, 0, 9, 9, 9</td>
  </tr>
  <tr nowrap> 
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9 </td>
	<td align="center">triang.Upper<br>
	  ==></td>
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  0, 9, 9<br>
	  0, 0, 9<br>
	  0, 0, 0<br>
	  0, 0, 0</td>
  </tr>
  <tr nowrap> 
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9 </td>
	<td align="center">triang.Lower<br>
	  ==></td>
	<td valign="top">3 x 5 matrix:<br>
	  1, 0, 0, 0, 0<br>
	  9, 1, 0, 0, 0<br>
	  9, 9, 1, 0, 0</td>
  </tr>
  <tr nowrap> 
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9 </td>
	<td align="center">triang.Lower<br>
	  ==></td>
	<td valign="top">5 x 3 matrix:<br>
	  1, 0, 0<br>
	  9, 1, 0<br>
	  9, 9, 1<br>
	  9, 9, 9<br>
	  9, 9, 9</td>
  </tr>
</table>

@return <tt>A</tt> (for convenience only).
@see #triangulateUpper(DoubleMatrix2D)
*/
protected DoubleMatrix2D lowerTriangular(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	int min = Math.min(rows,columns);
	for (int r = min; --r >= 0; ) {
		for (int c = min; --c >= 0; ) {
			if (r < c) A.setQuick(r,c, 0);
			else if (r == c) A.setQuick(r,c, 1);
		}
	}
	if (columns>rows) A.viewPart(0,min,rows,columns-min).assign(0);

	return A;
}
 

/**
Constructs and returns a sampling view with <tt>round(matrix.rows() * rowFraction)</tt> rows and <tt>round(matrix.columns() * columnFraction)</tt> columns.
Samples "without replacement".
Rows and columns are randomly chosen from the uniform distribution.
Examples: 
<table border="1" cellspacing="0">
  <tr valign="top" align="center"> 
	<td> 
	  <div align="left"><tt>matrix</tt></div>
	</td>
	<td> 
	  <div align="left"><tt>rowFraction=0.2<br>
		columnFraction=0.2</tt></div>
	</td>
	<td> 
	  <div align="left"><tt>rowFraction=0.2<br>
		columnFraction=1.0 </tt></div>
	</td>
	<td> 
	  <div align="left"><tt>rowFraction=1.0<br>
		columnFraction=0.2 </tt></div>
	</td>
  </tr>
  <tr valign="top"> 
	<td><tt> 10&nbsp;x&nbsp;10&nbsp;matrix<br>
	  &nbsp;1&nbsp;&nbsp;2&nbsp;&nbsp;3&nbsp;&nbsp;4&nbsp;&nbsp;5&nbsp;&nbsp;6&nbsp;&nbsp;7&nbsp;&nbsp;8&nbsp;&nbsp;9&nbsp;&nbsp;10<br>
	  11&nbsp;12&nbsp;13&nbsp;14&nbsp;15&nbsp;16&nbsp;17&nbsp;18&nbsp;19&nbsp;&nbsp;20<br>
	  21&nbsp;22&nbsp;23&nbsp;24&nbsp;25&nbsp;26&nbsp;27&nbsp;28&nbsp;29&nbsp;&nbsp;30<br>
	  31&nbsp;32&nbsp;33&nbsp;34&nbsp;35&nbsp;36&nbsp;37&nbsp;38&nbsp;39&nbsp;&nbsp;40<br>
	  41&nbsp;42&nbsp;43&nbsp;44&nbsp;45&nbsp;46&nbsp;47&nbsp;48&nbsp;49&nbsp;&nbsp;50<br>
	  51&nbsp;52&nbsp;53&nbsp;54&nbsp;55&nbsp;56&nbsp;57&nbsp;58&nbsp;59&nbsp;&nbsp;60<br>
	  61&nbsp;62&nbsp;63&nbsp;64&nbsp;65&nbsp;66&nbsp;67&nbsp;68&nbsp;69&nbsp;&nbsp;70<br>
	  71&nbsp;72&nbsp;73&nbsp;74&nbsp;75&nbsp;76&nbsp;77&nbsp;78&nbsp;79&nbsp;&nbsp;80<br>
	  81&nbsp;82&nbsp;83&nbsp;84&nbsp;85&nbsp;86&nbsp;87&nbsp;88&nbsp;89&nbsp;&nbsp;90<br>
	  91&nbsp;92&nbsp;93&nbsp;94&nbsp;95&nbsp;96&nbsp;97&nbsp;98&nbsp;99&nbsp;100 
	  </tt> </td>
	<td><tt> 2&nbsp;x&nbsp;2&nbsp;matrix<br>
	  43&nbsp;50<br>
	  53&nbsp;60 </tt></td>
	<td><tt> 2&nbsp;x&nbsp;10&nbsp;matrix<br>
	  41&nbsp;42&nbsp;43&nbsp;44&nbsp;45&nbsp;46&nbsp;47&nbsp;48&nbsp;49&nbsp;&nbsp;50<br>
	  91&nbsp;92&nbsp;93&nbsp;94&nbsp;95&nbsp;96&nbsp;97&nbsp;98&nbsp;99&nbsp;100 
	  </tt> </td>
	<td><tt> 10&nbsp;x&nbsp;2&nbsp;matrix<br>
	  &nbsp;4&nbsp;&nbsp;8<br>
	  14&nbsp;18<br>
	  24&nbsp;28<br>
	  34&nbsp;38<br>
	  44&nbsp;48<br>
	  54&nbsp;58<br>
	  64&nbsp;68<br>
	  74&nbsp;78<br>
	  84&nbsp;88<br>
	  94&nbsp;98 </tt> </td>
  </tr>
</table> 
@param matrix any matrix.
@param rowFraction the percentage of rows to be included in the view.
@param columnFraction the percentage of columns to be included in the view.
@param randomGenerator a uniform random number generator; set this parameter to <tt>null</tt> to use a default generator seeded with the current time.
@return the sampling view.
@throws IllegalArgumentException if <tt>! (0 <= rowFraction <= 1 && 0 <= columnFraction <= 1)</tt>.
@see cern.jet.random.sampling.RandomSampler
*/
public static DoubleMatrix2D viewSample(DoubleMatrix2D matrix, double rowFraction, double columnFraction, RandomEngine randomGenerator) {
	// check preconditions and allow for a little tolerance
	double epsilon = 1e-09;
	if (rowFraction < 0 - epsilon || rowFraction > 1 + epsilon) throw new IllegalArgumentException();
	if (rowFraction < 0) rowFraction = 0;
	if (rowFraction > 1) rowFraction = 1;

	if (columnFraction < 0 - epsilon || columnFraction > 1 + epsilon) throw new IllegalArgumentException();
	if (columnFraction < 0) columnFraction = 0;
	if (columnFraction > 1) columnFraction = 1;

	// random generator seeded with current time
	if (randomGenerator==null) randomGenerator = new cern.jet.random.engine.MersenneTwister((int) System.currentTimeMillis());

	int nrows = (int) Math.round(matrix.rows() * rowFraction);
	int ncols = (int) Math.round(matrix.columns() * columnFraction);
	int max = Math.max(nrows,ncols);
	long[] selected = new long[max]; // sampler works on long's, not int's

	// sample rows
	int n = nrows;
	int N = matrix.rows();
	cern.jet.random.sampling.RandomSampler.sample(n,N,n,0,selected,0,randomGenerator);
	int[] selectedRows = new int[n];
	for (int i=0; i<n; i++) selectedRows[i] = (int) selected[i];

	// sample columns
	n = ncols;
	N = matrix.columns();
	cern.jet.random.sampling.RandomSampler.sample(n,N,n,0,selected,0,randomGenerator);
	int[] selectedCols = new int[n];
	for (int i=0; i<n; i++) selectedCols[i] = (int) selected[i];

	return matrix.viewSelection(selectedRows, selectedCols);
}