vnl_qr.h

InsightToolkit-1.4.0(有大量的优化算法程序)

// This is core/vnl/algo/vnl_qr.h
#ifndef vnl_qr_h_
#define vnl_qr_h_
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma interface
#endif
//:
//  \file
//  \brief Calculate inverse of a matrix using QR
//  \author  Andrew W. Fitzgibbon, Oxford RRG
//  \date   08 Dec 96
//
// \verbatim
//  Modifications:
//  081296 AWF Temporarily abandoned as I realized my problem was symmetric...
//  080697 AWF Recovered, implemented solve().
//  200897 AWF Added determinant().
//  071097 AWF Added Q(), R().
//  Christian Stoecklin, ETH Zurich, added QtB(v)
//  31-mar-2000 fsm: templated
//  dac (Manchester) 28/03/2001: tidied up documentation
//  13 Jan.2003 - Peter Vanroose - added missing implementation for inverse(),
//                                 tinverse(), solve(matrix), extract_q_and_r().
// \endverbatim

#include <vnl/vnl_vector.h>
#include <vnl/vnl_matrix.h>

//: Extract the Q*R decomposition of matrix M.
//  The decomposition is stored in a compact and time-efficient
// packed form, which is most easily used via the "solve" and
// "determinant" methods.

export template <class T>
class vnl_qr
{
 public:
  vnl_qr(vnl_matrix<T> const & M);
 ~vnl_qr();

  //: return the inverse matrix of M
  vnl_matrix<T> inverse () const;
  //: return the transpose of the inverse matrix of M
  vnl_matrix<T> tinverse () const;
  //: return the original matrix M
  vnl_matrix<T> recompose () const;

  //: Solve equation M x = rhs for x using the computed decomposition.
  vnl_matrix<T> solve (const vnl_matrix<T>& rhs) const;
  //: Solve equation M x = rhs for x using the computed decomposition.
  vnl_vector<T> solve (const vnl_vector<T>& rhs) const;

  //: Return the determinant of M.  This is computed from M = Q R as follows:
  // |M| = |Q| |R|.
  // |R| is the product of the diagonal elements.
  // |Q| is (-1)^n as it is a product of Householder reflections.
  // So det = -prod(-r_ii).
  T determinant() const;
  //: Unpack and return unitary part Q.
  vnl_matrix<T> const& Q() const;
  //: Unpack and return R.
  vnl_matrix<T> const& R() const;
  //: Return residual vector d of M x = b -> d = Q'b
  vnl_vector<T> QtB(const vnl_vector<T>& b) const;

  void extract_q_and_r(vnl_matrix<T>* q, vnl_matrix<T>* r) const { *q = Q(); *r = R(); }

 private:
  vnl_matrix<T> qrdc_out_;
  vnl_vector<T> qraux_;
  vnl_vector<int> jpvt_;
  vnl_matrix<T>* Q_;
  vnl_matrix<T>* R_;

  // Disallow assignment.
  vnl_qr(const vnl_qr<T> &) { }
  void operator=(const vnl_qr<T> &) { }
};

//: Compute determinant of matrix "M" using QR.
template <class T>
inline T vnl_qr_determinant(vnl_matrix<T> const& m)
{
  return vnl_qr<T>(m).determinant();
}

export template <class T>
vcl_ostream& operator<<(vcl_ostream&, vnl_qr<T> const & qr);

#endif // vnl_qr_h_