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

static DoubleMatrix1D diagMult(DoubleMatrix2D A, DoubleMatrix2D B){
    //Return the diagonal of the product of two matrices
    //A^T and B should have the same dimensions
    int m = A.rows();
    int n = A.columns();
    if(B.rows()!=n || B.columns()!=m){
        throw new Error("Size mismatch in diagMult: A is "+m+"x"+n+
                ", B is "+B.rows()+"x"+B.columns());
    }

    DoubleMatrix1D ans = DoubleFactory1D.dense.make(m);

    for(int i=0; i<m; i++){
        double s = 0;
        for(int k=0; k<n; k++)
            s += A.getQuick(i, k)*B.getQuick(k, i);
        ans.setQuick(i,s);
    }

    return ans;
}
 

private 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;
}
 

/**
 * Modifies the given covariance matrix to be a correlation matrix (in-place).
 * The correlation matrix is a square, symmetric matrix consisting of nothing but correlation coefficients.
 * The rows and the columns represent the variables, the cells represent correlation coefficients. 
 * The diagonal cells (i.e. the correlation between a variable and itself) will equal 1, for the simple reason that the correlation coefficient of a variable with itself equals 1. 
 * The correlation of two column vectors x and y is given by <tt>corr(x,y) = cov(x,y) / (stdDev(x)*stdDev(y))</tt> (Pearson's correlation coefficient).
 * A correlation coefficient varies between -1 (for a perfect negative relationship) to +1 (for a perfect positive relationship). 
 * See the <A HREF="http://www.cquest.utoronto.ca/geog/ggr270y/notes/not05efg.html"> math definition</A>
 * and <A HREF="http://www.stat.berkeley.edu/users/stark/SticiGui/Text/gloss.htm#correlation_coef"> another def</A>.
 * Compares two column vectors at a time. Use dice views to compare two row vectors at a time.
 * 
 * @param covariance a covariance matrix, as, for example, returned by method {@link #covariance(DoubleMatrix2D)}.
 * @return the modified covariance, now correlation matrix (for convenience only).
 */
public static DoubleMatrix2D correlation(DoubleMatrix2D covariance) {
	for (int i=covariance.columns(); --i >= 0; ) {
		for (int j=i; --j >= 0; ) {
			double stdDev1 = Math.sqrt(covariance.getQuick(i,i));
			double stdDev2 = Math.sqrt(covariance.getQuick(j,j));
			double cov = covariance.getQuick(i,j);
			double corr = cov / (stdDev1*stdDev2);
			
			covariance.setQuick(i,j,corr);
			covariance.setQuick(j,i,corr); // symmetric
		}
	}
	for (int i=covariance.columns(); --i >= 0; ) covariance.setQuick(i,i,1);

	return covariance;	
}
 

/**
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>non-negative</i> if <tt>A[i,j] &gt;= 0</tt> holds for all cells.
 * <p>
 * Note: Ignores tolerance.
 */
public boolean isNonNegative(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	for (int row = rows; --row >=0; ) {
		for (int column = columns; --column >= 0; ) {
			if (! (A.getQuick(row,column) >= 0)) 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 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;
}
 

/**
 * Returns the one-norm of matrix <tt>A</tt>, which is the maximum absolute column sum.
 */
public double norm1(DoubleMatrix2D A) {
	double max = 0;
	for (int column = A.columns(); --column >=0; ) {
		max = Math.max(max, norm1(A.viewColumn(column)));
	}
	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;
}
 

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

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

/**
 * A matrix <tt>A</tt> is <i>lower 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 isLowerBidiagonal(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;
}
 

/**
 * A matrix <tt>A</tt> is <i>lower 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 isLowerBidiagonal(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;
}
 

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

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 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) {}
}
 

static DoubleMatrix2D invertX(DoubleMatrix2D X){
        
        if(X.rows()!=X.columns())
            throw new RuntimeException("ERROR: Can't invert square matrix");
        
        if(X.rows()==1){//easy case of 1-D
            if(Math.abs(X.get(0,0))<1e-8)
                return DoubleFactory2D.dense.make(1,1,0);
            else
                return DoubleFactory2D.dense.make(1,1,1./X.get(0,0));
        }
    
        /*DoubleMatrix2D invX = null;
    
        try{
            invX = Algebra.DEFAULT.inverse(X);
        }
        catch(Exception e){
            System.err.println("ERROR: Singular matrix in MinVolEllipse calculation");
        }*/
        
        //Invert X, but if singular, this probably indicates a lack of points
        //but we can still make a meaningful lower-dimensional ellipsoid for the dimensions we have data
        //so we construct a pseudoinverse by setting those eigencomponents to 0
        //X is assumed to be symmetric (if not would need SVD instead...)
    
        EigenvalueDecomposition edec = new EigenvalueDecomposition(X);
        
        DoubleMatrix2D newD = edec.getD().copy();
        for(int m=0; m<X.rows(); m++){
            if( Math.abs(newD.get(m,m)) < 1e-8 ){
                //System.out.println("Warning: singular X in MinVolEllipse calculation");
                //this may come up somewhat routinely, especially with discrete DOFs in place
                newD.set(m,m,0);
            }
            else
                newD.set(m,m,1./newD.get(m,m));
        }

        DoubleMatrix2D invX = Algebra.DEFAULT.mult(edec.getV(),newD).zMult(edec.getV(), null, 1, 0, false, true);
        
        return invX;
}
 

public RobustEigenDecomposition(DoubleMatrix2D A, int maxIterations) throws ArithmeticException {

	Property.DEFAULT.checkSquare(A);

    this.maxIterations = maxIterations;

	n = A.columns();
	V = new double[n][n];
	d = new double[n];
	e = new double[n];

	issymmetric = Property.DEFAULT.isSymmetric(A);

	if (issymmetric) {
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				V[i][j] = A.getQuick(i,j);
			}
		}

		// Tridiagonalize.
		tred2();

		// Diagonalize.
		tql2();

	}
	else {
		H = new double[n][n];
		ort = new double[n];

		for (int j = 0; j < n; j++) {
			for (int i = 0; i < n; i++) {
				H[i][j] = A.getQuick(i,j);
			}
		}

		// Reduce to Hessenberg form.
		orthes();

		// Reduce Hessenberg to real Schur form.
		hqr2();
	}
}
 

/**
 * A matrix <tt>A</tt> is <i>square</i> if it has the same number of rows and columns.
 */
public boolean isSquare(DoubleMatrix2D A) {
	return A.rows() == A.columns();
}
 

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