This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 1535 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/Eigenvalues/RealSchur.h (-5 / +91 lines)
Lines 83-88 Link Here
83
    explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
83
    explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
84
            : m_matT(size, size),
84
            : m_matT(size, size),
85
              m_matU(size, size),
85
              m_matU(size, size),
86
              m_eigenvalues(size),
86
              m_workspaceVector(size),
87
              m_workspaceVector(size),
87
              m_hess(size),
88
              m_hess(size),
88
              m_isInitialized(false),
89
              m_isInitialized(false),
Lines 93-101 Link Here
93
    /** \brief Constructor; computes real Schur decomposition of given matrix. 
94
    /** \brief Constructor; computes real Schur decomposition of given matrix. 
94
      * 
95
      * 
95
      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
96
      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
96
      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
97
      * \param[in]  computeU  If true, both T and U are computed; if false, only T and the eigenvalues are computed.
97
      *
98
      *
98
      * This constructor calls compute() to compute the Schur decomposition.
99
      * This constructor calls compute() to compute the Schur decomposition and the
100
      * eigenvalues of the matrix.
99
      *
101
      *
100
      * Example: \include RealSchur_RealSchur_MatrixType.cpp
102
      * Example: \include RealSchur_RealSchur_MatrixType.cpp
101
      * Output: \verbinclude RealSchur_RealSchur_MatrixType.out
103
      * Output: \verbinclude RealSchur_RealSchur_MatrixType.out
Lines 147-156 Link Here
147
      return m_matT;
149
      return m_matT;
148
    }
150
    }
149
  
151
  
152
    /** \brief Returns the eigenvalues of the given matrix.
153
      *
154
      * \returns A column vector containing the eigenvalues.
155
      *
156
      * \pre Either the constructor RealSchur(const MatrixType&, bool) or the
157
      * member function compute(const MatrixType&, bool) has been called before
158
      * to compute the Schur decomposition of a matrix.
159
      *
160
      * The eigenvalues of the matrix are the same as the eigenvalues of the 1x1
161
      * and 2x2 diagonal blocks of the matrix T and are extracted from that matrix
162
      * after a factorization has been computed.
163
      *
164
      * \sa RealSchur(const MatrixType&, bool) for an example
165
      */
166
    const EigenvalueType& eigenvalues() const
167
    {
168
      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
169
      return m_eigenvalues;
170
    }
171
150
    /** \brief Computes Schur decomposition of given matrix. 
172
    /** \brief Computes Schur decomposition of given matrix. 
151
      * 
173
      * 
152
      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
174
      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
153
      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
175
      * \param[in]  computeU  If true, both T and U are computed; if false, only T and the eigenvalues are computed.
154
      * \returns    Reference to \c *this
176
      * \returns    Reference to \c *this
155
      *
177
      *
156
      * The Schur decomposition is computed by first reducing the matrix to
178
      * The Schur decomposition is computed by first reducing the matrix to
Lines 178-185 Link Here
178
     *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
200
     *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
179
     *  using either the class HessenbergDecomposition or another mean. 
201
     *  using either the class HessenbergDecomposition or another mean. 
180
     *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
202
     *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
203
     *  and calculates the eigenvalues of T.
181
     *  When computeU is true, this routine computes the matrix U such that 
204
     *  When computeU is true, this routine computes the matrix U such that 
182
     *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
205
     *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix.
183
     * 
206
     * 
184
     * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
207
     * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
185
     * is not available, the user should give an identity matrix (Q.setIdentity())
208
     * is not available, the user should give an identity matrix (Q.setIdentity())
Lines 226-231 Link Here
226
    
249
    
227
    MatrixType m_matT;
250
    MatrixType m_matT;
228
    MatrixType m_matU;
251
    MatrixType m_matU;
252
    EigenvalueType m_eigenvalues;
229
    ColumnVectorType m_workspaceVector;
253
    ColumnVectorType m_workspaceVector;
230
    HessenbergDecomposition<MatrixType> m_hess;
254
    HessenbergDecomposition<MatrixType> m_hess;
231
    ComputationInfo m_info;
255
    ComputationInfo m_info;
Lines 241-246 Link Here
241
    void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
265
    void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
242
    void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
266
    void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
243
    void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
267
    void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
268
    /** \brief Resizes vector m_eigenvalues and populates it with the
269
     *         eigenvalues of quasi-triangular matrix m_T.
270
     *  \return NumericalIssue if a non-finite eigenvalue was detected, Success otherwise
271
     */
