Skip to content

Commit

Permalink
fixed files form Math #76
Browse files Browse the repository at this point in the history
  • Loading branch information
tdurieux committed Mar 7, 2017
1 parent d689957 commit 726d49f
Showing 1 changed file with 14 additions and 8 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -159,24 +159,27 @@ public RealMatrix getU()
if (m >= n) {
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
final RealMatrix e =
eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1);
final double[][] eData = e.getData();
final double[][] wData = new double[m][p];
double[] ei1 = eData[0];
for (int i = 0; i < p - 1; ++i) {
for (int i = 0; i < p; ++i) {
// compute W = B.E.S^(-1) where E is the eigenvectors matrix
final double mi = mainBidiagonal[i];
final double[] ei0 = ei1;
final double[] wi = wData[i];
if (i < n - 1) {
ei1 = eData[i + 1];
final double si = secondaryBidiagonal[i];
for (int j = 0; j < p; ++j) {
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
}
}
} else {
for (int j = 0; j < p; ++j) {
wData[p - 1][j] = ei1[j] * mainBidiagonal[p - 1] / singularValues[j];
wi[j] = mi * ei0[j] / singularValues[j];
}
}
}

for (int i = p; i < m; ++i) {
wData[i] = new double[p];
Expand Down Expand Up @@ -245,23 +248,26 @@ public RealMatrix getV()
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
// compute W = Bt.E.S^(-1) where E is the eigenvectors matrix
final RealMatrix e =
eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
final double[][] eData = e.getData();
final double[][] wData = new double[n][p];
double[] ei1 = eData[0];
for (int i = 0; i < p - 1; ++i) {
for (int i = 0; i < p; ++i) {
final double mi = mainBidiagonal[i];
final double[] ei0 = ei1;
final double[] wi = wData[i];
if (i < m - 1) {
ei1 = eData[i + 1];
final double si = secondaryBidiagonal[i];
for (int j = 0; j < p; ++j) {
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
}
}
} else {
for (int j = 0; j < p; ++j) {
wData[p - 1][j] = ei1[j] * mainBidiagonal[p - 1] / singularValues[j];
wi[j] = mi * ei0[j] / singularValues[j];
}
}
}
for (int i = p; i < n; ++i) {
wData[i] = new double[p];
}
Expand Down

0 comments on commit 726d49f

Please sign in to comment.