[Vxl-maintainers] [ vxl-Feature Requests-523333 ] some pbs wit vnl_svd

[<no...@so...>]

Feature Requests item #523333, was opened at 2002-02-27 06:03
You can respond by visiting: 
http://sourceforge.net/tracker/?func=detail&atid=381096&aid=523333&group_id=24293

Category: Interface Improvements (example)
Group: None
Status: Open
Priority: 5
Submitted By: lapresté jean-thierry (lapreste)
Assigned to: Nobody/Anonymous (nobody)
Summary: some pbs wit vnl_svd

Initial Comment:
I have found no real bugs in svd, but I tried

 #define C_Matrix vnl_matrix < double > 

 I have a non square X and I do

  vnl_svd < double > svd(X,  1.0E-8);
  C_Matrix v_sol = svd.solve(Y);

And that doesn't worked fine v_sol being filled
with inf 

I tried 

 C_Matrix inv =  svd.pinverse();
 C_Matrix v_sol = inv*Y; 

and that worked fine

So I have investigated the two files vnl_svd.h and
vnl_svd.txx and propose some modifications that
do not change the interface and a new function


Now I can write :
  
  vnl_svd < double > svd(X);
  C_Matrix v_sol = svd.solve(Y,  1.0E-8);

and it works. Perhaps I missed something, however
I join my versions of the files and would like to know
how to
participate to further devlopments or corrections if I
find
any. 

   sincerely

    J.T. Lapresté

lap...@la...

------------------------------------------------------------------
------------------------------------------------------------------
here vnl_svd.h 
------------------------------------------------------------------
#ifndef vnl_svd_h_
#define vnl_svd_h_
#ifdef __GNUC__
#pragma interface
#endif

//:
//  \file
//  \brief Holds the singular value decomposition of a
vnl_matrix.
//  \author Andrew W. Fitzgibbon, Oxford IERG
//  \date   15 Jul 96
//
// \verbatim
//  Modifications
// F. Schaffalitzky, Oxford IESRG, 26 Mar 1999
//     1. The singular values are now stored as reals
(not complexes) when T is complex.
//     2. Fixed bug : for complex T, matrices have to
be conjugated as well as transposed.
//   Feb.2002 - Peter Vanroose - brief doxygen comment
placed on single line
// J.T. Laprest, Lasmea Clermont-Ferrand,  27-Feb-2002
//     1. The typename singval_t must be used in solve,
it allows to add a supplementary 
//        tolerance parameter and mahe svdsolve
consistent with pinverse 
//     2. The modified interfaces are
//        vnl_matrix<T> vnl_svd<T>::solve(vnl_matrix<T>
const& B, singval_t tol)  const;
//        vnl_vector<T> solve (vnl_vector<T> const& y,
singval_t tol=singval_t(0)) const;
//        void          solve (T const *rhs, T *lhs, 
singval_t tol=singval_t(0)) const; 
//     3. The routine 
//        vnl_matrix<T>
vnl_svd<T>::solve_preinverted(vnl_matrix<T> const& B,
singval_t tol)  const;
//        has been added 
// \endverbatim

#include <vnl/vnl_numeric_traits.h>
#include <vnl/vnl_vector.h>
#include <vnl/vnl_matrix.h>
#include <vnl/vnl_diag_matrix.h>
#include <vcl_iosfwd.h>

//: Holds the singular value decomposition of a
vnl_matrix.
//
//  The class holds three matrices U, W, V such that
the original matrix
//  $M = U W V^\top$.  The DiagMatrix W stores the
singular values in decreasing
//  order.  The columns of U which correspond to the
nonzero singular values
//  form a basis for range of M, while the columns of V
corresponding to the
//  zero singular values are the nullspace.
//
//  The SVD is computed at construction time, and
enquiries may then be made
//  of the SVD.  In particular, this allows easy access
to multiple
//  right-hand-side solves without the bother of
putting all the RHS's into a
//  Matrix.
//
//  This class is supplied even though there is an
existing vnl_matrix method
//  for several reasons:
//
//  It is more convenient to use as it manages all the
storage for
//  the U,S,V matrices, allowing repeated queries of
the same SVD
//  results.
//
//  It avoids namespace clutter in the Matrix class.  
While svd()
//  is a perfectly reasonable method for a Matrix,
there are many other
//  decompositions that might be of interest, and
adding them all would
//  make for a very large Matrix class.
//
//  It demonstrates the holder model of compute class,
implementing an
//  algorithm on an object without adding a member that
may not be of
//  general interest.  A similar pattern can be used
for other
//  decompositions which are not defined as members of
the library Matrix
//  class.
//
//  It extends readily to n-ary operations, such as
generalized
//  eigensystems, which cannot be members of just one
matrix.

