10 #ifndef EIGEN_MATRIX_FUNCTION
11 #define EIGEN_MATRIX_FUNCTION
13 #include "StemFunction.h"
14 #include "MatrixFunctionAtomic.h"
34 template <
typename MatrixType,
36 int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
59 template <
typename ResultType>
60 void compute(ResultType &result);
67 template <
typename MatrixType,
typename AtomicType>
72 typedef internal::traits<MatrixType> Traits;
73 typedef typename Traits::Scalar Scalar;
74 static const int Rows = Traits::RowsAtCompileTime;
75 static const int Cols = Traits::ColsAtCompileTime;
76 static const int Options = MatrixType::Options;
77 static const int MaxRows = Traits::MaxRowsAtCompileTime;
78 static const int MaxCols = Traits::MaxColsAtCompileTime;
80 typedef std::complex<Scalar> ComplexScalar;
90 MatrixFunction(
const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
101 template <
typename ResultType>
102 void compute(ResultType& result)
104 ComplexMatrix CA = m_A.template cast<ComplexScalar>();
105 ComplexMatrix Cresult;
106 MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
108 result = Cresult.real();
112 typename internal::nested<MatrixType>::type m_A;
113 AtomicType& m_atomic;
122 template <
typename MatrixType,
typename AtomicType>
123 class MatrixFunction<MatrixType, AtomicType, 1>
127 typedef internal::traits<MatrixType> Traits;
128 typedef typename MatrixType::Scalar Scalar;
129 typedef typename MatrixType::Index Index;
130 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
131 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
132 static const int Options = MatrixType::Options;
133 typedef typename NumTraits<Scalar>::Real RealScalar;
134 typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
135 typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
136 typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
137 typedef std::list<Scalar> Cluster;
138 typedef std::list<Cluster> ListOfClusters;
139 typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
144 template <
typename ResultType>
void compute(ResultType& result);
148 void computeSchurDecomposition();
149 void partitionEigenvalues();
150 typename ListOfClusters::iterator findCluster(Scalar key);
151 void computeClusterSize();
152 void computeBlockStart();
153 void constructPermutation();
155 void swapEntriesInSchur(Index index);
156 void computeBlockAtomic();
157 Block<MatrixType> block(MatrixType& A, Index i, Index j);
158 void computeOffDiagonal();
159 DynMatrixType solveTriangularSylvester(
const DynMatrixType& A,
const DynMatrixType& B,
const DynMatrixType& C);
161 typename internal::nested<MatrixType>::type m_A;
162 AtomicType& m_atomic;
166 ListOfClusters m_clusters;
167 DynamicIntVectorType m_eivalToCluster;
168 DynamicIntVectorType m_clusterSize;
169 DynamicIntVectorType m_blockStart;
170 IntVectorType m_permutation;
178 static const RealScalar separation() {
return static_cast<RealScalar
>(0.1); }
188 template <
typename MatrixType,
typename AtomicType>
190 : m_A(A), m_atomic(atomic)
200 template <
typename MatrixType,
typename AtomicType>
201 template <
typename ResultType>
204 computeSchurDecomposition();
205 partitionEigenvalues();
206 computeClusterSize();
208 constructPermutation();
210 computeBlockAtomic();
211 computeOffDiagonal();
212 result = m_U * m_fT * m_U.adjoint();
216 template <
typename MatrixType,
typename AtomicType>
217 void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
219 const ComplexSchur<MatrixType> schurOfA(m_A);
220 m_T = schurOfA.matrixT();
221 m_U = schurOfA.matrixU();
235 template <
typename MatrixType,
typename AtomicType>
236 void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
238 const Index rows = m_T.rows();
239 VectorType diag = m_T.diagonal();
241 for (Index i=0; i<rows; ++i) {
243 typename ListOfClusters::iterator qi = findCluster(diag(i));
244 if (qi == m_clusters.end()) {
246 l.push_back(diag(i));
247 m_clusters.push_back(l);
248 qi = m_clusters.end();
253 for (Index j=i+1; j<rows; ++j) {
254 if (internal::abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
255 typename ListOfClusters::iterator qj = findCluster(diag(j));
256 if (qj == m_clusters.end()) {
257 qi->push_back(diag(j));
259 qi->insert(qi->end(), qj->begin(), qj->end());
260 m_clusters.erase(qj);
272 template <
typename MatrixType,
typename AtomicType>
273 typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
275 typename Cluster::iterator j;
276 for (
typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
277 j = std::find(i->begin(), i->end(), key);
281 return m_clusters.end();
285 template <
typename MatrixType,
typename AtomicType>
286 void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
288 const Index rows = m_T.rows();
289 VectorType diag = m_T.diagonal();
290 const Index numClusters =
static_cast<Index
>(m_clusters.size());
292 m_clusterSize.setZero(numClusters);
293 m_eivalToCluster.resize(rows);
294 Index clusterIndex = 0;
295 for (
typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
296 for (Index i = 0; i < diag.rows(); ++i) {
297 if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
298 ++m_clusterSize[clusterIndex];
299 m_eivalToCluster[i] = clusterIndex;
307 template <
typename MatrixType,
typename AtomicType>
308 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
310 m_blockStart.resize(m_clusterSize.rows());
312 for (Index i = 1; i < m_clusterSize.rows(); i++) {
313 m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
318 template <
typename MatrixType,
typename AtomicType>
319 void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
321 DynamicIntVectorType indexNextEntry = m_blockStart;
322 m_permutation.resize(m_T.rows());
323 for (Index i = 0; i < m_T.rows(); i++) {
324 Index cluster = m_eivalToCluster[i];
325 m_permutation[i] = indexNextEntry[cluster];
326 ++indexNextEntry[cluster];
331 template <
typename MatrixType,
typename AtomicType>
332 void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
334 IntVectorType p = m_permutation;
335 for (Index i = 0; i < p.rows() - 1; i++) {
337 for (j = i; j < p.rows(); j++) {
338 if (p(j) == i)
break;
340 eigen_assert(p(j) == i);
341 for (Index k = j-1; k >= i; k--) {
342 swapEntriesInSchur(k);
343 std::swap(p.coeffRef(k), p.coeffRef(k+1));
349 template <
typename MatrixType,
typename AtomicType>
350 void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
352 JacobiRotation<Scalar> rotation;
353 rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
354 m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
355 m_T.applyOnTheRight(index, index+1, rotation);
356 m_U.applyOnTheRight(index, index+1, rotation);
365 template <
typename MatrixType,
typename AtomicType>
366 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
368 m_fT.resize(m_T.rows(), m_T.cols());
370 for (Index i = 0; i < m_clusterSize.rows(); ++i) {
371 block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
376 template <
typename MatrixType,
typename AtomicType>
377 Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
379 return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
389 template <
typename MatrixType,
typename AtomicType>
390 void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
392 for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
393 for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
395 DynMatrixType A = block(m_T, blockIndex, blockIndex);
396 DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
397 DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
398 C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
399 for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
400 C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
401 C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
403 block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
431 template <
typename MatrixType,
typename AtomicType>
432 typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
433 const DynMatrixType& A,
434 const DynMatrixType& B,
435 const DynMatrixType& C)
437 eigen_assert(A.rows() == A.cols());
438 eigen_assert(A.isUpperTriangular());
439 eigen_assert(B.rows() == B.cols());
440 eigen_assert(B.isUpperTriangular());
441 eigen_assert(C.rows() == A.rows());
442 eigen_assert(C.cols() == B.rows());
446 DynMatrixType X(m, n);
448 for (Index i = m - 1; i >= 0; --i) {
449 for (Index j = 0; j < n; ++j) {
456 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
465 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
469 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
488 :
public ReturnByValue<MatrixFunctionReturnValue<Derived> >
492 typedef typename Derived::Scalar Scalar;
493 typedef typename Derived::Index Index;
494 typedef typename internal::stem_function<Scalar>::type StemFunction;
509 template <
typename ResultType>
510 inline void evalTo(ResultType& result)
const
512 typedef typename Derived::PlainObject PlainObject;
513 typedef internal::traits<PlainObject> Traits;
514 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
515 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
516 static const int Options = PlainObject::Options;
517 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
520 AtomicType atomic(m_f);
522 const PlainObject Aevaluated = m_A.eval();
527 Index rows()
const {
return m_A.rows(); }
528 Index cols()
const {
return m_A.cols(); }
531 typename internal::nested<Derived>::type m_A;
538 template<
typename Derived>
539 struct traits<MatrixFunctionReturnValue<Derived> >
541 typedef typename Derived::PlainObject ReturnType;
549 template <
typename Derived>
550 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(
typename internal::stem_function<
typename internal::traits<Derived>::Scalar>::type f)
const
552 eigen_assert(rows() == cols());
553 return MatrixFunctionReturnValue<Derived>(derived(), f);
556 template <
typename Derived>
557 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin()
const
559 eigen_assert(rows() == cols());
560 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
561 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
564 template <
typename Derived>
565 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos()
const
567 eigen_assert(rows() == cols());
568 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
569 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
572 template <
typename Derived>
573 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh()
const
575 eigen_assert(rows() == cols());
576 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
577 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
580 template <
typename Derived>
581 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh()
const
583 eigen_assert(rows() == cols());
584 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
585 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
590 #endif // EIGEN_MATRIX_FUNCTION