Attachment #964 for bug #1763

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

View | Details | Raw Unified | Return to bug 1763
Collapse All | Expand All




  : public traits<typename MatrixType::PlainObject>
{
  typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
  enum { Flags = 0 };
};

template<typename MatrixType, typename CoeffVectorType>
EIGEN_DEVICE_FUNC
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, typename MatrixType::Scalar* workspace);

EIGEN_DEVICE_FUNC
Index tridiagonalization_inplace_blocked_workspace_size(Index N, Index r);

}

/** \eigenvalues_module \ingroup Eigenvalues_Module
  *
  *
  * \class Tridiagonalization
  *
  * \brief Tridiagonal decomposition of a selfadjoint matrix

    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3

    enum {
      Size = MatrixType::RowsAtCompileTime,
      SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
      Options = MatrixType::Options,
      MaxSize = MatrixType::MaxRowsAtCompileTime,
      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1),
      BlockThreshold = 300,
      BlockSize = 60
    };

    typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
    typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
    typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
    typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
    typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
    typedef Matrix<Scalar, Dynamic, 1, Options & ~RowMajor> WorkspaceType;

    typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
              typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
              const Diagonal<const MatrixType>
            >::type DiagonalReturnType;

    typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
              typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,

      *
      * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
      * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
      */
    template<typename InputType>
    explicit Tridiagonalization(const EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()),
        m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
        m_workspace(matrix.cols() > BlockThreshold ? internal::tridiagonalization_inplace_blocked_workspace_size(matrix.cols(), BlockSize) : 0),
        m_isInitialized(false)
    {
      internal::tridiagonalization_inplace(m_matrix, m_hCoeffs, m_workspace.data());
      m_isInitialized = true;
    }

    /** \brief Computes tridiagonal decomposition of given matrix.
      *
      * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
      * is to be computed.
      * \returns    Reference to \c *this

      * Example: \include Tridiagonalization_compute.cpp
      * Output: \verbinclude Tridiagonalization_compute.out
      */
    template<typename InputType>
    Tridiagonalization& compute(const EigenBase<InputType>& matrix)
    {
      m_matrix = matrix.derived();
      m_hCoeffs.resize(matrix.rows()-1, 1);
      if(matrix.cols() > BlockThreshold){
        m_workspace.resize(internal::tridiagonalization_inplace_blocked_workspace_size(matrix.cols(), BlockSize));
      }
      internal::tridiagonalization_inplace(m_matrix, m_hCoeffs, m_workspace.data());
      m_isInitialized = true;
      return *this;
    }

    /** \brief Returns the Householder coefficients.
      *
      * \returns a const reference to the vector of Householder coefficients
      *

      * \sa diagonal() for an example, matrixT()
      */
    SubDiagonalReturnType subDiagonal() const;

  protected:

    MatrixType m_matrix;
    CoeffVectorType m_hCoeffs;
    WorkspaceType m_workspace;
    bool m_isInitialized;
};

template<typename MatrixType>
typename Tridiagonalization<MatrixType>::DiagonalReturnType
Tridiagonalization<MatrixType>::diagonal() const
{
  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");

typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
Tridiagonalization<MatrixType>::subDiagonal() const
{
  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
  return m_matrix.template diagonal<-1>().real();
}

namespace internal {

/** \internal
  * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
  *
  * \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.
  *                     On output, the strict upper part is left unchanged, and the lower triangular part
  *                     represents the T and Q matrices in packed format has detailed below.
  * \param[out]    hCoeffs returned Householder coefficients (see below)
  *

  *       \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.
  *
  * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
  *
  * \sa Tridiagonalization::packedMatrix()
  */
template<typename MatrixType, typename CoeffVectorType>
EIGEN_DEVICE_FUNC
void tridiagonalization_inplace_unblocked(MatrixType& matA, CoeffVectorType& hCoeffs)
{
  using numext::conj;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  Index n = matA.rows();
  eigen_assert(n==matA.cols());
  eigen_assert(n==hCoeffs.size()+1 || n==1);
  

    matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
      .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));

    matA.col(i).coeffRef(i+1) = beta;
    hCoeffs.coeffRef(i) = h;
  }
}