export template <class T>
class vnl_svd {
public:
  //: The singular values of a matrix of complex<T> are
of type T, not complex<T>
  typedef typename vnl_numeric_traits<T>::abs_t
singval_t;

  //:
  // Construct an vnl_svd<T> object from $m \times n$
matrix $M$.  The
  // vnl_svd<T> object contains matrices $U$, $W$, $V$
such that
  // $U W V^\top = M$.
  //
  // Uses linpack routine DSVDC to calculate an
``economy-size'' SVD
  // where the returned $U$ is the same size as $M$,
while $W$
  // and $V$ are both $n \times n$.  This is efficient
for
  // large rectangular solves where $m > n$, typical in
least squares.
  //
  // The optional argument zero_out_tol is used to mark
the zero singular
  // values: If nonnegative, any s.v. smaller than
zero_out_tol in
// absolute value is set to zero.  If zero_out_tol is
negative, the
  // zeroing is relative to |zero_out_tol| *
sigma_max();

  vnl_svd(vnl_matrix<T> const &M, double zero_out_tol =
0.0);
 ~vnl_svd() {}

  // Data
Access---------------------------------------------------------------

  //: find weights below threshold tol, zero them out,
and update W_ and Winverse_
  void            zero_out_absolute(double tol = 1e-8);
//sqrt(machine epsilon)

  //: find weights below tol*max(w) and zero them out
  void            zero_out_relative(double tol = 1e-8);
//sqrt(machine epsilon)
  int             singularities () const { return
W_.n() - rank(); }
  int             rank () const { return rank_; }
  singval_t       well_condition () const { return
sigma_min()/sigma_max(); }

    //: Calculate determinant as product of diagonals
in W.
  singval_t       determinant_magnitude () const;
  singval_t       norm() const;

  //: Return the matrix U.
  vnl_matrix<T>      & U()       { return U_; }

  //: Return the matrix U.
  vnl_matrix<T> const& U() const { return U_; }

  //: Return the matrix U's (i,j)th entry (to avoid
svd.U()(i,j); ).
  T U(int i, int j) { return U_(i,j); }

  //: Get at DiagMatrix (q.v.) of singular values,
sorted from largest to smallest
  vnl_diag_matrix<singval_t>       & W()             {
return W_; }

  //: Get at DiagMatrix (q.v.) of singular values,
sorted from largest to smallest
  vnl_diag_matrix<singval_t> const & W() const       {
return W_; }
  vnl_diag_matrix<singval_t>       & Winverse()      {
return Winverse_; }
  vnl_diag_matrix<singval_t> const & Winverse() const {
return Winverse_; }
  singval_t                   & W(int i, int j) {
return W_(i,j); }
  singval_t                   & W(int i)        {
return W_(i,i); }
  singval_t     sigma_max() const { return W_(0,0);
}       // largest
  singval_t     sigma_min() const { return
W_(n_-1,n_-1); } // smallest

  //: Return the matrix V.
  vnl_matrix<T>      & V()       { return V_; }

  //: Return the matrix V.
  vnl_matrix<T> const& V() const { return V_; }

  //: Return the matrix V's (i,j)th entry (to avoid
svd.V()(i,j); ).
  T V(int i, int j) { return V_(i,j); }

  //:
  vnl_matrix<T> inverse () const;

  //: pseudo-inverse (for non-square matrix).
  vnl_matrix<T> pinverse () const;

  //: Calculate inverse of transpose.
  vnl_matrix<T> tinverse () const;

  //: Recompose SVD to U*W*V'
  vnl_matrix<T> recompose () const;

  //: Solve the matrix equation M X = B, returning X
  vnl_matrix<T> solve (vnl_matrix<T> const& B,
singval_t tol =  singval_t(0)) const;

  //: Solve the matrix equation M X = B.
  // Assuming that the singular values W have been
preinverted by the caller.
  vnl_matrix<T>
