luRows[i] = lu.viewRow(i);
}
IntArrayList nonZeroIndexes =
new IntArrayList(); // sparsity
DoubleMatrix1D luColj = lu.viewColumn(0).like(); // blocked column j
Mult multFunction = Mult.mult(0);
// Outer loop.
int cutOff = 10;
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 / cutOff; // == 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;
}
}