Java Code Examples for cern.colt.matrix.DoubleMatrix1D#getQuick()

public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	if (isUpperTriangular) A = A.viewDice();
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size() || size!=y.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
	}
	DoubleMatrix1D tmp = x.like();
	for (int i = 0; i < size; i++) {
		double sum = 0;
		for (int j = 0; j <= i; j++) {
			sum += A.getQuick(i,j) * x.getQuick(j);
		}
		for (int j = i + 1; j < size; j++) {
			sum += A.getQuick(j,i) * x.getQuick(j);
		}
		tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
	}
	y.assign(tmp);
}
 

public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	if (isUpperTriangular) A = A.viewDice();
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size() || size!=y.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
	}
	DoubleMatrix1D tmp = x.like();
	for (int i = 0; i < size; i++) {
		double sum = 0;
		for (int j = 0; j <= i; j++) {
			sum += A.getQuick(i,j) * x.getQuick(j);
		}
		for (int j = i + 1; j < size; j++) {
			sum += A.getQuick(j,i) * x.getQuick(j);
		}
		tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
	}
	y.assign(tmp);
}
 

/**
 * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
 * The result is <tt>true</tt> if <tt>A==B</tt>. 
 * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>, 
 * have the same size and 
 * <tt>! (Math.abs(A[i] - B[i]) > tolerance())</tt> holds for all indexes.
 * @param   A   the first matrix to compare.
 * @param   B   the second matrix to compare.
 * @return  <tt>true</tt> if both matrices are equal;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, DoubleMatrix1D B) {
	if (A==B) return true;
	if (! (A != null && B != null)) return false;
	int size = A.size();
	if (size != B.size()) return false;

	double epsilon = tolerance();
	for (int i=size; --i >= 0;) {
		//if (!(getQuick(i) == B.getQuick(i))) return false;
		//if (Math.abs(A.getQuick(i) - B.getQuick(i)) > epsilon) return false;
		double x = A.getQuick(i);
		double value = B.getQuick(i);
		double diff = Math.abs(value - x);
		if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
		if (!(diff <= epsilon)) return false;
	}
	return true;
}
 

/**
 * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
 * The result is <tt>true</tt> if <tt>A==B</tt>. 
 * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>, 
 * have the same size and 
 * <tt>! (Math.abs(A[i] - B[i]) > tolerance())</tt> holds for all indexes.
 * @param   A   the first matrix to compare.
 * @param   B   the second matrix to compare.
 * @return  <tt>true</tt> if both matrices are equal;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, DoubleMatrix1D B) {
	if (A==B) return true;
	if (! (A != null && B != null)) return false;
	int size = A.size();
	if (size != B.size()) return false;

	double epsilon = tolerance();
	for (int i=size; --i >= 0;) {
		//if (!(getQuick(i) == B.getQuick(i))) return false;
		//if (Math.abs(A.getQuick(i) - B.getQuick(i)) > epsilon) return false;
		double x = A.getQuick(i);
		double value = B.getQuick(i);
		double diff = Math.abs(value - x);
		if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
		if (!(diff <= epsilon)) return false;
	}
	return true;
}
 

/**
Sorts the matrix slices into ascending order, according to the <i>natural ordering</i> of the matrix values in the given <tt>[row,column]</tt> position.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort by other dimensions, use dice views. To sort descending, use flip views ...
<p>
The algorithm compares two 2-d slices at a time, determinining whether one is smaller, equal or larger than the other.
Comparison is based on the cell <tt>[row,column]</tt> within a slice.
Let <tt>A</tt> and <tt>B</tt> be two 2-d slices. Then we have the following rules
<ul>
<li><tt>A &lt;  B  iff A.get(row,column) &lt;  B.get(row,column)</tt>
<li><tt>A == B iff A.get(row,column) == B.get(row,column)</tt>
<li><tt>A &gt;  B  iff A.get(row,column) &gt;  B.get(row,column)</tt>
</ul>

@param matrix the matrix to be sorted.
@param row the index of the row inducing the order.
@param column the index of the column inducing the order.
@return a new matrix view having slices sorted by the values of the slice view <tt>matrix.viewRow(row).viewColumn(column)</tt>.
		<b>Note that the original matrix is left unaffected.</b>
@throws IndexOutOfBoundsException if <tt>row < 0 || row >= matrix.rows() || column < 0 || column >= matrix.columns()</tt>.
*/
public DoubleMatrix3D sort(DoubleMatrix3D matrix, int row, int column) {
	if (row < 0 || row >= matrix.rows()) throw new IndexOutOfBoundsException("row="+row+", matrix="+Formatter.shape(matrix));
	if (column < 0 || column >= matrix.columns()) throw new IndexOutOfBoundsException("column="+column+", matrix="+Formatter.shape(matrix));

	int[] sliceIndexes = new int[matrix.slices()]; // indexes to reorder instead of matrix itself
	for (int i=sliceIndexes.length; --i >= 0; ) sliceIndexes[i] = i;

	final DoubleMatrix1D sliceView = matrix.viewRow(row).viewColumn(column);
	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = sliceView.getQuick(a);
			double bv = sliceView.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(sliceIndexes,0,sliceIndexes.length,comp);

	// view the matrix according to the reordered slice indexes
	// take all rows and columns in the original order
	return matrix.viewSelection(sliceIndexes,null,null);
}
 