vnl_svd<T>::solve_preinverted(vnl_matrix<T> const& B,
singval_t tol)  const;

  //: Solve the matrix-vector system M x = y, returning
x.
  vnl_vector<T> solve (vnl_vector<T> const& y,
singval_t tol=singval_t(0)) const;
  void          solve (T const *rhs, T *lhs,  singval_t
tol=singval_t(0)) const; // min ||A*lhs - rhs||

  //: Solve the matrix-vector system M x = y.
  // Assuming that the singular values W have been
preinverted by the caller.
  void solve_preinverted(vnl_vector<T> const& rhs,
vnl_vector<T>* out) const;

  //: Return N such that M * N = 0
  vnl_matrix<T> nullspace() const;

  //: Return N such that M' * N = 0
  vnl_matrix<T> left_nullspace() const;

  //: Return N such that M * N = 0
  vnl_matrix<T> nullspace(int
required_nullspace_dimension) const;

  //: Implementation to be done yet; currently returns
left_nullspace(). - PVR.
  vnl_matrix<T> left_nullspace(int
required_nullspace_dimension) const;

  //: Return the rightmost column of V.
  //  Does not check to see whether or not the matrix
actually was rank-deficient -
  // the caller is assumed to have examined W and
decided that to his or her satisfaction.
  vnl_vector<T> nullvector() const;

  //: Return the rightmost column of U.
  //  Does not check to see whether or not the matrix
actually was rank-deficient.
  vnl_vector<T> left_nullvector() const;

private:

  int m_, n_;              // Size of M, local cache.
  vnl_matrix<T> U_;        // Columns Ui are basis for
range of M for Wi != 0
  vnl_diag_matrix<singval_t> W_;// Singular values,
sorted in decreasing order
  vnl_diag_matrix<singval_t> Winverse_;
  vnl_matrix<T> V_;       // Columns Vi are basis for
nullspace of M for Wi = 0
  unsigned rank_;
  bool have_max_;
  singval_t max_;
  bool have_min_;
  singval_t min_;
  double last_tol_;

  // Disallow assignment.
  vnl_svd(vnl_svd<T> const &) { }
  vnl_svd<T>& operator=(vnl_svd<T> const &) { return
*this; }
};

template <class T>
inline
vnl_matrix<T> vnl_svd_inverse(vnl_matrix<T> const& m)
{
  return vnl_svd<T>(m).inverse();
}

// this aint no friend.
export template <class T>
vcl_ostream& operator<<(vcl_ostream&, vnl_svd<T> const&
svd);

#endif // vnl_svd_h_


--------------------------------------------------------------------
here vnl_svd.txx
--------------------------------------------------------------------
#ifndef vnl_svd_txx_
#define vnl_svd_txx_

//:
// \file

#include "vnl_svd.h"

#include <vcl_cstdlib.h> // vcl_abort()
#include <vcl_cassert.h>
#include <vcl_complex.h>
#include <vcl_fstream.h>
#include <vcl_algorithm.h> // min

#include <vnl/vnl_math.h>
#include <vnl/vnl_fortran_copy.h>
#include <vnl/algo/vnl_netlib.h>

