Document Tridiagonalization class, remove unused types.

1 files changed, 206 insertions, 44 deletions

diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index a41a5f670..024590f3c 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -30,14 +30,32 @@
   *
   * \class Tridiagonalization
   *
-  * \brief Trigiagonal decomposition of a selfadjoint matrix
+  * \brief Tridiagonal decomposition of a selfadjoint matrix
   *
-  * \param MatrixType the type of the matrix of which we are performing the tridiagonalization
+  * \tparam _MatrixType the type of the matrix of which we are computing the
+  * tridiagonal decomposition; this is expected to be an instantiation of the
+  * Matrix class template.
   *
   * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
   * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix.
   *
-  * \sa MatrixBase::tridiagonalize()
+  * A tridiagonal matrix is a matrix which has nonzero elements only on the
+  * main diagonal and the first diagonal below and above it. The Hessenberg
+  * decomposition of a selfadjoint matrix is in fact a tridiagonal
+  * decomposition. This class is used in SelfAdjointEigenSolver to compute the
+  * eigenvalues and eigenvectors of a selfadjoint matrix.
+  *
+  * Call the function compute() to compute the tridiagonal decomposition of a
+  * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&)
+  * constructor which computes the tridiagonal Schur decomposition at
+  * construction time. Once the decomposition is computed, you can use the
+  * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the
+  * decomposition.
+  *
+  * The documentation of Tridiagonalization(const MatrixType&) contains an
+  * example of the typical use of this class.
+  *
+  * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver
   */
 template<typename _MatrixType> class Tridiagonalization
 {
@@ -46,21 +64,18 @@ template<typename _MatrixType> class Tridiagonalization
     typedef _MatrixType MatrixType;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef typename ei_packet_traits<Scalar>::type Packet;
 
     enum {
       Size = MatrixType::RowsAtCompileTime,
       SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
       Options = MatrixType::Options,
       MaxSize = MatrixType::MaxRowsAtCompileTime,
-      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1,
-      PacketSize = ei_packet_traits<Scalar>::size
+      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
     };
 
     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
     typedef typename ei_plain_col_type<MatrixType, RealScalar>::type DiagonalType;
     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
-    typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
     
     typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
               typename Diagonal<MatrixType,0>::RealReturnType,
@@ -74,22 +89,53 @@ template<typename _MatrixType> class Tridiagonalization
                 Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >
             >::ret SubDiagonalReturnType;
 
-    /** This constructor initializes a Tridiagonalization object for
-      * further use with Tridiagonalization::compute()
+    /** \brief Default constructor.
+      *
+      * \param [in]  size  Positive integer, size of the matrix whose tridiagonal
+      * decomposition will be computed.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via compute().  The \p size parameter is only
+      * used as a hint. It is not an error to give a wrong \p size, but it may
+      * impair performance.
+      *
+      * \sa compute() for an example.
       */
     Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
       : m_matrix(size,size), m_hCoeffs(size-1)
     {}
 
+    /** \brief Constructor; computes tridiagonal decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
+      * is to be computed.
+      *
+      * This constructor calls compute() to compute the tridiagonal decomposition.
+      *
+      * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
+      * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
+      */
     Tridiagonalization(const MatrixType& matrix)
       : m_matrix(matrix), m_hCoeffs(matrix.cols()-1)
     {
       _compute(m_matrix, m_hCoeffs);
     }
 
-    /** Computes or re-compute the tridiagonalization for the matrix \a matrix.
+    /** \brief Computes tridiagonal decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
+      * is to be computed.
+      *
+      * The tridiagonal decomposition is computed by bringing the columns of
+      * the matrix successively in the required form using Householder
+      * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes
+      * the size of the given matrix.
+      *
+      * This method reuses of the allocated data in the Tridiagonalization
+      * object, if the size of the matrix does not change.
       *
-      * This method allows to re-use the allocated data.
+      * Example: \include Tridiagonalization_compute.cpp
+      * Output: \verbinclude Tridiagonalization_compute.out
       */
     void compute(const MatrixType& matrix)
     {
@@ -98,74 +144,191 @@ template<typename _MatrixType> class Tridiagonalization
       _compute(m_matrix, m_hCoeffs);
     }
 
-    /** \returns the householder coefficients allowing to
-      * reconstruct the matrix Q from the packed data.
+    /** \brief Returns the Householder coefficients.
+      *
+      * \returns a const reference to the vector of Householder coefficients
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
+      *
+      * The Householder coefficients allow the reconstruction of the matrix 
+      * \f$ Q \f$ in the tridiagonal decomposition from the packed data.
       *
-      * \sa packedMatrix()
+      * Example: \include Tridiagonalization_householderCoefficients.cpp
+      * Output: \verbinclude Tridiagonalization_householderCoefficients.out
+      *
+      * \sa packedMatrix(), \ref Householder_Module "Householder module"
       */
-    inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
+    inline CoeffVectorType householderCoefficients() const { return m_hCoeffs; }
 
-    /** \returns the internal result of the decomposition.
+    /** \brief Returns the internal representation of the decomposition 
+      *
+      *	\returns a const reference to a matrix with the internal representation
+      *	         of the decomposition.
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
       *
       * The returned matrix contains the following information:
-      *  - the strict upper part is equal to the input matrix A
-      *  - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real).
-      *  - the rest of the lower part contains the Householder vectors that, combined with
-      *    Householder coefficients returned by householderCoefficients(),
-      *    allows to reconstruct the matrix Q as follow:
-      *       Q = H_{N-1} ... H_1 H_0
-      *    where the matrices H are the Householder transformations:
-      *       H_i = (I - h_i * v_i * v_i')
-      *    where h_i == householderCoefficients()[i] and v_i is a Householder vector:
-      *       v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
+      *  - the strict upper triangular part is equal to the input matrix A.
+      *  - the diagonal and lower sub-diagonal represent the real tridiagonal
+      *    symmetric matrix T. 
+      *  - the rest of the lower part contains the Householder vectors that,
+      *    combined with Householder coefficients returned by
+      *    householderCoefficients(), allows to reconstruct the matrix Q as
+      *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
+      *    Here, the matrices \f$ H_i \f$ are the Householder transformations 
+      *       \f$ H_i = (I - h_i v_i v_i^T) \f$
+      *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and 
+      *    \f$ v_i \f$ is the Householder vector defined by
+      *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
+      *    with M the matrix returned by this function.
       *
       * See LAPACK for further details on this packed storage.
+      *
+      * Example: \include Tridiagonalization_packedMatrix.cpp
+      * Output: \verbinclude Tridiagonalization_packedMatrix.out
+      *
+      * \sa householderCoefficients()
       */
-    inline const MatrixType& packedMatrix(void) const { return m_matrix; }
+    inline const MatrixType& packedMatrix() const { return m_matrix; }
 
+    /** \brief Reconstructs the unitary matrix Q in the decomposition 
+      *
+      * \returns the matrix Q
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
+      *
+      * This function reconstructs the matrix Q from the Householder
+      * coefficients and the packed matrix stored internally. This
+      * reconstruction requires \f$ 4n^3 / 3 \f$ flops.
+      *
+      * \sa Tridiagonalization(const MatrixType&) for an example,
+      * matrixT(), matrixQInPlace()
+      */
     MatrixType matrixQ() const;
+
+    /** \brief Reconstructs the unitary matrix Q in the decomposition 
+      *
+      * This is an in-place variant of matrixQ() which avoids the copy. 
+      * This function will probably be deleted soon.
+      */
     template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
+
+    /** \brief Constructs the tridiagonal matrix T in the decomposition
+      *
+      * \returns the matrix T
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
+      *
+      * This function copies the matrix T from internal data. The diagonal and
+      * subdiagonal of the packed matrix as returned by packedMatrix()
+      * represents the matrix T. It may sometimes be sufficient to directly use
+      * the packed matrix or the vector expressions returned by diagonal()
+      * and subDiagonal() instead of creating a new matrix with this function.
+      *
+      * \sa Tridiagonalization(const MatrixType&) for an example,
+      * matrixQ(), packedMatrix(), diagonal(), subDiagonal()
+      */
     MatrixType matrixT() const;
-    const DiagonalReturnType diagonal(void) const;
-    const SubDiagonalReturnType subDiagonal(void) const;
 
-    static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
+    /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition.
+      *
+      * \returns expression representing the diagonal of T
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
+      *
+      * Example: \include Tridiagonalization_diagonal.cpp
+      * Output: \verbinclude Tridiagonalization_diagonal.out
+      *
+      * \sa matrixT(), subDiagonal()
+      */
+    const DiagonalReturnType diagonal() const;
 
-    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+    /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
+      *
+      * \returns expression representing the subdiagonal of T
+      *
+      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called before
+      * to compute the tridiagonal decomposition of a matrix.
+      *
+      * \sa diagonal() for an example, matrixT()
+      */
+    const SubDiagonalReturnType subDiagonal() const;
+
+    /** \brief Performs a full decomposition in place 
+      *
+      * \param[in,out]  mat  On input, the selfadjoint matrix whose tridiagonal
+      *    decomposition is to be computed. On output, the orthogonal matrix Q
+      *    in the decomposition if \p extractQ is true.
+      * \param[out]  diag  The diagonal of the tridiagonal matrix T in the
+      *    decomposition.
+      * \param[out]  subdiag  The subdiagonal of the tridiagonal matrix T in
+      *    the decomposition. 
+      * \param[in]  extractQ  If true, the orthogonal matrix Q in the
+      *    decomposition is computed and stored in \p mat.
+      *
+      * Compute the tridiagonal matrix of \p mat in place. The tridiagonal
+      * matrix T is passed to the output parameters \p diag and \p subdiag. If
+      * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat.
+      *
+      * The vectors \p diag and \p subdiag are not resized. The function
+      * assumes that they are already of the correct size. The length of the
+      * vector \p diag should equal the number of rows in \p mat, and the
+      * length of the vector \p subdiag should be one left.
+      *
+      * This implementation contains an optimized path for real 3-by-3 matrices
+      * which is especially useful for plane fitting.
+      *
+      * \note Notwithstanding the name, the current implementation copies
+      * \p mat to a temporary matrix and uses that matrix to compute the
+      * decomposition. 
+      *
+      * Example (this uses the same matrix as the example in
+      *    Tridiagonalization(const MatrixType&)): 
+      *    \include Tridiagonalization_decomposeInPlace.cpp
+      * Output: \verbinclude Tridiagonalization_decomposeInPlace.out
+      *
+      * \sa Tridiagonalization(const MatrixType&), compute()
+      */
+    static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
 
   protected:
 
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
     static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
 
     MatrixType m_matrix;
     CoeffVectorType m_hCoeffs;
 };
 