/**
 * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
 * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
 * <tt>! (Math.abs(value - A[i]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   value   the value to compare against.
 * @return  <tt>true</tt> if the matrix is equal to the value;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, double value) {
	if (A==null) return false;
	double epsilon = tolerance();
	for (int i = A.size(); --i >= 0;) {
		//if (!(A.getQuick(i) == value)) return false;
		//if (Math.abs(value - A.getQuick(i)) > epsilon) return false;
		double x = A.getQuick(i);
		double diff = Math.abs(value - x);
		if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
		if (!(diff <= epsilon)) return false;
	}
	return true;
}
 

public void dger(double alpha, DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) {
	cern.jet.math.PlusMult fun = cern.jet.math.PlusMult.plusMult(0);
	for (int i=A.rows(); --i >= 0; ) {
		fun.multiplicator = alpha * x.getQuick(i);
 		A.viewRow(i).assign(y,fun);
		
	}
}
 

public FactorizationUserIntentModel(DoubleMatrix1D userVector) {
    Set<Integer> nonZeroFidx = new HashSet<>();
    for (int i = 0; i < userVector.size(); i++) {
        if (userVector.getQuick(i) > 0) {
            nonZeroFidx.add(i);
        }
    }
    this.userVector = userVector;
    this.nonZeroFactors = nonZeroFidx.stream()
            .map(featureData::fidx2feature)
            .collect(Collectors.toSet());
}
 

public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
	if (transposeA) {
		A = A.viewDice();
		isUpperTriangular = !isUpperTriangular;
	}
	
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
	}
	    
	DoubleMatrix1D b = x.like();
	DoubleMatrix1D y = x.like();
	if (isUnitTriangular) {
		y.assign(1);
	}
	else {
		for (int i = 0; i < size; i++) {
			y.setQuick(i, A.getQuick(i,i));
		}
	}
	
	for (int i = 0; i < size; i++) {
		double sum = 0;
		if (!isUpperTriangular) {
			for (int j = 0; j < i; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);
			}
			sum += y.getQuick(i) * x.getQuick(i);
		}
		else {
			sum += y.getQuick(i) * x.getQuick(i);
			for (int j = i + 1; j < size; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);			}
		}
		b.setQuick(i,sum);
	}
	x.assign(b);
}
 

public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
	if (transposeA) {
		A = A.viewDice();
		isUpperTriangular = !isUpperTriangular;
	}
	
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
	}
	    
	DoubleMatrix1D b = x.like();
	DoubleMatrix1D y = x.like();
	if (isUnitTriangular) {
		y.assign(1);
	}
	else {
		for (int i = 0; i < size; i++) {
			y.setQuick(i, A.getQuick(i,i));
		}
	}
	
	for (int i = 0; i < size; i++) {
		double sum = 0;
		if (!isUpperTriangular) {
			for (int j = 0; j < i; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);
			}
			sum += y.getQuick(i) * x.getQuick(i);
		}
		else {
			sum += y.getQuick(i) * x.getQuick(i);
			for (int j = i + 1; j < size; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);			}
		}
		b.setQuick(i,sum);
	}
	x.assign(b);
}
 

/**
 * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
 * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
 * <tt>! (Math.abs(value - A[i]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   value   the value to compare against.
 * @return  <tt>true</tt> if the matrix is equal to the value;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, double value) {
	if (A==null) return false;
	double epsilon = tolerance();
	for (int i = A.size(); --i >= 0;) {
		//if (!(A.getQuick(i) == value)) return false;
		//if (Math.abs(value - A.getQuick(i)) > epsilon) return false;
		double x = A.getQuick(i);
		double diff = Math.abs(value - x);
		if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
		if (!(diff <= epsilon)) return false;
	}
	return true;
}
 

private static void doRR1(int L, DoubleMatrix1D w, double[][] x, double[] y, double[] c, double lambda) {
    int N = x.length;
    int K = x[0].length;
    
    double[] e = new double[N];
    for (int i = 0; i < N; i++) {
        double pred = 0.0;
        for (int k = 0; k < K; k++) {
            pred += w.getQuick(k) * x[i][k];
        }
        e[i] = y[i] - pred;
    }

    for (int l = 0; l < L; l++) {
        for (int k = 0; k < K; k++) {
            for (int i = 0; i < N; i++) {
                e[i] += w.getQuick(k) * x[i][k];
            }
            double a = 0.0;
            double d = 0.0;
            for (int i = 0; i < N; i++) {
                a += c[i] * x[i][k] * x[i][k];
                d += c[i] * x[i][k] * e[i];
            }
            w.setQuick(k, d / (lambda + a));
            for (int i = 0; i < N; i++) {
                e[i] -= w.getQuick(k) * x[i][k];
            }
        }
    }

}
 

public void dger(double alpha, DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) {
	cern.jet.math.PlusMult fun = cern.jet.math.PlusMult.plusMult(0);
	for (int i=A.rows(); --i >= 0; ) {
		fun.multiplicator = alpha * x.getQuick(i);
 		A.viewRow(i).assign(y,fun);
		
	}
}
 

public PLSAUserIntentModel(DoubleMatrix1D userVector) {
    this.nonZeroFactors = new HashSet<>();
    for (int i = 0; i < userVector.size(); i++) {
        if (userVector.getQuick(i) > 0) {
            nonZeroFactors.add(i);
        }
    }
    this.userVector = userVector;
}
 

/**
Sorts the matrix rows into ascending order, according to the <i>natural ordering</i> of the matrix values in the given column.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort columns by rows, use dice views. To sort descending, use flip views ...
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt>4 x 2 matrix: <br>
	  7, 6<br>
	  5, 4<br>
	  3, 2<br>
	  1, 0 <br>
	  </tt></td>
	<td align="left" valign="top"> 
	  <p><tt>column = 0;<br>
		view = quickSort(matrix,column);<br>
		System.out.println(view); </tt><tt><br>
		==> </tt></p>
	  </td>
	<td valign="top"> 
	  <p><tt>4 x 2 matrix:<br>
		1, 0<br>
		3, 2<br>
		5, 4<br>
		7, 6</tt><br>
		The matrix IS NOT SORTED.<br>
		The new VIEW IS SORTED.</p>
	  </td>
  </tr>
</table>

@param matrix the matrix to be sorted.
@param column the index of the column inducing the order.
@return a new matrix view having rows sorted by the given column.
		<b>Note that the original matrix is left unaffected.</b>
@throws IndexOutOfBoundsException if <tt>column < 0 || column >= matrix.columns()</tt>.
*/
public DoubleMatrix2D sort(DoubleMatrix2D matrix, int column) {
	if (column < 0 || column >= matrix.columns()) throw new IndexOutOfBoundsException("column="+column+", matrix="+Formatter.shape(matrix));

	int[] rowIndexes = new int[matrix.rows()]; // row indexes to reorder instead of matrix itself
	for (int i=rowIndexes.length; --i >= 0; ) rowIndexes[i] = i;

	final DoubleMatrix1D col = matrix.viewColumn(column);
	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = col.getQuick(a);
			double bv = col.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(rowIndexes,0,rowIndexes.length,comp);

	// view the matrix according to the reordered row indexes
	// take all columns in the original order
	return matrix.viewSelection(rowIndexes,null);
}
 

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