/**
 * Calculates required workspace size for block householder tridiagonalization
 * of a matrix of size n * n, by using block size of r.
 * @param N matrix dimenstion
 * @param r block size
 * @return minimum workspace size
 */
EIGEN_DEVICE_FUNC
Index tridiagonalization_inplace_blocked_workspace_size(Index N, Index r){
  Index first_actual_r = std::min({r, static_cast<Index>(N - 1)});
  Index V_size = (N - 1) * r;
  Index partial_update_size = (N - first_actual_r) * (N - first_actual_r);
  return V_size + partial_update_size;
}

/**
 * Tridiagonalize a selfadjoint matrix using block Housholder algorithm.
 * A = Q * T * Q^T, where T is tridiagonal and Q is unitary.
 * @param[in,out] packed On input the input matrix. On output packed form of the
 * tridiagonal matrix. Elements of the resulting selfadjoint tridiagonal matrix
 * T are in the diagonal and first subdiagonal, which contains subdiagonal
 * entries of T. Columns bellow diagonal contain householder vectors that can be
 * used to construct unitary matrix Q.
 * @param[out] hCoeffs householder coefficients
 * @param workspace Workspace. Minimal required size can be determined by
 * calling `tridiagonalization_inplace_blocked_workspace_size()`
 * @param r Block size. Affects only performance of the algorithm. Optimal value
 * depends on the size of A and cache of the processor. For larger matrices or
 * larger cache sizes a larger value is optimal.
 */
template<typename MatrixType, typename CoeffVectorType>
EIGEN_DEVICE_FUNC
void tridiagonalization_inplace_blocked(
    MatrixType& packed,
    CoeffVectorType& hCoeffs,
    typename MatrixType::Scalar* workspace, const Index r) {
  using Scalar = typename MatrixType::Scalar;
  using Real = typename NumTraits<Scalar>::Real;
  using MapType = typename Matrix<Scalar, Dynamic, Dynamic>::MapType;
  using MapTypeVec = typename Matrix<Scalar, Dynamic, 1>::MapType;
  Index N = packed.rows();
  for (Index k = 0; k < N - 1; k += r) {
    const Index actual_r = std::min({r, static_cast<Index>(N - k - 1)});
    MapType V(workspace, N - k - 1, actual_r);
    for (Index j = 0; j < actual_r; j++) {
      if (j != 0) {
        Index householder_whole_size = N - k - j;
        Scalar temp = packed(k + j, k + j - 1);
        packed(k + j, k + j - 1) = Scalar(1);
        packed.col(k + j).tail(householder_whole_size) -= packed.block(k + j, k, householder_whole_size, j) * V.block(j - 1, 0, 1, j).adjoint();
        packed.col(k + j).tail(householder_whole_size) -= V.block(j - 1, 0, householder_whole_size, j) * packed.block(k + j, k, 1, j).adjoint();
        packed(k + j, k + j - 1) = temp;
      }
      typename Matrix<Scalar, Dynamic, Dynamic>::ColXpr::SegmentReturnType householder = packed.col(k + j).tail(N- k - j - 1);
      Scalar tau;
      Real beta;
      householder.makeHouseholderInPlace(tau, beta);
      householder[0] = Scalar(1);

      typename MapType::ColXpr::SegmentReturnType v = V.col(j).tail(householder.size());
      v.noalias() = packed.bottomRightCorner(N - k - j - 1,N - k - j - 1)
                        .template selfadjointView<Lower>() * householder * numext::conj(tau);
      MapTypeVec tmp(workspace + V.size(), j);

      //Reminder of packed is not transformed by current block yet - v needs some fixes
      tmp.noalias() = V.bottomLeftCorner(householder.size(), j).adjoint() * householder * numext::conj(tau);
      v.noalias() -= packed.block(k + j + 1, k, householder.size(), j) * tmp;
      tmp.noalias() = packed.block(k + j + 1, k, householder.size(), j).adjoint() * householder * numext::conj(tau);
      v.noalias() -= V.bottomLeftCorner(householder.size(), j) * tmp;

      const Scalar cnst = (v.adjoint() * householder)[0];
      v -= Real(0.5)  * numext::conj(tau) * cnst * householder;

      //store householder transformation scaling and subdiagonal of T
      hCoeffs(k + j) = tau;
      packed(k + j + 1, k + j) = beta;
    }
    //update reminder of packed with the last block
    MapType partial_update(workspace + V.size(), V.rows() - actual_r + 1,V.rows() - actual_r + 1);
    Scalar tmp = packed(k + actual_r, k+actual_r-1);
    packed(k + actual_r, k+actual_r-1) = Scalar(1);
    partial_update.noalias() = packed.block(k + actual_r, k, N - k - actual_r, actual_r) * V.bottomRows(V.rows() - actual_r + 1).adjoint();
    packed(k + actual_r, k+actual_r-1) = tmp;
    packed.block(k + actual_r, k + actual_r, N - k - actual_r, N - k - actual_r).template triangularView<Lower>()
        -= partial_update + partial_update.adjoint();
  }
}

