1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2009-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
---|
5 | //
|
---|
6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
9 |
|
---|
10 | #ifndef EIGEN_MATRIX_FUNCTION
|
---|
11 | #define EIGEN_MATRIX_FUNCTION
|
---|
12 |
|
---|
13 | #include "StemFunction.h"
|
---|
14 | #include "MatrixFunctionAtomic.h"
|
---|
15 |
|
---|
16 |
|
---|
17 | namespace Eigen {
|
---|
18 |
|
---|
19 | /** \ingroup MatrixFunctions_Module
|
---|
20 | * \brief Class for computing matrix functions.
|
---|
21 | * \tparam MatrixType type of the argument of the matrix function,
|
---|
22 | * expected to be an instantiation of the Matrix class template.
|
---|
23 | * \tparam AtomicType type for computing matrix function of atomic blocks.
|
---|
24 | * \tparam IsComplex used internally to select correct specialization.
|
---|
25 | *
|
---|
26 | * This class implements the Schur-Parlett algorithm for computing matrix functions. The spectrum of the
|
---|
27 | * matrix is divided in clustered of eigenvalues that lies close together. This class delegates the
|
---|
28 | * computation of the matrix function on every block corresponding to these clusters to an object of type
|
---|
29 | * \p AtomicType and uses these results to compute the matrix function of the whole matrix. The class
|
---|
30 | * \p AtomicType should have a \p compute() member function for computing the matrix function of a block.
|
---|
31 | *
|
---|
32 | * \sa class MatrixFunctionAtomic, class MatrixLogarithmAtomic
|
---|
33 | */
|
---|
34 | template <typename MatrixType,
|
---|
35 | typename AtomicType,
|
---|
36 | int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
|
---|
37 | class MatrixFunction
|
---|
38 | {
|
---|
39 | public:
|
---|
40 |
|
---|
41 | /** \brief Constructor.
|
---|
42 | *
|
---|
43 | * \param[in] A argument of matrix function, should be a square matrix.
|
---|
44 | * \param[in] atomic class for computing matrix function of atomic blocks.
|
---|
45 | *
|
---|
46 | * The class stores references to \p A and \p atomic, so they should not be
|
---|
47 | * changed (or destroyed) before compute() is called.
|
---|
48 | */
|
---|
49 | MatrixFunction(const MatrixType& A, AtomicType& atomic);
|
---|
50 |
|
---|
51 | /** \brief Compute the matrix function.
|
---|
52 | *
|
---|
53 | * \param[out] result the function \p f applied to \p A, as
|
---|
54 | * specified in the constructor.
|
---|
55 | *
|
---|
56 | * See MatrixBase::matrixFunction() for details on how this computation
|
---|
57 | * is implemented.
|
---|
58 | */
|
---|
59 | template <typename ResultType>
|
---|
60 | void compute(ResultType &result);
|
---|
61 | };
|
---|
62 |
|
---|
63 |
|
---|
64 | /** \internal \ingroup MatrixFunctions_Module
|
---|
65 | * \brief Partial specialization of MatrixFunction for real matrices
|
---|
66 | */
|
---|
67 | template <typename MatrixType, typename AtomicType>
|
---|
68 | class MatrixFunction<MatrixType, AtomicType, 0>
|
---|
69 | {
|
---|
70 | private:
|
---|
71 |
|
---|
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;
|
---|
79 |
|
---|
80 | typedef std::complex<Scalar> ComplexScalar;
|
---|
81 | typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
|
---|
82 |
|
---|
83 | public:
|
---|
84 |
|
---|
85 | /** \brief Constructor.
|
---|
86 | *
|
---|
87 | * \param[in] A argument of matrix function, should be a square matrix.
|
---|
88 | * \param[in] atomic class for computing matrix function of atomic blocks.
|
---|
89 | */
|
---|
90 | MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
|
---|
91 |
|
---|
92 | /** \brief Compute the matrix function.
|
---|
93 | *
|
---|
94 | * \param[out] result the function \p f applied to \p A, as
|
---|
95 | * specified in the constructor.
|
---|
96 | *
|
---|
97 | * This function converts the real matrix \c A to a complex matrix,
|
---|
98 | * uses MatrixFunction<MatrixType,1> and then converts the result back to
|
---|
99 | * a real matrix.
|
---|
100 | */
|
---|
101 | template <typename ResultType>
|
---|
102 | void compute(ResultType& result)
|
---|
103 | {
|
---|
104 | ComplexMatrix CA = m_A.template cast<ComplexScalar>();
|
---|
105 | ComplexMatrix Cresult;
|
---|
106 | MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
|
---|
107 | mf.compute(Cresult);
|
---|
108 | result = Cresult.real();
|
---|
109 | }
|
---|
110 |
|
---|
111 | private:
|
---|
112 | typename internal::nested<MatrixType>::type m_A; /**< \brief Reference to argument of matrix function. */
|
---|
113 | AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */
|
---|
114 |
|
---|
115 | MatrixFunction& operator=(const MatrixFunction&);
|
---|
116 | };
|
---|
117 |
|
---|
118 |
|
---|
119 | /** \internal \ingroup MatrixFunctions_Module
|
---|
120 | * \brief Partial specialization of MatrixFunction for complex matrices
|
---|
121 | */
|
---|
122 | template <typename MatrixType, typename AtomicType>
|
---|
123 | class MatrixFunction<MatrixType, AtomicType, 1>
|
---|
124 | {
|
---|
125 | private:
|
---|
126 |
|
---|
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;
|
---|
140 |
|
---|
141 | public:
|
---|
142 |
|
---|
143 | MatrixFunction(const MatrixType& A, AtomicType& atomic);
|
---|
144 | template <typename ResultType> void compute(ResultType& result);
|
---|
145 |
|
---|
146 | private:
|
---|
147 |
|
---|
148 | void computeSchurDecomposition();
|
---|
149 | void partitionEigenvalues();
|
---|
150 | typename ListOfClusters::iterator findCluster(Scalar key);
|
---|
151 | void computeClusterSize();
|
---|
152 | void computeBlockStart();
|
---|
153 | void constructPermutation();
|
---|
154 | void permuteSchur();
|
---|
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);
|
---|
160 |
|
---|
161 | typename internal::nested<MatrixType>::type m_A; /**< \brief Reference to argument of matrix function. */
|
---|
162 | AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */
|
---|
163 | MatrixType m_T; /**< \brief Triangular part of Schur decomposition */
|
---|
164 | MatrixType m_U; /**< \brief Unitary part of Schur decomposition */
|
---|
165 | MatrixType m_fT; /**< \brief %Matrix function applied to #m_T */
|
---|
166 | ListOfClusters m_clusters; /**< \brief Partition of eigenvalues into clusters of ei'vals "close" to each other */
|
---|
167 | DynamicIntVectorType m_eivalToCluster; /**< \brief m_eivalToCluster[i] = j means i-th ei'val is in j-th cluster */
|
---|
168 | DynamicIntVectorType m_clusterSize; /**< \brief Number of eigenvalues in each clusters */
|
---|
169 | DynamicIntVectorType m_blockStart; /**< \brief Row index at which block corresponding to i-th cluster starts */
|
---|
170 | IntVectorType m_permutation; /**< \brief Permutation which groups ei'vals in the same cluster together */
|
---|
171 |
|
---|
172 | /** \brief Maximum distance allowed between eigenvalues to be considered "close".
|
---|
173 | *
|
---|
174 | * This is morally a \c static \c const \c Scalar, but only
|
---|
175 | * integers can be static constant class members in C++. The
|
---|
176 | * separation constant is set to 0.1, a value taken from the
|
---|
177 | * paper by Davies and Higham. */
|
---|
178 | static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
|
---|
179 |
|
---|
180 | MatrixFunction& operator=(const MatrixFunction&);
|
---|
181 | };
|
---|
182 |
|
---|
183 | /** \brief Constructor.
|
---|
184 | *
|
---|
185 | * \param[in] A argument of matrix function, should be a square matrix.
|
---|
186 | * \param[in] atomic class for computing matrix function of atomic blocks.
|
---|
187 | */
|
---|
188 | template <typename MatrixType, typename AtomicType>
|
---|
189 | MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
|
---|
190 | : m_A(A), m_atomic(atomic)
|
---|
191 | {
|
---|
192 | /* empty body */
|
---|
193 | }
|
---|
194 |
|
---|
195 | /** \brief Compute the matrix function.
|
---|
196 | *
|
---|
197 | * \param[out] result the function \p f applied to \p A, as
|
---|
198 | * specified in the constructor.
|
---|
199 | */
|
---|
200 | template <typename MatrixType, typename AtomicType>
|
---|
201 | template <typename ResultType>
|
---|
202 | void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result)
|
---|
203 | {
|
---|
204 | computeSchurDecomposition();
|
---|
205 | partitionEigenvalues();
|
---|
206 | computeClusterSize();
|
---|
207 | computeBlockStart();
|
---|
208 | constructPermutation();
|
---|
209 | permuteSchur();
|
---|
210 | computeBlockAtomic();
|
---|
211 | computeOffDiagonal();
|
---|
212 | result = m_U * (m_fT.template triangularView<Upper>() * m_U.adjoint());
|
---|
213 | }
|
---|
214 |
|
---|
215 | /** \brief Store the Schur decomposition of #m_A in #m_T and #m_U */
|
---|
216 | template <typename MatrixType, typename AtomicType>
|
---|
217 | void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
|
---|
218 | {
|
---|
219 | const ComplexSchur<MatrixType> schurOfA(m_A);
|
---|
220 | m_T = schurOfA.matrixT();
|
---|
221 | m_U = schurOfA.matrixU();
|
---|
222 | }
|
---|
223 |
|
---|
224 | /** \brief Partition eigenvalues in clusters of ei'vals close to each other
|
---|
225 | *
|
---|
226 | * This function computes #m_clusters. This is a partition of the
|
---|
227 | * eigenvalues of #m_T in clusters, such that
|
---|
228 | * # Any eigenvalue in a certain cluster is at most separation() away
|
---|
229 | * from another eigenvalue in the same cluster.
|
---|
230 | * # The distance between two eigenvalues in different clusters is
|
---|
231 | * more than separation().
|
---|
232 | * The implementation follows Algorithm 4.1 in the paper of Davies
|
---|
233 | * and Higham.
|
---|
234 | */
|
---|
235 | template <typename MatrixType, typename AtomicType>
|
---|
236 | void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
|
---|
237 | {
|
---|
238 | using std::abs;
|
---|
239 | const Index rows = m_T.rows();
|
---|
240 | VectorType diag = m_T.diagonal(); // contains eigenvalues of A
|
---|
241 |
|
---|
242 | for (Index i=0; i<rows; ++i) {
|
---|
243 | // Find set containing diag(i), adding a new set if necessary
|
---|
244 | typename ListOfClusters::iterator qi = findCluster(diag(i));
|
---|
245 | if (qi == m_clusters.end()) {
|
---|
246 | Cluster l;
|
---|
247 | l.push_back(diag(i));
|
---|
248 | m_clusters.push_back(l);
|
---|
249 | qi = m_clusters.end();
|
---|
250 | --qi;
|
---|
251 | }
|
---|
252 |
|
---|
253 | // Look for other element to add to the set
|
---|
254 | for (Index j=i+1; j<rows; ++j) {
|
---|
255 | if (abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
|
---|
256 | typename ListOfClusters::iterator qj = findCluster(diag(j));
|
---|
257 | if (qj == m_clusters.end()) {
|
---|
258 | qi->push_back(diag(j));
|
---|
259 | } else {
|
---|
260 | qi->insert(qi->end(), qj->begin(), qj->end());
|
---|
261 | m_clusters.erase(qj);
|
---|
262 | }
|
---|
263 | }
|
---|
264 | }
|
---|
265 | }
|
---|
266 | }
|
---|
267 |
|
---|
268 | /** \brief Find cluster in #m_clusters containing some value
|
---|
269 | * \param[in] key Value to find
|
---|
270 | * \returns Iterator to cluster containing \c key, or
|
---|
271 | * \c m_clusters.end() if no cluster in m_clusters contains \c key.
|
---|
272 | */
|
---|
273 | template <typename MatrixType, typename AtomicType>
|
---|
274 | typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
|
---|
275 | {
|
---|
276 | typename Cluster::iterator j;
|
---|
277 | for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
|
---|
278 | j = std::find(i->begin(), i->end(), key);
|
---|
279 | if (j != i->end())
|
---|
280 | return i;
|
---|
281 | }
|
---|
282 | return m_clusters.end();
|
---|
283 | }
|
---|
284 |
|
---|
285 | /** \brief Compute #m_clusterSize and #m_eivalToCluster using #m_clusters */
|
---|
286 | template <typename MatrixType, typename AtomicType>
|
---|
287 | void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
|
---|
288 | {
|
---|
289 | const Index rows = m_T.rows();
|
---|
290 | VectorType diag = m_T.diagonal();
|
---|
291 | const Index numClusters = static_cast<Index>(m_clusters.size());
|
---|
292 |
|
---|
293 | m_clusterSize.setZero(numClusters);
|
---|
294 | m_eivalToCluster.resize(rows);
|
---|
295 | Index clusterIndex = 0;
|
---|
296 | for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
|
---|
297 | for (Index i = 0; i < diag.rows(); ++i) {
|
---|
298 | if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
|
---|
299 | ++m_clusterSize[clusterIndex];
|
---|
300 | m_eivalToCluster[i] = clusterIndex;
|
---|
301 | }
|
---|
302 | }
|
---|
303 | ++clusterIndex;
|
---|
304 | }
|
---|
305 | }
|
---|
306 |
|
---|
307 | /** \brief Compute #m_blockStart using #m_clusterSize */
|
---|
308 | template <typename MatrixType, typename AtomicType>
|
---|
309 | void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
|
---|
310 | {
|
---|
311 | m_blockStart.resize(m_clusterSize.rows());
|
---|
312 | m_blockStart(0) = 0;
|
---|
313 | for (Index i = 1; i < m_clusterSize.rows(); i++) {
|
---|
314 | m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
|
---|
315 | }
|
---|
316 | }
|
---|
317 |
|
---|
318 | /** \brief Compute #m_permutation using #m_eivalToCluster and #m_blockStart */
|
---|
319 | template <typename MatrixType, typename AtomicType>
|
---|
320 | void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
|
---|
321 | {
|
---|
322 | DynamicIntVectorType indexNextEntry = m_blockStart;
|
---|
323 | m_permutation.resize(m_T.rows());
|
---|
324 | for (Index i = 0; i < m_T.rows(); i++) {
|
---|
325 | Index cluster = m_eivalToCluster[i];
|
---|
326 | m_permutation[i] = indexNextEntry[cluster];
|
---|
327 | ++indexNextEntry[cluster];
|
---|
328 | }
|
---|
329 | }
|
---|
330 |
|
---|
331 | /** \brief Permute Schur decomposition in #m_U and #m_T according to #m_permutation */
|
---|
332 | template <typename MatrixType, typename AtomicType>
|
---|
333 | void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
|
---|
334 | {
|
---|
335 | IntVectorType p = m_permutation;
|
---|
336 | for (Index i = 0; i < p.rows() - 1; i++) {
|
---|
337 | Index j;
|
---|
338 | for (j = i; j < p.rows(); j++) {
|
---|
339 | if (p(j) == i) break;
|
---|
340 | }
|
---|
341 | eigen_assert(p(j) == i);
|
---|
342 | for (Index k = j-1; k >= i; k--) {
|
---|
343 | swapEntriesInSchur(k);
|
---|
344 | std::swap(p.coeffRef(k), p.coeffRef(k+1));
|
---|
345 | }
|
---|
346 | }
|
---|
347 | }
|
---|
348 |
|
---|
349 | /** \brief Swap rows \a index and \a index+1 in Schur decomposition in #m_U and #m_T */
|
---|
350 | template <typename MatrixType, typename AtomicType>
|
---|
351 | void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
|
---|
352 | {
|
---|
353 | JacobiRotation<Scalar> rotation;
|
---|
354 | rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
|
---|
355 | m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
|
---|
356 | m_T.applyOnTheRight(index, index+1, rotation);
|
---|
357 | m_U.applyOnTheRight(index, index+1, rotation);
|
---|
358 | }
|
---|
359 |
|
---|
360 | /** \brief Compute block diagonal part of #m_fT.
|
---|
361 | *
|
---|
362 | * This routine computes the matrix function applied to the block diagonal part of #m_T, with the blocking
|
---|
363 | * given by #m_blockStart. The matrix function of each diagonal block is computed by #m_atomic. The
|
---|
364 | * off-diagonal parts of #m_fT are set to zero.
|
---|
365 | */
|
---|
366 | template <typename MatrixType, typename AtomicType>
|
---|
367 | void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
|
---|
368 | {
|
---|
369 | m_fT.resize(m_T.rows(), m_T.cols());
|
---|
370 | m_fT.setZero();
|
---|
371 | for (Index i = 0; i < m_clusterSize.rows(); ++i) {
|
---|
372 | block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
|
---|
373 | }
|
---|
374 | }
|
---|
375 |
|
---|
376 | /** \brief Return block of matrix according to blocking given by #m_blockStart */
|
---|
377 | template <typename MatrixType, typename AtomicType>
|
---|
378 | Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
|
---|
379 | {
|
---|
380 | return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
|
---|
381 | }
|
---|
382 |
|
---|
383 | /** \brief Compute part of #m_fT above block diagonal.
|
---|
384 | *
|
---|
385 | * This routine assumes that the block diagonal part of #m_fT (which
|
---|
386 | * equals the matrix function applied to #m_T) has already been computed and computes
|
---|
387 | * the part above the block diagonal. The part below the diagonal is
|
---|
388 | * zero, because #m_T is upper triangular.
|
---|
389 | */
|
---|
390 | template <typename MatrixType, typename AtomicType>
|
---|
391 | void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
|
---|
392 | {
|
---|
393 | for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
|
---|
394 | for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
|
---|
395 | // compute (blockIndex, blockIndex+diagIndex) block
|
---|
396 | DynMatrixType A = block(m_T, blockIndex, blockIndex);
|
---|
397 | DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
|
---|
398 | DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
|
---|
399 | C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
|
---|
400 | for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
|
---|
401 | C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
|
---|
402 | C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
|
---|
403 | }
|
---|
404 | block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
|
---|
405 | }
|
---|
406 | }
|
---|
407 | }
|
---|
408 |
|
---|
409 | /** \brief Solve a triangular Sylvester equation AX + XB = C
|
---|
410 | *
|
---|
411 | * \param[in] A the matrix A; should be square and upper triangular
|
---|
412 | * \param[in] B the matrix B; should be square and upper triangular
|
---|
413 | * \param[in] C the matrix C; should have correct size.
|
---|
414 | *
|
---|
415 | * \returns the solution X.
|
---|
416 | *
|
---|
417 | * If A is m-by-m and B is n-by-n, then both C and X are m-by-n.
|
---|
418 | * The (i,j)-th component of the Sylvester equation is
|
---|
419 | * \f[
|
---|
420 | * \sum_{k=i}^m A_{ik} X_{kj} + \sum_{k=1}^j X_{ik} B_{kj} = C_{ij}.
|
---|
421 | * \f]
|
---|
422 | * This can be re-arranged to yield:
|
---|
423 | * \f[
|
---|
424 | * X_{ij} = \frac{1}{A_{ii} + B_{jj}} \Bigl( C_{ij}
|
---|
425 | * - \sum_{k=i+1}^m A_{ik} X_{kj} - \sum_{k=1}^{j-1} X_{ik} B_{kj} \Bigr).
|
---|
426 | * \f]
|
---|
427 | * It is assumed that A and B are such that the numerator is never
|
---|
428 | * zero (otherwise the Sylvester equation does not have a unique
|
---|
429 | * solution). In that case, these equations can be evaluated in the
|
---|
430 | * order \f$ i=m,\ldots,1 \f$ and \f$ j=1,\ldots,n \f$.
|
---|
431 | */
|
---|
432 | template <typename MatrixType, typename AtomicType>
|
---|
433 | typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
|
---|
434 | const DynMatrixType& A,
|
---|
435 | const DynMatrixType& B,
|
---|
436 | const DynMatrixType& C)
|
---|
437 | {
|
---|
438 | eigen_assert(A.rows() == A.cols());
|
---|
439 | eigen_assert(A.isUpperTriangular());
|
---|
440 | eigen_assert(B.rows() == B.cols());
|
---|
441 | eigen_assert(B.isUpperTriangular());
|
---|
442 | eigen_assert(C.rows() == A.rows());
|
---|
443 | eigen_assert(C.cols() == B.rows());
|
---|
444 |
|
---|
445 | Index m = A.rows();
|
---|
446 | Index n = B.rows();
|
---|
447 | DynMatrixType X(m, n);
|
---|
448 |
|
---|
449 | for (Index i = m - 1; i >= 0; --i) {
|
---|
450 | for (Index j = 0; j < n; ++j) {
|
---|
451 |
|
---|
452 | // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
|
---|
453 | Scalar AX;
|
---|
454 | if (i == m - 1) {
|
---|
455 | AX = 0;
|
---|
456 | } else {
|
---|
457 | Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
|
---|
458 | AX = AXmatrix(0,0);
|
---|
459 | }
|
---|
460 |
|
---|
461 | // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
|
---|
462 | Scalar XB;
|
---|
463 | if (j == 0) {
|
---|
464 | XB = 0;
|
---|
465 | } else {
|
---|
466 | Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
|
---|
467 | XB = XBmatrix(0,0);
|
---|
468 | }
|
---|
469 |
|
---|
470 | X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
|
---|
471 | }
|
---|
472 | }
|
---|
473 | return X;
|
---|
474 | }
|
---|
475 |
|
---|
476 | /** \ingroup MatrixFunctions_Module
|
---|
477 | *
|
---|
478 | * \brief Proxy for the matrix function of some matrix (expression).
|
---|
479 | *
|
---|
480 | * \tparam Derived Type of the argument to the matrix function.
|
---|
481 | *
|
---|
482 | * This class holds the argument to the matrix function until it is
|
---|
483 | * assigned or evaluated for some other reason (so the argument
|
---|
484 | * should not be changed in the meantime). It is the return type of
|
---|
485 | * matrixBase::matrixFunction() and related functions and most of the
|
---|
486 | * time this is the only way it is used.
|
---|
487 | */
|
---|
488 | template<typename Derived> class MatrixFunctionReturnValue
|
---|
489 | : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
|
---|
490 | {
|
---|
491 | public:
|
---|
492 |
|
---|
493 | typedef typename Derived::Scalar Scalar;
|
---|
494 | typedef typename Derived::Index Index;
|
---|
495 | typedef typename internal::stem_function<Scalar>::type StemFunction;
|
---|
496 |
|
---|
497 | /** \brief Constructor.
|
---|
498 | *
|
---|
499 | * \param[in] A %Matrix (expression) forming the argument of the
|
---|
500 | * matrix function.
|
---|
501 | * \param[in] f Stem function for matrix function under consideration.
|
---|
502 | */
|
---|
503 | MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
|
---|
504 |
|
---|
505 | /** \brief Compute the matrix function.
|
---|
506 | *
|
---|
507 | * \param[out] result \p f applied to \p A, where \p f and \p A
|
---|
508 | * are as in the constructor.
|
---|
509 | */
|
---|
510 | template <typename ResultType>
|
---|
511 | inline void evalTo(ResultType& result) const
|
---|
512 | {
|
---|
513 | typedef typename Derived::PlainObject PlainObject;
|
---|
514 | typedef internal::traits<PlainObject> Traits;
|
---|
515 | static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
|
---|
516 | static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
|
---|
517 | static const int Options = PlainObject::Options;
|
---|
518 | typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
|
---|
519 | typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
|
---|
520 | typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
|
---|
521 | AtomicType atomic(m_f);
|
---|
522 |
|
---|
523 | const PlainObject Aevaluated = m_A.eval();
|
---|
524 | MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
|
---|
525 | mf.compute(result);
|
---|
526 | }
|
---|
527 |
|
---|
528 | Index rows() const { return m_A.rows(); }
|
---|
529 | Index cols() const { return m_A.cols(); }
|
---|
530 |
|
---|
531 | private:
|
---|
532 | typename internal::nested<Derived>::type m_A;
|
---|
533 | StemFunction *m_f;
|
---|
534 |
|
---|
535 | MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&);
|
---|
536 | };
|
---|
537 |
|
---|
538 | namespace internal {
|
---|
539 | template<typename Derived>
|
---|
540 | struct traits<MatrixFunctionReturnValue<Derived> >
|
---|
541 | {
|
---|
542 | typedef typename Derived::PlainObject ReturnType;
|
---|
543 | };
|
---|
544 | }
|
---|
545 |
|
---|
546 |
|
---|
547 | /********** MatrixBase methods **********/
|
---|
548 |
|
---|
549 |
|
---|
550 | template <typename Derived>
|
---|
551 | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
|
---|
552 | {
|
---|
553 | eigen_assert(rows() == cols());
|
---|
554 | return MatrixFunctionReturnValue<Derived>(derived(), f);
|
---|
555 | }
|
---|
556 |
|
---|
557 | template <typename Derived>
|
---|
558 | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
|
---|
559 | {
|
---|
560 | eigen_assert(rows() == cols());
|
---|
561 | typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
|
---|
562 | return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
|
---|
563 | }
|
---|
564 |
|
---|
565 | template <typename Derived>
|
---|
566 | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
|
---|
567 | {
|
---|
568 | eigen_assert(rows() == cols());
|
---|
569 | typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
|
---|
570 | return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
|
---|
571 | }
|
---|
572 |
|
---|
573 | template <typename Derived>
|
---|
574 | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
|
---|
575 | {
|
---|
576 | eigen_assert(rows() == cols());
|
---|
577 | typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
|
---|
578 | return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
|
---|
579 | }
|
---|
580 |
|
---|
581 | template <typename Derived>
|
---|
582 | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
|
---|
583 | {
|
---|
584 | eigen_assert(rows() == cols());
|
---|
585 | typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
|
---|
586 | return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
|
---|
587 | }
|
---|
588 |
|
---|
589 | } // end namespace Eigen
|
---|
590 |
|
---|
591 | #endif // EIGEN_MATRIX_FUNCTION
|
---|