272
    ComputationInfo extractEigenvalues();
244
};
273
};
245
274
246
275
Lines 259-264 Link Here
259
  if(scale<considerAsZero)
288
  if(scale<considerAsZero)
260
  {
289
  {
261
    m_matT.setZero(matrix.rows(),matrix.cols());
290
    m_matT.setZero(matrix.rows(),matrix.cols());
291
    m_eigenvalues.setZero(matrix.rows());
262
    if(computeU)
292
    if(computeU)
263
      m_matU.setIdentity(matrix.rows(),matrix.cols());
293
      m_matU.setIdentity(matrix.rows(),matrix.cols());
264
    m_info = Success;
294
    m_info = Success;
Lines 274-282 Link Here
274
  computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
304
  computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
275
305
276
  m_matT *= scale;
306
  m_matT *= scale;
277
  
307
  m_eigenvalues *= scale;
308
278
  return *this;
309
  return *this;
279
}
310
}
311
280
template<typename MatrixType>
312
template<typename MatrixType>
281
template<typename HessMatrixType, typename OrthMatrixType>
313
template<typename HessMatrixType, typename OrthMatrixType>
282
RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU)
314
RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU)
Lines 343-348 Link Here
343
  else
375
  else
344
    m_info = NoConvergence;
376
    m_info = NoConvergence;
345
377
378
  ComputationInfo ev_info = extractEigenvalues();
379
  // If the factorization failed, we keep the original message (NoConvergence...)
380
  // Otherwise, set m_info to the result of extractEigenvalues(),
381
  // which may return NumericalIssue if a NaN/infinity is found.
382
  if (m_info == Success)
383
    m_info = ev_info;
384
346
  m_isInitialized = true;
385
  m_isInitialized = true;
347
  m_matUisUptodate = computeU;
386
  m_matUisUptodate = computeU;
348
  return *this;
387
  return *this;
Lines 541-546 Link Here
541
  }
580
  }
542
}
581
}
543
582
583
template<typename MatrixType>
584
ComputationInfo RealSchur<MatrixType>::extractEigenvalues()
585
{
586
  using std::sqrt;
587
  using std::abs;
588
  using numext::isfinite;
589
  m_eigenvalues.setZero(m_matT.cols());
590
  Index i = 0;
591
  while (i < m_matT.cols())
592
  {
593
    if (i == m_matT.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
594
    {
595
      m_eigenvalues.coeffRef(i) = m_matT.coeff(i, i);
596
      if(!(isfinite)(m_eigenvalues.coeffRef(i)))
597
      {
598
        return NumericalIssue;
599
      }
600
      ++i;
601
    }
602
    else
603
    {
604
      Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
605
      Scalar z;
606
      // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
607
      // without overflow
608
      {
609
        Scalar t0 = m_matT.coeff(i+1, i);
610
        Scalar t1 = m_matT.coeff(i, i+1);
611
        Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
612
        t0 /= maxval;
613
        t1 /= maxval;
614
        Scalar p0 = p/maxval;
615
        z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
616
      }
617
      m_eigenvalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
618
      m_eigenvalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
619
      if(!((isfinite)(m_eigenvalues.coeffRef(i)) && (isfinite)(m_eigenvalues.coeffRef(i+1))))
620
      {
621
        return NumericalIssue;
622
      }
623
      i += 2;
624
    }
625
  }
626
  return Success;
627
}
628
629
544
} // end namespace Eigen
630
} // end namespace Eigen
545
631
546
#endif // EIGEN_REAL_SCHUR_H
632
#endif // EIGEN_REAL_SCHUR_H
(-)a/doc/snippets/RealSchur_RealSchur_MatrixType.cpp (+1 lines)
Lines 4-9 Link Here
4
RealSchur<MatrixXd> schur(A);
4
RealSchur<MatrixXd> schur(A);
5
cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl;
5
cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl;
6
cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl;
6
cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl;
7
cout << "The eigenvalues are:" << endl << schur.eigenvalues() << endl << endl;
7
8
8
MatrixXd U = schur.matrixU();
9
MatrixXd U = schur.matrixU();
9
MatrixXd T = schur.matrixT();
10
MatrixXd T = schur.matrixT();

Return to bug 1535