template<typename MatrixType, typename CoeffVectorType>
EIGEN_DEVICE_FUNC
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, typename MatrixType::Scalar* workspace)
{
  if(matA.rows() > Tridiagonalization<MatrixType>::BlockThreshold)
    tridiagonalization_inplace_blocked(matA, hCoeffs, workspace, Tridiagonalization<MatrixType>::BlockSize);
  else {
    tridiagonalization_inplace_unblocked(matA, hCoeffs);
  }
}

// forward declaration, implementation at the end of this file
template<typename MatrixType,
         int Size=MatrixType::ColsAtCompileTime,
         bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
struct tridiagonalization_inplace_selector;

/** \brief Performs a full tridiagonalization in place
  *

/** \internal
  * General full tridiagonalization
  */
template<typename MatrixType, int Size, bool IsComplex>
struct tridiagonalization_inplace_selector
{
  typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
  typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
  typedef typename Tridiagonalization<MatrixType>::WorkspaceType WorkspaceType;
  template<typename DiagonalType, typename SubDiagonalType>
  static EIGEN_DEVICE_FUNC
  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
  {
    CoeffVectorType hCoeffs(mat.cols()-1);
    if(mat.cols()>Tridiagonalization<MatrixType>::BlockThreshold){
      WorkspaceType workspace(tridiagonalization_inplace_blocked_workspace_size(mat.cols(), Tridiagonalization<MatrixType>::BlockSize));
      tridiagonalization_inplace_blocked(mat, hCoeffs, workspace.data());
    }
    else {
      tridiagonalization_inplace_unblocked(mat, hCoeffs);
    }
    diag = mat.diagonal().real();
    subdiag = mat.template diagonal<-1>().real();
    if(extractQ)
      mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
            .setLength(mat.rows() - 1)
            .setShift(1);
  }
};

Return to bug 1763

Lines 21-37 struct traits<TridiagonalizationMatrixTR Link Here

(-)a/Eigen/src/Eigenvalues/Tridiagonalization.h (-7 / +124 lines)
21	: public traits<typename MatrixType::PlainObject>	21	: public traits<typename MatrixType::PlainObject>
22	{	22	{
23	typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?	23	typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
24	enum { Flags = 0 };	24	enum { Flags = 0 };
25	};	25	};
26		26
27	template<typename MatrixType, typename CoeffVectorType>	27	template<typename MatrixType, typename CoeffVectorType>
28	EIGEN_DEVICE_FUNC	28	EIGEN_DEVICE_FUNC
29	void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);	29	void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, typename MatrixType::Scalar* workspace);
		30
		31	EIGEN_DEVICE_FUNC
		32	Index tridiagonalization_inplace_blocked_workspace_size(Index N, Index r);
		33