/**
Sorts the vector into ascending order, according to the <i>natural ordering</i>.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort descending, use flip views ...
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt> 7, 1, 3, 1<br>
	  </tt></td>
	<td valign="top"> 
	  <p><tt> ==&gt; 1, 1, 3, 7<br>
		The vector IS NOT SORTED.<br>
		The new VIEW IS SORTED.</tt></p>
	</td>
  </tr>
</table>

@param vector the vector to be sorted.
@return a new sorted vector (matrix) view. 
		<b>Note that the original matrix is left unaffected.</b>
*/
public DoubleMatrix1D sort(final DoubleMatrix1D vector) {
	int[] indexes = new int[vector.size()]; // row indexes to reorder instead of matrix itself
	for (int i=indexes.length; --i >= 0; ) indexes[i] = i;

	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = vector.getQuick(a);
			double bv = vector.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(indexes,0,indexes.length,comp);

	return vector.viewSelection(indexes);
}
 

/**
Sorts the vector into ascending order, according to the <i>natural ordering</i>.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort descending, use flip views ...
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt> 7, 1, 3, 1<br>
	  </tt></td>
	<td valign="top"> 
	  <p><tt> ==&gt; 1, 1, 3, 7<br>
		The vector IS NOT SORTED.<br>
		The new VIEW IS SORTED.</tt></p>
	</td>
  </tr>
</table>

@param vector the vector to be sorted.
@return a new sorted vector (matrix) view. 
		<b>Note that the original matrix is left unaffected.</b>
*/
public DoubleMatrix1D sort(final DoubleMatrix1D vector) {
	int[] indexes = new int[vector.size()]; // row indexes to reorder instead of matrix itself
	for (int i=indexes.length; --i >= 0; ) indexes[i] = i;

	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = vector.getQuick(a);
			double bv = vector.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(indexes,0,indexes.length,comp);

	return vector.viewSelection(indexes);
}
 