#ifdef HAS_FSM_PACK
template <typename T> int
fsm_svdc_cxx(vnl_netlib_svd_proto(T));
# define vnl_linpack_svdc fsm_svdc_cxx
#else
// use C++ overloading to call the right linpack
routine from the template code :
#define macro(p, T) \
inline int vnl_linpack_svdc(vnl_netlib_svd_proto(T)) \
{ return p##svdc_(vnl_netlib_svd_params); }
macro(s, float);
macro(d, double);
macro(c, vcl_complex<float>);
macro(z, vcl_complex<double>);
#undef macro
#endif

//--------------------------------------------------------------------------------

static bool test_heavily = false;

template <class T>
vnl_svd<T>::vnl_svd(vnl_matrix<T> const& M, double
zero_out_tol):
  m_(M.rows()),
  n_(M.columns()),
  U_(m_, n_),
  W_(n_),
  Winverse_(n_),
  V_(n_, n_)
{
  assert(m_ > 0);
  assert(n_ > 0);

  {
    int n = M.rows();
    int p = M.columns();
    int mm = vcl_min(n+1,p);

    // Copy source matrix into fortran storage
    // SVD is slow, don't worry about the cost of this
transpose.
    vnl_fortran_copy<T> X(M);

    // Make workspace vectors.
    vnl_vector<T> work(n, T(0));
    vnl_vector<T> uspace(n*p, T(0));
    vnl_vector<T> vspace(p*p, T(0));
    vnl_vector<T> wspace(mm, T(0)); // complex fortran
routine actually _wants_ complex W!
    vnl_vector<T> espace(p, T(0));

    // Call Linpack SVD
    int info = 0;
    vnl_linpack_svdc((T*)X, n, n, p,
                     wspace.data_block(),
                     espace.data_block(),
                     uspace.data_block(), n,
                     vspace.data_block(), p,
                     work.data_block(),
                     21, &info);

    // Error return?
    if (info != 0) {
      // If info is non-zero, it contains the number of
singular values
      // for this the SVD algorithm failed to converge.
The condition is
      // not bogus. Even if the returned singular
values are sensible, 
      // the singular vectors can be utterly wrong.

      // It is possible the failure was due to NaNs or
infinities in the
      // matrix. Check for that now.
      M.assert_finite();

      // If we get here it might be because the scalar
type has such 
      // extreme precision that too few iterations were
performed to
      // converge to within machine precision (that is
the svdc criterion).
      // The only solution to that is to increase the
maximum number of
      // iterations in the netlib code. Diagnose the
problem here by
      // printing a warning message.
      vcl_cerr << __FILE__ ": suspicious return value
(" << info << ") from SVDC\n"
               << __FILE__ ": M is " << M.rows() << 'x'
<< M.cols() << vcl_endl;
    }

    // Copy fortran outputs into our storage
    {
      const T *d = uspace.data_block();
      for(int j = 0; j < p; ++j)
        for(int i = 0; i < n; ++i)
          U_(i,j) = *d++;
    }

    for(int j = 0; j < mm; ++j)
      W_(j, j) = vcl_abs(wspace(j)); // we get rid of
complexness here.

    for(int j = mm; j < n_; ++j)
      W_(j, j) = 0;

    {
      const T *d = vspace.data_block();
      for(int j = 0; j < p; ++j)
        for(int i = 0; i < p; ++i)
          V_(i,j) = *d++;
    }
  }

  if (test_heavily) {
    // Test that recomposed matrix == M
    typedef typename vnl_numeric_traits<T>::abs_t
abs_t;
    abs_t recomposition_residual = vcl_abs((recompose()
- M).fro_norm());
    abs_t n = vcl_abs(M.fro_norm());
    abs_t thresh = m_ * vnl_math::eps * n;
    if (recomposition_residual > thresh) {
      vcl_cerr << "vnl_svd<T>::vnl_svd<T>() -- Warning,
recomposition_residual = "
               << recomposition_residual << vcl_endl
               << "fro_norm(M) = " << n << vcl_endl
               << "eps*fro_norm(M) = " << thresh <<
vcl_endl
               << "Press return to continue\n";
      char x;
      vcl_cin.get(&x, 1, '\n');
    }
  }

  if (zero_out_tol >= 0)
    // Zero out small sv's and update rank count.
    zero_out_absolute(double(+zero_out_tol));
  else
    // negative tolerance implies relative to max elt.
    zero_out_relative(double(-zero_out_tol));
}

#if 0
// Assignment
template <class T>
vnl_svd<T>& vnl_svd<T>::operator=(vnl_svd<T> const&
that)
{
  U_ = that.U_;
  W_ = that.W_;
  Winverse_ = that.Winverse_;
  V_ = that.V_;
  rank_ = that.rank_;
  return *this;
}
#endif

template <class T>
vcl_ostream& operator<<(vcl_ostream& s, const
vnl_svd<T>& svd)
{
  s << "vnl_svd<T>:\n"
//  << "M = [\n" << M << "]\n"
    << "U = [\n" << svd.U() << "]\n"
    << "W = " << svd.W() << "\n"
    << "V = [\n" << svd.V() << "]\n"
    << "rank = " << svd.rank() << vcl_endl;
  return s;
}

//-----------------------------------------------------------------------------
// Chunky bits.

