--- vnl_svd.txx Thu Feb 21 04:45:27 2002 +++ /users/visics/vanroose/tmp/vnl_svd.txx Wed Mar 6 13:33:15 2002 @@ -190,8 +171,8 @@ last_tol_ = tol; rank_ = W_.n(); for (unsigned k = 0; k < W_.n(); k++) { singval_t& weight = W_(k, k); - if (vcl_abs(weight) <= tol) { + if (weight <= tol) { // no need for vcl_abs(weight) Winverse_(k,k) = 0; weight = 0; --rank_; @@ -277,37 +259,61 @@ return U_ * Winverse * V_.conjugate_transpose(); } +//: Solve the matrix equation M X = B, returning X +template +vnl_matrix vnl_svd::solve(vnl_matrix const& B, singval_t tol) const +{ + vnl_matrix x; // solution matrix + if (U_.rows() < U_.columns()) { // augment y with extra rows of + vnl_matrix yy(U_.rows(), B.columns(), T(0)); // zeros, so that it matches + yy.update(B); // cols of u.transpose. ??? + x = U_.conjugate_transpose() * yy; + } else + x = U_.conjugate_transpose() * B; + unsigned long i, j; + for (i = 0; i < x.rows(); i++) { // multiply with diagonal 1/W + singval_t weight = W_(i, i); + if (weight >= tol) + weight = singval_t(1) / weight; + else + weight = singval_t(0); + for (j = 0; j < x.columns(); j++) + x(i, j) *= weight; + } + x = V_ * x; // premultiply with v. + return x; +} //: Solve the matrix equation M X = B, returning X +// Assuming that the singular values W have been preinverted by the caller. template -vnl_matrix vnl_svd::solve(vnl_matrix const& B) const +vnl_matrix vnl_svd::solve_preinverted(vnl_matrix const& B, singval_t tol) const { vnl_matrix x; // solution matrix if (U_.rows() < U_.columns()) { // augment y with extra rows of vnl_matrix yy(U_.rows(), B.columns(), T(0)); // zeros, so that it matches yy.update(B); // cols of u.transpose. ??? x = U_.conjugate_transpose() * yy; } else x = U_.conjugate_transpose() * B; unsigned long i, j; - for (i = 0; i < x.rows(); i++) { // multiply with diagonal 1/W - T weight = W_(i, i); - if (weight != T(0)); //vnl_numeric_traits::zero) - weight = T(1) / weight; + for (i = 0; i < x.rows(); i++) { // multiply with diagonal W assumed inverted + // WAS: T weight = W_(i, i); if (weight != T(0)) weight = T(1) / weight; + singval_t weight = W_(i, i); for (j = 0; j < x.columns(); j++) x(i, j) *= weight; } x = V_ * x; // premultiply with v. return x; } //: Solve the matrix-vector system M x = y, returning x. template -vnl_vector vnl_svd::solve(vnl_vector const& y) const +vnl_vector vnl_svd::solve(vnl_vector const& y, singval_t tol) const { // fsm sanity check : if (y.size() != U_.rows()) { vcl_cerr << __FILE__ << ": size of rhs is incompatible with no. of rows in U_\n" << "y =" << y << "\n" << "m_=" << m_ << "\n" << "n_=" << n_ << "\n" @@ -331,17 +349,18 @@ x = U_.conjugate_transpose() * y; for (unsigned i = 0; i < x.size(); i++) { // multiply with diagonal 1/W - T weight = W_(i, i), zero_(0); - if (weight != zero_) + singval_t weight = W_(i, i), zero_(0); + if (weight <= tol) x[i] /= weight; else x[i] = zero_; } return V_ * x; // premultiply with v. } + template // FIXME. this should implement the above, not the other way round. -void vnl_svd::solve(T const *y, T *x) const { - solve(vnl_vector(y, m_)).copy_out(x); +void vnl_svd::solve(T const *y, T *x, singval_t tol) const { + solve(vnl_vector(y, m_), tol).copy_out(x); } //: Solve the matrix-vector system M x = y.