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 |
|