if (state.getScaleFactor() <= 0) {
state.setScaleFactor(calculateScaleFactor(nextVector));
}
nextVector.assign(new Scale(1.0 / state.getScaleFactor()));
if (previousVector != null) {
nextVector.assign(previousVector, new PlusMult(-beta));
}
// now orthogonalize
double alpha = currentVector.dot(nextVector);
nextVector.assign(currentVector, new PlusMult(-alpha));
endTime(TimingSection.ITERATE);
startTime(TimingSection.ORTHOGANLIZE);
orthoganalizeAgainstAllButLast(nextVector, state);
endTime(TimingSection.ORTHOGANLIZE);
// and normalize
beta = nextVector.norm(2);
if (outOfRange(beta) || outOfRange(alpha)) {
log.warn("Lanczos parameters out of range: alpha = {}, beta = {}. Bailing out early!",
alpha, beta);
break;
}
nextVector.assign(new Scale(1 / beta));
state.setBasisVector(i, nextVector);
previousVector = currentVector;
currentVector = nextVector;
// save the projections and norms!
triDiag.set(i - 1, i - 1, alpha);
if (i < desiredRank - 1) {
triDiag.set(i - 1, i, beta);
triDiag.set(i, i - 1, beta);
}
state.setIterationNumber(++i);
}
startTime(TimingSection.TRIDIAG_DECOMP);
log.info("Lanczos iteration complete - now to diagonalize the tri-diagonal auxiliary matrix.");
// at this point, have tridiag all filled out, and basis is all filled out, and orthonormalized
EigenDecomposition decomp = new EigenDecomposition(triDiag);
Matrix eigenVects = decomp.getV();
Vector eigenVals = decomp.getRealEigenvalues();
endTime(TimingSection.TRIDIAG_DECOMP);
startTime(TimingSection.FINAL_EIGEN_CREATE);
for (int row = 0; row < i; row++) {
Vector realEigen = null;
// the eigenvectors live as columns of V, in reverse order. Weird but true.
Vector ejCol = eigenVects.viewColumn(i - row - 1);
int size = Math.min(ejCol.size(), state.getBasisSize());
for (int j = 0; j < size; j++) {
double d = ejCol.get(j);
Vector rowJ = state.getBasisVector(j);
if (realEigen == null) {
realEigen = rowJ.like();
}
realEigen.assign(rowJ, new PlusMult(d));
}
realEigen = realEigen.normalize();
state.setRightSingularVector(row, realEigen);
double e = eigenVals.get(row) * state.getScaleFactor();
if (!isSymmetric) {