-/** \returns an expression of the diagonal vector */
 template<typename MatrixType>
 const typename Tridiagonalization<MatrixType>::DiagonalReturnType
-Tridiagonalization<MatrixType>::diagonal(void) const
+Tridiagonalization<MatrixType>::diagonal() const
 {
   return m_matrix.diagonal();
 }
 
-/** \returns an expression of the sub-diagonal vector */
 template<typename MatrixType>
 const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
-Tridiagonalization<MatrixType>::subDiagonal(void) const
+Tridiagonalization<MatrixType>::subDiagonal() const
 {
   int n = m_matrix.rows();
   return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
 }
 
-/** constructs and returns the tridiagonal matrix T.
-  * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
-  * Therefore, it might be often sufficient to directly use the packed matrix, or the vector
-  * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
-  */
 template<typename MatrixType>
 typename Tridiagonalization<MatrixType>::MatrixType
-Tridiagonalization<MatrixType>::matrixT(void) const
+Tridiagonalization<MatrixType>::matrixT() const
 {
   // FIXME should this function (and other similar ones) rather take a matrix as argument
   // and fill it ? (to avoid temporaries)
@@ -223,10 +386,9 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
   }
 }
 
-/** reconstructs and returns the matrix Q */
 template<typename MatrixType>
 typename Tridiagonalization<MatrixType>::MatrixType
-Tridiagonalization<MatrixType>::matrixQ(void) const
+Tridiagonalization<MatrixType>::matrixQ() const
 {
   MatrixType matQ;
   matrixQInPlace(&matQ);
@@ -240,6 +402,7 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
   QDerived& matQ = q->derived();
   int n = m_matrix.rows();
   matQ = MatrixType::Identity(n,n);
+  typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
   RowVectorType aux(n);
   for (int i = n-2; i>=0; i--)
   {
@@ -248,7 +411,6 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
   }
 }
 
-/** Performs a full decomposition in place */
 template<typename MatrixType>
 void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
 {