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/EigenSolver.h (-54 / +3 lines)
Lines 380-388 Link Here
380
{
380
{
381
  check_template_parameters();
381
  check_template_parameters();
382
  
382
  
383
  using std::sqrt;
384
  using std::abs;
385
  using numext::isfinite;
386
  eigen_assert(matrix.cols() == matrix.rows());
383
  eigen_assert(matrix.cols() == matrix.rows());
387
384
388
  // Reduce to real Schur form.
385
  // Reduce to real Schur form.
Lines 393-455 Link Here
393
  if (m_info == Success)
390
  if (m_info == Success)
394
  {
391
  {
395
    m_matT = m_realSchur.matrixT();
392
    m_matT = m_realSchur.matrixT();
393
    m_eivalues = m_realSchur.eigenvalues();
396
    if (computeEigenvectors)
394
    if (computeEigenvectors)
397
      m_eivec = m_realSchur.matrixU();
398
  
399
    // Compute eigenvalues from matT
400
    m_eivalues.resize(matrix.cols());
401
    Index i = 0;
402
    while (i < matrix.cols()) 
403
    {
395
    {
404
      if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
396
      m_eivec = m_realSchur.matrixU();
405
      {
397
      doComputeEigenvectors();
406
        m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
407
        if(!(isfinite)(m_eivalues.coeffRef(i)))
408
        {
409
          m_isInitialized = true;
410
          m_eigenvectorsOk = false;
411
          m_info = NumericalIssue;
412
          return *this;
413
        }
414
        ++i;
415
      }
416
      else
417
      {
418
        Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
419
        Scalar z;
420
        // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
421
        // without overflow
422
        {
423
          Scalar t0 = m_matT.coeff(i+1, i);
424
          Scalar t1 = m_matT.coeff(i, i+1);
425
          Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
426
          t0 /= maxval;
427
          t1 /= maxval;
428
          Scalar p0 = p/maxval;
429
          z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
430
        }
431
        
432
        m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
433
        m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
434
        if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
435
        {
436
          m_isInitialized = true;
437
          m_eigenvectorsOk = false;
438
          m_info = NumericalIssue;
439
          return *this;
440
        }
441
        i += 2;
442
      }
443
    }
398
    }
444
    
445
    // Compute eigenvectors.
446
    if (computeEigenvectors)
447
      doComputeEigenvectors();
448
  }
399
  }
449
450
  m_isInitialized = true;
400
  m_isInitialized = true;
451
  m_eigenvectorsOk = computeEigenvectors;
401
  m_eigenvectorsOk = computeEigenvectors;
452
453
  return *this;
402
  return *this;
454
}
403
}
455
404

Return to bug 1535