30	}	34	}
31		35
32	/** \eigenvalues_module \ingroup Eigenvalues_Module	36	/** \eigenvalues_module \ingroup Eigenvalues_Module
33	*	37	*
34	*	38	*
35	* \class Tridiagonalization	39	* \class Tridiagonalization
36	*	40	*
37	* \brief Tridiagonal decomposition of a selfadjoint matrix	41	* \brief Tridiagonal decomposition of a selfadjoint matrix
Lines 72-95 template<typename _MatrixType> class Tri Link Here
72	typedef typename NumTraits<Scalar>::Real RealScalar;	76	typedef typename NumTraits<Scalar>::Real RealScalar;
73	typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3	77	typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
74		78
75	enum {	79	enum {
76	Size = MatrixType::RowsAtCompileTime,	80	Size = MatrixType::RowsAtCompileTime,
77	SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),	81	SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
78	Options = MatrixType::Options,	82	Options = MatrixType::Options,
79	MaxSize = MatrixType::MaxRowsAtCompileTime,	83	MaxSize = MatrixType::MaxRowsAtCompileTime,
80	MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)	84	MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1),
		85	BlockThreshold = 300,
		86	BlockSize = 60
81	};	87	};
82		88
83	typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;	89	typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
84	typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;	90	typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
85	typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;	91	typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
86	typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;	92	typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
87	typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;	93	typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
		94	typedef Matrix<Scalar, Dynamic, 1, Options & ~RowMajor> WorkspaceType;
88		95
89	typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,	96	typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
90	typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,	97	typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
91	const Diagonal<const MatrixType>	98	const Diagonal<const MatrixType>
92	>::type DiagonalReturnType;	99	>::type DiagonalReturnType;
93		100
94	typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,	101	typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
95	typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,	102	typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
Lines 126-144 template<typename _MatrixType> class Tri Link Here
126	*	133	*
127	* Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp	134	* Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
128	* Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out	135	* Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
129	*/	136	*/
130	template<typename InputType>	137	template<typename InputType>
131	explicit Tridiagonalization(const EigenBase<InputType>& matrix)	138	explicit Tridiagonalization(const EigenBase<InputType>& matrix)
132	: m_matrix(matrix.derived()),	139	: m_matrix(matrix.derived()),
133	m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),	140	m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
		141	m_workspace(matrix.cols() > BlockThreshold ? internal::tridiagonalization_inplace_blocked_workspace_size(matrix.cols(), BlockSize) : 0),
134	m_isInitialized(false)	142	m_isInitialized(false)
135	{	143	{
136	internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);	144	internal::tridiagonalization_inplace(m_matrix, m_hCoeffs, m_workspace.data());
137	m_isInitialized = true;	145	m_isInitialized = true;
138	}	146	}
139		147
140	/** \brief Computes tridiagonal decomposition of given matrix.	148	/** \brief Computes tridiagonal decomposition of given matrix.
141	*	149	*
142	* \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition	150	* \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition
143	* is to be computed.	151	* is to be computed.
144	* \returns Reference to \c *this	152	* \returns Reference to \c *this
Lines 154-170 template<typename _MatrixType> class Tri Link Here
154	* Example: \include Tridiagonalization_compute.cpp	162	* Example: \include Tridiagonalization_compute.cpp
155	* Output: \verbinclude Tridiagonalization_compute.out	163	* Output: \verbinclude Tridiagonalization_compute.out
156	*/	164	*/
157	template<typename InputType>	165	template<typename InputType>
158	Tridiagonalization& compute(const EigenBase<InputType>& matrix)	166	Tridiagonalization& compute(const EigenBase<InputType>& matrix)
159	{	167	{
160	m_matrix = matrix.derived();	168	m_matrix = matrix.derived();
161	m_hCoeffs.resize(matrix.rows()-1, 1);	169	m_hCoeffs.resize(matrix.rows()-1, 1);
162	internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);	170	if(matrix.cols() > BlockThreshold){
		171	m_workspace.resize(internal::tridiagonalization_inplace_blocked_workspace_size(matrix.cols(), BlockSize));
		172	}
		173	internal::tridiagonalization_inplace(m_matrix, m_hCoeffs, m_workspace.data());
163	m_isInitialized = true;	174	m_isInitialized = true;
164	return *this;	175	return *this;
165	}	176	}
166		177
167	/** \brief Returns the Householder coefficients.	178	/** \brief Returns the Householder coefficients.
168	*	179	*
169	* \returns a const reference to the vector of Householder coefficients	180	* \returns a const reference to the vector of Householder coefficients
170	*	181	*
Lines 295-310 template<typename _MatrixType> class Tri Link Here
295	* \sa diagonal() for an example, matrixT()	306	* \sa diagonal() for an example, matrixT()
296	*/	307	*/
297	SubDiagonalReturnType subDiagonal() const;	308	SubDiagonalReturnType subDiagonal() const;
298		309
299	protected:	310	protected:
300		311
301	MatrixType m_matrix;	312	MatrixType m_matrix;
302	CoeffVectorType m_hCoeffs;	313	CoeffVectorType m_hCoeffs;
		314	WorkspaceType m_workspace;
