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];
}
}

Dear take2rec,

You were kind enough not to call this a bug, which technically it isn't, because the routine is in the book subsection, "SVD of a Square Matrix". But it is clear that our intent in later subsections was that this routine should work for non-square matrices A, which, as written, it didn't!

Your fixed routine is absolutely correct. Thanks for submitting it.

Cheers,
Bill P.