00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_MATRIX_FUNCTION
00026 #define EIGEN_MATRIX_FUNCTION
00027
00028 #include "StemFunction.h"
00029 #include "MatrixFunctionAtomic.h"
00030
00031
00032 namespace Eigen {
00033
00049 template <typename MatrixType,
00050 typename AtomicType,
00051 int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
00052 class MatrixFunction
00053 {
00054 public:
00055
00064 MatrixFunction(const MatrixType& A, AtomicType& atomic);
00065
00074 template <typename ResultType>
00075 void compute(ResultType &result);
00076 };
00077
00078
00082 template <typename MatrixType, typename AtomicType>
00083 class MatrixFunction<MatrixType, AtomicType, 0>
00084 {
00085 private:
00086
00087 typedef internal::traits<MatrixType> Traits;
00088 typedef typename Traits::Scalar Scalar;
00089 static const int Rows = Traits::RowsAtCompileTime;
00090 static const int Cols = Traits::ColsAtCompileTime;
00091 static const int Options = MatrixType::Options;
00092 static const int MaxRows = Traits::MaxRowsAtCompileTime;
00093 static const int MaxCols = Traits::MaxColsAtCompileTime;
00094
00095 typedef std::complex<Scalar> ComplexScalar;
00096 typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
00097
00098 public:
00099
00105 MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
00106
00116 template <typename ResultType>
00117 void compute(ResultType& result)
00118 {
00119 ComplexMatrix CA = m_A.template cast<ComplexScalar>();
00120 ComplexMatrix Cresult;
00121 MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
00122 mf.compute(Cresult);
00123 result = Cresult.real();
00124 }
00125
00126 private:
00127 typename internal::nested<MatrixType>::type m_A;
00128 AtomicType& m_atomic;
00130 MatrixFunction& operator=(const MatrixFunction&);
00131 };
00132
00133
00137 template <typename MatrixType, typename AtomicType>
00138 class MatrixFunction<MatrixType, AtomicType, 1>
00139 {
00140 private:
00141
00142 typedef internal::traits<MatrixType> Traits;
00143 typedef typename MatrixType::Scalar Scalar;
00144 typedef typename MatrixType::Index Index;
00145 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00146 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00147 static const int Options = MatrixType::Options;
00148 typedef typename NumTraits<Scalar>::Real RealScalar;
00149 typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
00150 typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
00151 typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
00152 typedef std::list<Scalar> Cluster;
00153 typedef std::list<Cluster> ListOfClusters;
00154 typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00155
00156 public:
00157
00158 MatrixFunction(const MatrixType& A, AtomicType& atomic);
00159 template <typename ResultType> void compute(ResultType& result);
00160
00161 private:
00162
00163 void computeSchurDecomposition();
00164 void partitionEigenvalues();
00165 typename ListOfClusters::iterator findCluster(Scalar key);
00166 void computeClusterSize();
00167 void computeBlockStart();
00168 void constructPermutation();
00169 void permuteSchur();
00170 void swapEntriesInSchur(Index index);
00171 void computeBlockAtomic();
00172 Block<MatrixType> block(MatrixType& A, Index i, Index j);
00173 void computeOffDiagonal();
00174 DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C);
00175
00176 typename internal::nested<MatrixType>::type m_A;
00177 AtomicType& m_atomic;
00178 MatrixType m_T;
00179 MatrixType m_U;
00180 MatrixType m_fT;
00181 ListOfClusters m_clusters;
00182 DynamicIntVectorType m_eivalToCluster;
00183 DynamicIntVectorType m_clusterSize;
00184 DynamicIntVectorType m_blockStart;
00185 IntVectorType m_permutation;
00193 static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
00194
00195 MatrixFunction& operator=(const MatrixFunction&);
00196 };
00197
00203 template <typename MatrixType, typename AtomicType>
00204 MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
00205 : m_A(A), m_atomic(atomic)
00206 {
00207
00208 }
00209
00215 template <typename MatrixType, typename AtomicType>
00216 template <typename ResultType>
00217 void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result)
00218 {
00219 computeSchurDecomposition();
00220 partitionEigenvalues();
00221 computeClusterSize();
00222 computeBlockStart();
00223 constructPermutation();
00224 permuteSchur();
00225 computeBlockAtomic();
00226 computeOffDiagonal();
00227 result = m_U * m_fT * m_U.adjoint();
00228 }
00229
00231 template <typename MatrixType, typename AtomicType>
00232 void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
00233 {
00234 const ComplexSchur<MatrixType> schurOfA(m_A);
00235 m_T = schurOfA.matrixT();
00236 m_U = schurOfA.matrixU();
00237 }
00238
00250 template <typename MatrixType, typename AtomicType>
00251 void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
00252 {
00253 const Index rows = m_T.rows();
00254 VectorType diag = m_T.diagonal();
00255
00256 for (Index i=0; i<rows; ++i) {
00257
00258 typename ListOfClusters::iterator qi = findCluster(diag(i));
00259 if (qi == m_clusters.end()) {
00260 Cluster l;
00261 l.push_back(diag(i));
00262 m_clusters.push_back(l);
00263 qi = m_clusters.end();
00264 --qi;
00265 }
00266
00267
00268 for (Index j=i+1; j<rows; ++j) {
00269 if (internal::abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
00270 typename ListOfClusters::iterator qj = findCluster(diag(j));
00271 if (qj == m_clusters.end()) {
00272 qi->push_back(diag(j));
00273 } else {
00274 qi->insert(qi->end(), qj->begin(), qj->end());
00275 m_clusters.erase(qj);
00276 }
00277 }
00278 }
00279 }
00280 }
00281
00287 template <typename MatrixType, typename AtomicType>
00288 typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
00289 {
00290 typename Cluster::iterator j;
00291 for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
00292 j = std::find(i->begin(), i->end(), key);
00293 if (j != i->end())
00294 return i;
00295 }
00296 return m_clusters.end();
00297 }
00298
00300 template <typename MatrixType, typename AtomicType>
00301 void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
00302 {
00303 const Index rows = m_T.rows();
00304 VectorType diag = m_T.diagonal();
00305 const Index numClusters = static_cast<Index>(m_clusters.size());
00306
00307 m_clusterSize.setZero(numClusters);
00308 m_eivalToCluster.resize(rows);
00309 Index clusterIndex = 0;
00310 for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
00311 for (Index i = 0; i < diag.rows(); ++i) {
00312 if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
00313 ++m_clusterSize[clusterIndex];
00314 m_eivalToCluster[i] = clusterIndex;
00315 }
00316 }
00317 ++clusterIndex;
00318 }
00319 }
00320
00322 template <typename MatrixType, typename AtomicType>
00323 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
00324 {
00325 m_blockStart.resize(m_clusterSize.rows());
00326 m_blockStart(0) = 0;
00327 for (Index i = 1; i < m_clusterSize.rows(); i++) {
00328 m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
00329 }
00330 }
00331
00333 template <typename MatrixType, typename AtomicType>
00334 void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
00335 {
00336 DynamicIntVectorType indexNextEntry = m_blockStart;
00337 m_permutation.resize(m_T.rows());
00338 for (Index i = 0; i < m_T.rows(); i++) {
00339 Index cluster = m_eivalToCluster[i];
00340 m_permutation[i] = indexNextEntry[cluster];
00341 ++indexNextEntry[cluster];
00342 }
00343 }
00344
00346 template <typename MatrixType, typename AtomicType>
00347 void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
00348 {
00349 IntVectorType p = m_permutation;
00350 for (Index i = 0; i < p.rows() - 1; i++) {
00351 Index j;
00352 for (j = i; j < p.rows(); j++) {
00353 if (p(j) == i) break;
00354 }
00355 eigen_assert(p(j) == i);
00356 for (Index k = j-1; k >= i; k--) {
00357 swapEntriesInSchur(k);
00358 std::swap(p.coeffRef(k), p.coeffRef(k+1));
00359 }
00360 }
00361 }
00362
00364 template <typename MatrixType, typename AtomicType>
00365 void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
00366 {
00367 JacobiRotation<Scalar> rotation;
00368 rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
00369 m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
00370 m_T.applyOnTheRight(index, index+1, rotation);
00371 m_U.applyOnTheRight(index, index+1, rotation);
00372 }
00373
00380 template <typename MatrixType, typename AtomicType>
00381 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
00382 {
00383 m_fT.resize(m_T.rows(), m_T.cols());
00384 m_fT.setZero();
00385 for (Index i = 0; i < m_clusterSize.rows(); ++i) {
00386 block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
00387 }
00388 }
00389
00391 template <typename MatrixType, typename AtomicType>
00392 Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
00393 {
00394 return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
00395 }
00396
00404 template <typename MatrixType, typename AtomicType>
00405 void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
00406 {
00407 for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
00408 for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
00409
00410 DynMatrixType A = block(m_T, blockIndex, blockIndex);
00411 DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
00412 DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
00413 C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
00414 for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
00415 C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
00416 C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
00417 }
00418 block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
00419 }
00420 }
00421 }
00422
00446 template <typename MatrixType, typename AtomicType>
00447 typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
00448 const DynMatrixType& A,
00449 const DynMatrixType& B,
00450 const DynMatrixType& C)
00451 {
00452 eigen_assert(A.rows() == A.cols());
00453 eigen_assert(A.isUpperTriangular());
00454 eigen_assert(B.rows() == B.cols());
00455 eigen_assert(B.isUpperTriangular());
00456 eigen_assert(C.rows() == A.rows());
00457 eigen_assert(C.cols() == B.rows());
00458
00459 Index m = A.rows();
00460 Index n = B.rows();
00461 DynMatrixType X(m, n);
00462
00463 for (Index i = m - 1; i >= 0; --i) {
00464 for (Index j = 0; j < n; ++j) {
00465
00466
00467 Scalar AX;
00468 if (i == m - 1) {
00469 AX = 0;
00470 } else {
00471 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
00472 AX = AXmatrix(0,0);
00473 }
00474
00475
00476 Scalar XB;
00477 if (j == 0) {
00478 XB = 0;
00479 } else {
00480 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
00481 XB = XBmatrix(0,0);
00482 }
00483
00484 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
00485 }
00486 }
00487 return X;
00488 }
00489
00502 template<typename Derived> class MatrixFunctionReturnValue
00503 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
00504 {
00505 public:
00506
00507 typedef typename Derived::Scalar Scalar;
00508 typedef typename Derived::Index Index;
00509 typedef typename internal::stem_function<Scalar>::type StemFunction;
00510
00517 MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
00518
00524 template <typename ResultType>
00525 inline void evalTo(ResultType& result) const
00526 {
00527 typedef typename Derived::PlainObject PlainObject;
00528 typedef internal::traits<PlainObject> Traits;
00529 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00530 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00531 static const int Options = PlainObject::Options;
00532 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
00533 typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00534 typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
00535 AtomicType atomic(m_f);
00536
00537 const PlainObject Aevaluated = m_A.eval();
00538 MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
00539 mf.compute(result);
00540 }
00541
00542 Index rows() const { return m_A.rows(); }
00543 Index cols() const { return m_A.cols(); }
00544
00545 private:
00546 typename internal::nested<Derived>::type m_A;
00547 StemFunction *m_f;
00548
00549 MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&);
00550 };
00551
00552 namespace internal {
00553 template<typename Derived>
00554 struct traits<MatrixFunctionReturnValue<Derived> >
00555 {
00556 typedef typename Derived::PlainObject ReturnType;
00557 };
00558 }
00559
00560
00561
00562
00563
00564 template <typename Derived>
00565 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
00566 {
00567 eigen_assert(rows() == cols());
00568 return MatrixFunctionReturnValue<Derived>(derived(), f);
00569 }
00570
00571 template <typename Derived>
00572 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
00573 {
00574 eigen_assert(rows() == cols());
00575 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00576 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
00577 }
00578
00579 template <typename Derived>
00580 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
00581 {
00582 eigen_assert(rows() == cols());
00583 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00584 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
00585 }
00586
00587 template <typename Derived>
00588 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
00589 {
00590 eigen_assert(rows() == cols());
00591 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00592 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
00593 }
00594
00595 template <typename Derived>
00596 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
00597 {
00598 eigen_assert(rows() == cols());
00599 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00600 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
00601 }
00602
00603 }
00604
00605 #endif // EIGEN_MATRIX_FUNCTION