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

Collapse All | Expand All

(-)a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h (-4 / +20 lines)
Lines 294-299 Link Here
294
template<typename MatrixType>
294
template<typename MatrixType>
295
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
295
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
296
{
296
{
297
  res.resizeLike(m_A);
297
  switch (m_A.rows()) {
298
  switch (m_A.rows()) {
298
    case 0:
299
    case 0:
299
      break;
300
      break;
Lines 320-325 Link Here
320
  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
321
  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
321
}
322
}
322
323
324
// this function assumes that res has the correct size (see bug 614)
323
template<typename MatrixType>
325
template<typename MatrixType>
324
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
326
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
325
{
327
{
Lines 332-348 Link Here
332
  for (Index i=1; i < m_A.cols(); ++i) {
334
  for (Index i=1; i < m_A.cols(); ++i) {
333
    res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
335
    res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
334
    if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
336
    if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
335
      res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1);
337
      res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
336
    }
338
    }
337
    else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
339
    else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
338
      res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
340
      res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
339
    }
341
    }
340
    else {
342
    else {
341
      int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
343
      int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
342
      Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
344
      Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
343
      res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) *
345
      res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
344
	  std::sinh(p * w) / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
346
          / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
345
    }
347
    }
348
    res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
346
  }
349
  }
347
}
350
}
348
351
Lines 360-365 Link Here
360
  int degree, degree2, numberOfSquareRoots = 0;
363
  int degree, degree2, numberOfSquareRoots = 0;
361
  bool hasExtraSquareRoot = false;
364
  bool hasExtraSquareRoot = false;
362
365
366
  /* FIXME
367
   * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
368
   * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
369
   * worry about tiny values (e.g. 1e-300) because they will reach 1 if
370
   * repetitively sqrt'ed.
371
   *
372
   * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
373
   * [      ]   = [                   ]
374
   * [ 0  0 ]     [  0         0      ]
375
   */
376
  for (Index i=0; i < m_A.cols(); ++i)
377
    eigen_assert(m_A(i,i) != RealScalar(0))
378
363
  while (true) {
379
  while (true) {
364
    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
380
    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
365
    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
381
    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
366
382
(-)a/unsupported/test/matrix_power.cpp (-1 / +1 lines)
Lines 181-186 Link Here
181
  CALL_SUBTEST_1(testMatrixVector(Matrix2f(),         Vector2f(),    1e-4));
181
  CALL_SUBTEST_1(testMatrixVector(Matrix2f(),         Vector2f(),    1e-4));
182
  CALL_SUBTEST_5(testMatrixVector(Matrix3cf(),        Vector3cf(),   1e-4));
182
  CALL_SUBTEST_5(testMatrixVector(Matrix3cf(),        Vector3cf(),   1e-4));
183
  CALL_SUBTEST_8(testMatrixVector(Matrix4f(),         Vector4f(),    1e-4));
183
  CALL_SUBTEST_8(testMatrixVector(Matrix4f(),         Vector4f(),    1e-4));
184
  CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8),      VectorXf(8),   1e-3));
184
  CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2),      VectorXf(2),   1e-3)); // see bug 614
185
  CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7),      VectorXe(7),   1e-13));
185
  CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7),      VectorXe(7),   1e-13));
186
}
186
}

Return to bug 614