303	bool m_isInitialized;	315	bool m_isInitialized;
304	};	316	};
305		317
306	template<typename MatrixType>	318	template<typename MatrixType>
307	typename Tridiagonalization<MatrixType>::DiagonalReturnType	319	typename Tridiagonalization<MatrixType>::DiagonalReturnType
308	Tridiagonalization<MatrixType>::diagonal() const	320	Tridiagonalization<MatrixType>::diagonal() const
309	{	321	{
310	eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");	322	eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
Lines 315-331 template<typename MatrixType> Link Here
315	typename Tridiagonalization<MatrixType>::SubDiagonalReturnType	327	typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
316	Tridiagonalization<MatrixType>::subDiagonal() const	328	Tridiagonalization<MatrixType>::subDiagonal() const
317	{	329	{
318	eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");	330	eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
319	return m_matrix.template diagonal<-1>().real();	331	return m_matrix.template diagonal<-1>().real();
320	}	332	}
321		333
322	namespace internal {	334	namespace internal {
323
324	/** \internal	335	/** \internal
325	* Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.	336	* Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
326	*	337	*
327	* \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.	338	* \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.
328	* On output, the strict upper part is left unchanged, and the lower triangular part	339	* On output, the strict upper part is left unchanged, and the lower triangular part
329	* represents the T and Q matrices in packed format has detailed below.	340	* represents the T and Q matrices in packed format has detailed below.
330	* \param[out] hCoeffs returned Householder coefficients (see below)	341	* \param[out] hCoeffs returned Householder coefficients (see below)
331	*	342	*
Lines 341-357 namespace internal { Link Here
341	* \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.	352	* \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.
342	*	353	*
343	* Implemented from Golub's "Matrix Computations", algorithm 8.3.1.	354	* Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
344	*	355	*
345	* \sa Tridiagonalization::packedMatrix()	356	* \sa Tridiagonalization::packedMatrix()
346	*/	357	*/
347	template<typename MatrixType, typename CoeffVectorType>	358	template<typename MatrixType, typename CoeffVectorType>
348	EIGEN_DEVICE_FUNC	359	EIGEN_DEVICE_FUNC
349	void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)	360	void tridiagonalization_inplace_unblocked(MatrixType& matA, CoeffVectorType& hCoeffs)
350	{	361	{
351	using numext::conj;	362	using numext::conj;
352	typedef typename MatrixType::Scalar Scalar;	363	typedef typename MatrixType::Scalar Scalar;
353	typedef typename MatrixType::RealScalar RealScalar;	364	typedef typename MatrixType::RealScalar RealScalar;
354	Index n = matA.rows();	365	Index n = matA.rows();
355	eigen_assert(n==matA.cols());	366	eigen_assert(n==matA.cols());
356	eigen_assert(n==hCoeffs.size()+1 \|\| n==1);	367	eigen_assert(n==hCoeffs.size()+1 \|\| n==1);
357		368
Lines 374-389 void tridiagonalization_inplace(MatrixTy Link Here
374	matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()	385	matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
375	.rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));	386	.rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
376		387
377	matA.col(i).coeffRef(i+1) = beta;	388	matA.col(i).coeffRef(i+1) = beta;
378	hCoeffs.coeffRef(i) = h;	389	hCoeffs.coeffRef(i) = h;
379	}	390	}
380	}	391	}
381		392
		393	/**
		394	* Calculates required workspace size for block householder tridiagonalization
		395	* of a matrix of size n * n, by using block size of r.
		396	* @param N matrix dimenstion
		397	* @param r block size
		398	* @return minimum workspace size
		399	*/
		400	EIGEN_DEVICE_FUNC
		401	Index tridiagonalization_inplace_blocked_workspace_size(Index N, Index r){
		402	Index first_actual_r = std::min({r, static_cast<Index>(N - 1)});
		403	Index V_size = (N - 1) * r;
		404	Index partial_update_size = (N - first_actual_r) * (N - first_actual_r);
		405	return V_size + partial_update_size;
		406	}
		407
		408	/**
		409	* Tridiagonalize a selfadjoint matrix using block Housholder algorithm.
		410	* A = Q * T * Q^T, where T is tridiagonal and Q is unitary.
		411	* @param[in,out] packed On input the input matrix. On output packed form of the
		412	* tridiagonal matrix. Elements of the resulting selfadjoint tridiagonal matrix
		413	* T are in the diagonal and first subdiagonal, which contains subdiagonal
		414	* entries of T. Columns bellow diagonal contain householder vectors that can be
		415	* used to construct unitary matrix Q.
		416	* @param[out] hCoeffs householder coefficients
		417	* @param workspace Workspace. Minimal required size can be determined by
		418	* calling `tridiagonalization_inplace_blocked_workspace_size()`
		419	* @param r Block size. Affects only performance of the algorithm. Optimal value
		420	* depends on the size of A and cache of the processor. For larger matrices or
		421	* larger cache sizes a larger value is optimal.
		422	*/
		423	template<typename MatrixType, typename CoeffVectorType>
		424	EIGEN_DEVICE_FUNC
		425	void tridiagonalization_inplace_blocked(
		426	MatrixType& packed,
		427	CoeffVectorType& hCoeffs,
		428	typename MatrixType::Scalar* workspace, const Index r) {
		429	using Scalar = typename MatrixType::Scalar;
		430	using Real = typename NumTraits<Scalar>::Real;
		431	using MapType = typename Matrix<Scalar, Dynamic, Dynamic>::MapType;
		432	using MapTypeVec = typename Matrix<Scalar, Dynamic, 1>::MapType;
		433	Index N = packed.rows();
		434	for (Index k = 0; k < N - 1; k += r) {
		435	const Index actual_r = std::min({r, static_cast<Index>(N - k - 1)});
		436	MapType V(workspace, N - k - 1, actual_r);
		437	for (Index j = 0; j < actual_r; j++) {
		438	if (j != 0) {
		439	Index householder_whole_size = N - k - j;
		440	Scalar temp = packed(k + j, k + j - 1);
		441	packed(k + j, k + j - 1) = Scalar(1);
		442	packed.col(k + j).tail(householder_whole_size) -= packed.block(k + j, k, householder_whole_size, j) * V.block(j - 1, 0, 1, j).adjoint();
		443	packed.col(k + j).tail(householder_whole_size) -= V.block(j - 1, 0, householder_whole_size, j) * packed.block(k + j, k, 1, j).adjoint();
		444	packed(k + j, k + j - 1) = temp;
		445	}
		446	typename Matrix<Scalar, Dynamic, Dynamic>::ColXpr::SegmentReturnType householder = packed.col(k + j).tail(N- k - j - 1);
		447	Scalar tau;
		448	Real beta;
		449	householder.makeHouseholderInPlace(tau, beta);
		450	householder[0] = Scalar(1);
		451
		452	typename MapType::ColXpr::SegmentReturnType v = V.col(j).tail(householder.size());
		453	v.noalias() = packed.bottomRightCorner(N - k - j - 1,N - k - j - 1)
		454	.template selfadjointView<Lower>() * householder * numext::conj(tau);
		455	MapTypeVec tmp(workspace + V.size(), j);
		456
457	//Reminder of packed is not transformed by current block yet - v needs some fixes
458	tmp.noalias() = V.bottomLeftCorner(householder.size(), j).adjoint() * householder * numext::conj(tau);
459	v.noalias() -= packed.block(k + j + 1, k, householder.size(), j) * tmp;
460	tmp.noalias() = packed.block(k + j + 1, k, householder.size(), j).adjoint() * householder * numext::conj(tau);
461	v.noalias() -= V.bottomLeftCorner(householder.size(), j) * tmp;
462
463	const Scalar cnst = (v.adjoint() * householder)[0];
464	v -= Real(0.5) * numext::conj(tau) * cnst * householder;
465
466	//store householder transformation scaling and subdiagonal of T
467	hCoeffs(k + j) = tau;
468	packed(k + j + 1, k + j) = beta;
469	}
470	//update reminder of packed with the last block
471	MapType partial_update(workspace + V.size(), V.rows() - actual_r + 1,V.rows() - actual_r + 1);
472	Scalar tmp = packed(k + actual_r, k+actual_r-1);
473	packed(k + actual_r, k+actual_r-1) = Scalar(1);
474	partial_update.noalias() = packed.block(k + actual_r, k, N - k - actual_r, actual_r) * V.bottomRows(V.rows() - actual_r + 1).adjoint();
475	packed(k + actual_r, k+actual_r-1) = tmp;
476	packed.block(k + actual_r, k + actual_r, N - k - actual_r, N - k - actual_r).template triangularView<Lower>()
477	-= partial_update + partial_update.adjoint();
478	}
479	}
480
481	template<typename MatrixType, typename CoeffVectorType>
482	EIGEN_DEVICE_FUNC
483	void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, typename MatrixType::Scalar* workspace)
484	{
485	if(matA.rows() > Tridiagonalization<MatrixType>::BlockThreshold)
486	tridiagonalization_inplace_blocked(matA, hCoeffs, workspace, Tridiagonalization<MatrixType>::BlockSize);
487	else {
488	tridiagonalization_inplace_unblocked(matA, hCoeffs);
489	}
490	}
491
382	// forward declaration, implementation at the end of this file	492	// forward declaration, implementation at the end of this file
383	template<typename MatrixType,	493	template<typename MatrixType,
384	int Size=MatrixType::ColsAtCompileTime,	494	int Size=MatrixType::ColsAtCompileTime,
385	bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>	495	bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
386	struct tridiagonalization_inplace_selector;	496	struct tridiagonalization_inplace_selector;
387		497
388	/** \brief Performs a full tridiagonalization in place	498	/** \brief Performs a full tridiagonalization in place
389	*	499	*
Lines 436-457 void tridiagonalization_inplace(MatrixTy Link Here
436	/** \internal	546	/** \internal
437	* General full tridiagonalization	547	* General full tridiagonalization
438	*/	548	*/
439	template<typename MatrixType, int Size, bool IsComplex>	549	template<typename MatrixType, int Size, bool IsComplex>
440	struct tridiagonalization_inplace_selector	550	struct tridiagonalization_inplace_selector
441	{	551	{
442	typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;	552	typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
443	typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;	553	typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
		554	typedef typename Tridiagonalization<MatrixType>::WorkspaceType WorkspaceType;
444	template<typename DiagonalType, typename SubDiagonalType>	555	template<typename DiagonalType, typename SubDiagonalType>
445	static EIGEN_DEVICE_FUNC	556	static EIGEN_DEVICE_FUNC
446	void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)	557	void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
447	{	558	{
448	CoeffVectorType hCoeffs(mat.cols()-1);	559	CoeffVectorType hCoeffs(mat.cols()-1);
449	tridiagonalization_inplace(mat,hCoeffs);	560	if(mat.cols()>Tridiagonalization<MatrixType>::BlockThreshold){
		561	WorkspaceType workspace(tridiagonalization_inplace_blocked_workspace_size(mat.cols(), Tridiagonalization<MatrixType>::BlockSize));
		562	tridiagonalization_inplace_blocked(mat, hCoeffs, workspace.data());
		563	}
		564	else {
		565	tridiagonalization_inplace_unblocked(mat, hCoeffs);
		566	}
450	diag = mat.diagonal().real();	567	diag = mat.diagonal().real();
451	subdiag = mat.template diagonal<-1>().real();	568	subdiag = mat.template diagonal<-1>().real();
452	if(extractQ)	569	if(extractQ)
453	mat = HouseholderSequenceType(mat, hCoeffs.conjugate())	570	mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
454	.setLength(mat.rows() - 1)	571	.setLength(mat.rows() - 1)
455	.setShift(1);	572	.setShift(1);
456	}	573	}
457	};	574	};