/**
Sorts the matrix rows into ascending order, according to the <i>natural ordering</i> of the matrix values in the given column.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort columns by rows, use dice views. To sort descending, use flip views ...
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt>4 x 2 matrix: <br>
	  7, 6<br>
	  5, 4<br>
	  3, 2<br>
	  1, 0 <br>
	  </tt></td>
	<td align="left" valign="top"> 
	  <p><tt>column = 0;<br>
		view = quickSort(matrix,column);<br>
		System.out.println(view); </tt><tt><br>
		==> </tt></p>
	  </td>
	<td valign="top"> 
	  <p><tt>4 x 2 matrix:<br>
		1, 0<br>
		3, 2<br>
		5, 4<br>
		7, 6</tt><br>
		The matrix IS NOT SORTED.<br>
		The new VIEW IS SORTED.</p>
	  </td>
  </tr>
</table>

@param matrix the matrix to be sorted.
@param column the index of the column inducing the order.
@return a new matrix view having rows sorted by the given column.
		<b>Note that the original matrix is left unaffected.</b>
@throws IndexOutOfBoundsException if <tt>column < 0 || column >= matrix.columns()</tt>.
*/
public DoubleMatrix2D sort(DoubleMatrix2D matrix, int column) {
	if (column < 0 || column >= matrix.columns()) throw new IndexOutOfBoundsException("column="+column+", matrix="+Formatter.shape(matrix));

	int[] rowIndexes = new int[matrix.rows()]; // row indexes to reorder instead of matrix itself
	for (int i=rowIndexes.length; --i >= 0; ) rowIndexes[i] = i;

	final DoubleMatrix1D col = matrix.viewColumn(column);
	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = col.getQuick(a);
			double bv = col.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(rowIndexes,0,rowIndexes.length,comp);

	// view the matrix according to the reordered row indexes
	// take all columns in the original order
	return matrix.viewSelection(rowIndexes,null);
}
 

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