U = new DenseDoubleMatrix2D(A.rows(), k);
V = new DenseDoubleMatrix2D(A.columns(), k);
seedingStrategy.seed(A, U, V);
// Temporary matrices
DoubleMatrix2D T = new DenseDoubleMatrix2D(k, k);
DoubleMatrix2D UT1 = new DenseDoubleMatrix2D(A.rows(), k);
DoubleMatrix2D UT2 = new DenseDoubleMatrix2D(A.rows(), k);
DoubleMatrix2D VT1 = new DenseDoubleMatrix2D(A.columns(), k);
DoubleMatrix2D VT2 = new DenseDoubleMatrix2D(A.columns(), k);
DoubleFunction plusEps = Functions.plus(eps);
if (stopThreshold >= 0)
{
updateApproximationError();
}
for (int i = 0; i < maxIterations; i++)
{
// Update V
U.zMult(U, T, 1, 0, true, false); // T <- U'U
A.zMult(U, VT1, 1, 0, true, false); // VT1 <- A'U
V.zMult(T, VT2, 1, 0, false, false); // VT2 <- VT
VT1.assign(plusEps); // TODO: shift this to the dividing function?
VT2.assign(plusEps);
VT1.assign(VT2, Functions.DIV); // VT1 <- VT1 ./ VT2
V.assign(VT1, Functions.MULT); // V <- V .* VT1
// Update U
V.zMult(V, T, 1, 0, true, false); // T <- V'V