//: find weights below threshold tol, zero them out,
and update W_ and Winverse_
template <class T>
void
vnl_svd<T>::zero_out_absolute(double tol)
{
  last_tol_ = tol;
  rank_ = W_.n();
  for (unsigned k = 0; k < W_.n(); k++) {
    singval_t & weight = W_(k, k);
    if (weight <= tol) {                    // no need
for vcl_abs(weight)
      Winverse_(k,k) = 0;
      weight = 0;
      --rank_;
    } else {
      Winverse_(k,k) = singval_t(1.0)/weight;
    }
  }
}

//: find weights below tol*max(w) and zero them out
template <class T> void
vnl_svd<T>::zero_out_relative(double tol) //
sqrt(machine epsilon)
{
  zero_out_absolute(tol * vcl_abs(sigma_max()));
}


//: Calculate determinant as product of diagonals in W.
template <class T>
vnl_svd<T>::singval_t
vnl_svd<T>::determinant_magnitude() const
{
  {
    static bool warned = false;
    if (!warned && m_ != n_) {
      vcl_cerr << __FILE__ ": called
determinant_magnitude() on SVD of non-square matrix" <<
vcl_endl;
      warned = true;
    }
  }
  singval_t product = W_(0, 0);
  for (unsigned long k = 1; k < W_.columns(); k++)
    product *= W_(k, k);

  return product;
}

template <class T>
vnl_svd<T>::singval_t vnl_svd<T>::norm() const
{
  return vcl_abs(sigma_max());
}

//: Recompose SVD to U*W*V'
template <class T>
vnl_matrix<T> vnl_svd<T>::recompose() const
{
  vnl_matrix<T> W(W_.rows(),W_.columns());
  W.fill(T(0));
  for (unsigned i=0;i<rank_;i++)
    W(i,i)=W_(i,i);

  return U_*W*V_.conjugate_transpose();
}


template <class T>
vnl_matrix<T> vnl_svd<T>::inverse() const
{
  return pinverse();
}


//: Calculate pseudo-inverse.
template <class T>
vnl_matrix<T> vnl_svd<T>::pinverse()  const
{
  vnl_matrix<T>
Winverse(Winverse_.rows(),Winverse_.columns());
  Winverse.fill(T(0));
  for (unsigned i=0;i<rank_;i++)
    Winverse(i,i)=Winverse_(i,i);

  return V_ * Winverse * U_.conjugate_transpose();
}


//: Calculate inverse of transpose.
template <class T>
vnl_matrix<T> vnl_svd<T>::tinverse()  const
{
  vnl_matrix<T>
Winverse(Winverse_.rows(),Winverse_.columns());
  Winverse.fill(T(0));
  for (unsigned i=0;i<rank_;i++)
    Winverse(i,i)=Winverse_(i,i);

  return U_ * Winverse * V_.conjugate_transpose();
}


