take2rec
12-23-2007, 11:04 AM
The code at page 71 for SVD::solve (second version for matrix parameters) only works for square matrices.
Below is a proposed enhanced version that works for non-square matrices. Comments welcome.
void SVD::solve(MatDoub_I &b, MatDoub_O &x, Doub thresh = -1.)
{
int i, j, p = b.ncols();
if (b.nrows() != m || x.nrows() != n || b.ncols() != x.ncols())
throw("SVD::solve bad sizes");
VecDoub xx(n), bcol(m);
for (j=0; j<p; j++) {
for (i=0; i<m; i++) bcol[i] = b[i][j];
solve(bcol, xx, thresh);
for (i=0; i<n; i++) x[i][j] = xx[i];
}
}
Below is a proposed enhanced version that works for non-square matrices. Comments welcome.
void SVD::solve(MatDoub_I &b, MatDoub_O &x, Doub thresh = -1.)
{
int i, j, p = b.ncols();
if (b.nrows() != m || x.nrows() != n || b.ncols() != x.ncols())
throw("SVD::solve bad sizes");
VecDoub xx(n), bcol(m);
for (j=0; j<p; j++) {
for (i=0; i<m; i++) bcol[i] = b[i][j];
solve(bcol, xx, thresh);
for (i=0; i<n; i++) x[i][j] = xx[i];
}
}