//: Solve the matrix equation M X = B, returning X
template <class T>
vnl_matrix<T> vnl_svd<T>::solve(vnl_matrix<T> const& B,
singval_t tol)  const
{
  vnl_matrix<T> x;                                     
// solution matrix
  if (U_.rows() < U_.columns()) {                      
// augment y with extra rows of
    vnl_matrix<T> 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){; //vnl_numeric_traits<T>::zero)
      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 <class T>
vnl_matrix<T>
vnl_svd<T>::solve_preinverted(vnl_matrix<T> const& B,
singval_t tol)  const
{
  vnl_matrix<T> x;                                     
// solution matrix
  if (U_.rows() < U_.columns()) {                      
// augment y with extra rows of
    vnl_matrix<T> 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 W assumed inverted
    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 <class T>
vnl_vector<T> vnl_svd<T>::solve(vnl_vector<T> 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"
             << "U_=\n" << U_
             << "V_=\n" << V_
             << "W_=\n" << W_;
  }

  vnl_vector<T> x(V_.rows());                   //
Solution matrix.
  if (U_.rows() < U_.columns()) {               //
Augment y with extra rows of
    vnl_vector<T> yy(U_.rows(), T(0));          //
zeros, so that it matches
    if (yy.size()<y.size()) { // fsm
      vcl_cerr << "yy=" << yy << vcl_endl
               << "y =" << y  << vcl_endl;
      // the update() call on the next line will
abort...
    }
    yy.update(y);                               // cols
of u.transpose.
    x = U_.conjugate_transpose() * yy;
  }
  else
    x = U_.conjugate_transpose() * y;

  for (unsigned i = 0; i < x.size(); i++) {        //
multiply with diagonal 1/W
    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 <class T> // FIXME. this should implement the
above, not the other way round.
void vnl_svd<T>::solve(T const *y, T *x,  singval_t
tol) const {
  solve(vnl_vector<T>(y, m_), tol).copy_out(x);
}

//: Solve the matrix-vector system M x = y.
// Assume that the singular values W have been
preinverted by the caller.
template <class T>
void vnl_svd<T>::solve_preinverted(vnl_vector<T> const&
y, vnl_vector<T>* x_out)  const
{
  vnl_vector<T> x;              // solution matrix
  if (U_.rows() < U_.columns()) {               //
augment y with extra rows of
    vcl_cout << "vnl_svd<T>::solve_preinverted() --
Augmenting y\n";
    vnl_vector<T> yy(U_.rows(), T(0));     // zeros, so
that it match
    yy.update(y);                               // cols
of u.transpose. ??
    x = U_.conjugate_transpose() * yy;
  } else
    x = U_.conjugate_transpose() * y;
  for (unsigned i = 0; i < x.size(); i++)  // multiply
with diagonal W, assumed inverted
    x[i] *= W_(i, i);

  *x_out = V_ * x;                                     
// premultiply with v.
}

//-----------------------------------------------------------------------------
//: Return N s.t. M * N = 0
template <class T>
vnl_matrix <T> vnl_svd<T>::nullspace()  const
{
  int k = rank();
  if (k == n_)
    vcl_cerr << "vnl_svd<T>::nullspace() -- Matrix is
full rank." << last_tol_ << vcl_endl;
  return nullspace(n_-k);
}

//-----------------------------------------------------------------------------
//: Return N s.t. M * N = 0
template <class T>
vnl_matrix <T> vnl_svd<T>::nullspace(int
required_nullspace_dimension)  const
{
  return V_.extract(V_.rows(),
required_nullspace_dimension, 0, n_ -
required_nullspace_dimension);
}

//-----------------------------------------------------------------------------
//: Return N s.t. M' * N = 0
template <class T>
vnl_matrix <T> vnl_svd<T>::left_nullspace()  const
{
  int k = rank();
  if (k == n_)
    vcl_cerr << "vnl_svd<T>::left_nullspace() -- Matrix
is full rank." << last_tol_ << vcl_endl;
  return U_.extract(U_.rows(), n_-k, 0, k);
}

//: Implementation to be done yet; currently returns
left_nullspace(). - PVR.
template <class T>
vnl_matrix<T> vnl_svd<T>::left_nullspace(int
/*required_nullspace_dimension*/) const
{
  return left_nullspace();
}


//-----------------------------------------------------------------------------
//: Return the rightmost column of V.
//  Does not check to see whether or not the matrix
actually was rank-deficient -
//  the caller is assumed to have examined W and
decided that to his or her satisfaction.
template <class T>
vnl_vector <T> vnl_svd<T>::nullvector()  const
{
  vnl_vector<T> ret(n_);
  for(int i = 0; i < n_; ++i)
    ret(i) = V_(i, n_-1);
  return ret;
}

//-----------------------------------------------------------------------------
//: Return the rightmost column of U.
//  Does not check to see whether or not the matrix
actually was rank-deficient.
template <class T>
vnl_vector <T> vnl_svd<T>::left_nullvector()  const
{
  vnl_vector<T> ret(m_);
  int col = vcl_min(m_, n_) - 1;
  for(int i = 0; i < m_; ++i)
    ret(i) = U_(i, col);
  return ret;
}

//--------------------------------------------------------------------------------

#undef VNL_SVD_INSTANTIATE
#define VNL_SVD_INSTANTIATE(T) \
template class vnl_svd<T >; \
template vcl_ostream& operator<<(vcl_ostream &,
vnl_svd<T > const &)

#endif // vnl_svd_txx_



----------------------------------------------------------------------

>Comment By: Amitha Perera (amithaperera)
Date: 2002-03-06 12:56

Message:
Logged In: YES 
user_id=237910

It would be more helpful for us if you just posted the
changes you made (instead of the whole thing), and describe
in more detail what you are trying to acheive with the changes.

Cheers,
Amitha.


----------------------------------------------------------------------

You can respond by visiting: 
http://sourceforge.net/tracker/?func=detail&atid=381096&aid=523333&group_id=24293