New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
View | Details | Raw Unified | Return to bug 1014 | Differences between
and this patch

Collapse All | Expand All

 Lines 586-723 template struct dir (-)a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h (-82 / +101 lines) 586 ` {` 586 ` {` 587 ` numext::swap(roots(1),roots(2));` 587 ` numext::swap(roots(1),roots(2));` 588 ` if (roots(0) >= roots(1))` 588 ` if (roots(0) >= roots(1))` 589 ` numext::swap(roots(0),roots(1));` 589 ` numext::swap(roots(0),roots(1));` 590 ` }` 590 ` }` 591 ` }` 591 ` }` 592 ` ` 592 ` ` 593 ` EIGEN_DEVICE_FUNC` 593 ` EIGEN_DEVICE_FUNC` 594 ` static inline bool extract_kernel(MatrixType& mat, Ref res, Ref representative)` 595 ` {` 596 ` using std::abs;` 597 ` // Compute the kernel of 'tmp' using a LU factorization of tmp with symmetric pivoting.` 598 ` // Note that we only need to compute the factor U, which is of the form:` 599 ` // x ? ?` 600 ` // 0 x ?` 601 ` // 0 0 0` 602 ` // With 'x' != 0.` 603 ` Index i0, i1, i2;` 604 ` Scalar p0, p1, p2, f;` 605 ` // Find the first pivot of index i0 (by construction, there must exist a non zero coefficient on the diagonal):` 606 ` mat.diagonal().cwiseAbs().maxCoeff(&i0);` 607 ` // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector,` 608 ` // so let's save it:` 609 ` representative = mat.col(i0);` 610 ` p0 = Scalar(1)/mat.coeff(i0,i0);` 611 ` // The other two rows/columns are:` 612 ` i1 = (i0+1)%3;` 613 ` i2 = (i0+2)%3;` 614 ` // Update the first row of U:` 615 ` mat.coeffRef(i0,i1) *= p0;` 616 ` mat.coeffRef(i0,i2) *= p0;` 617 ` ` 618 ` // Thanks to symmetry, the remaining 2x2 block is of the form:` 619 ` // p1 f` 620 ` // f p2` 621 ` // with:` 622 ` p1 = mat.coeffRef(i1,i1) - mat.coeff(i1,i0) * mat.coeff(i0,i1);` 623 ` p2 = mat.coeffRef(i2,i2) - mat.coeff(i2,i0) * mat.coeff(i0,i2);` 624 ` f = mat.coeffRef(i1,i2) - mat.coeff(i1,i0) * mat.coeff(i0,i2);` 625 ` ` 626 ` // Second pivoting:` 627 ` if(abs(p1) < abs(p2))` 628 ` {` 629 ` std::swap(i1,i2);` 630 ` p1 = p2;` 631 ` }` 632 ` // Make sure we have a meaningful pivot:` 633 ` if(abs(p1)<=abs(p0)*Eigen::NumTraits::epsilon())` 634 ` {` 635 ` // `mat` appears to be rank 1 only` 636 ` return false;` 637 ` }` 638 ` ` 639 ` // Back-substitutions to form the 1D kernel of 'mat':` 640 ` res.coeffRef(i2) = 1;` 641 ` res.coeffRef(i1) = - f / p1;` 642 ` res.coeffRef(i0) = - (mat.coeff(i0,i2) + mat.coeff(i0,i1) * res.coeff(i1) );` 643 ` res.normalize();` 644 ` ` 645 ` return true;` 646 ` }` 647 ` ` 648 ` EIGEN_DEVICE_FUNC` 594 ` static inline void run(SolverType& solver, const MatrixType& mat, int options)` 649 ` static inline void run(SolverType& solver, const MatrixType& mat, int options)` 595 ` {` 650 ` {` 596 ` using std::sqrt;` 597 ` eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());` 651 ` eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());` 598 ` eigen_assert((options&~(EigVecMask|GenEigMask))==0` 652 ` eigen_assert((options&~(EigVecMask|GenEigMask))==0` 599 ` && (options&EigVecMask)!=EigVecMask` 653 ` && (options&EigVecMask)!=EigVecMask` 600 ` && "invalid option parameter");` 654 ` && "invalid option parameter");` 601 ` bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;` 655 ` bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;` 602 ` ` 656 ` ` 603 ` MatrixType& eivecs = solver.m_eivec;` 657 ` MatrixType& eivecs = solver.m_eivec;` 604 ` VectorType& eivals = solver.m_eivalues;` 658 ` VectorType& eivals = solver.m_eivalues;` 605 ` ` 659 ` ` 606 ` // map the matrix coefficients to [-1:1] to avoid over- and underflow.` 660 ` // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.` 607 ` Scalar scale = mat.cwiseAbs().maxCoeff();` 661 ` Scalar shift = mat.trace() / Scalar(3);` 608 ` MatrixType scaledMat = mat / scale;` 662 ` MatrixType scaledMat = mat;` 663 ` scaledMat.diagonal().array() -= shift;` 664 ` Scalar scale = scaledMat.cwiseAbs().maxCoeff();` 665 ` if(scale > 0) scaledMat /= scale;` 609 666 610 ` // compute the eigenvalues` 667 ` // compute the eigenvalues` 611 ` computeRoots(scaledMat,eivals);` 668 ` computeRoots(scaledMat,eivals);` 612 669 613 ` // compute the eigen vectors` 670 ` // compute the eigenvectors` 614 ` if(computeEigenvectors)` 671 ` if(computeEigenvectors)` 615 ` {` 672 ` {` 616 ` Scalar safeNorm2 = Eigen::NumTraits::epsilon();` 617 ` if((eivals(2)-eivals(0))<=Eigen::NumTraits::epsilon())` 673 ` if((eivals(2)-eivals(0))<=Eigen::NumTraits::epsilon())` 618 ` {` 674 ` {` 675 ` // All three eigenvalues are numerically the same` 619 ` eivecs.setIdentity();` 676 ` eivecs.setIdentity();` 620 ` }` 677 ` }` 621 ` else` 678 ` else` 622 ` {` 679 ` {` 623 ` scaledMat = scaledMat.template selfadjointView();` 680 ` scaledMat = scaledMat.template selfadjointView();` 624 ` MatrixType tmp;` 681 ` MatrixType tmp;` 625 ` tmp = scaledMat;` 682 ` tmp = scaledMat;` 626 683 684 ` // Compute the eigenvector of the most distinct eigenvalue` 627 ` Scalar d0 = eivals(2) - eivals(1);` 685 ` Scalar d0 = eivals(2) - eivals(1);` 628 ` Scalar d1 = eivals(1) - eivals(0);` 686 ` Scalar d1 = eivals(1) - eivals(0);` 629 ` int k = d0 > d1 ? 2 : 0;` 687 ` Index k(0), l(2);` 630 ` d0 = d0 > d1 ? d0 : d1;` 688 ` if(d0 > d1)` 631 632 ` tmp.diagonal().array () -= eivals(k);` 633 ` VectorType cross;` 634 ` Scalar n;` 635 ` n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();` 636 637 ` if(n>safeNorm2)` 638 ` {` 689 ` {` 639 ` eivecs.col(k) = cross / sqrt(n);` 690 ` std::swap(k,l);` 691 ` d0 = d1;` 692 ` }` 693 ` ` 694 ` // Compute the eigenvector of index k` 695 ` {` 696 ` tmp.diagonal().array () -= eivals(k);` 697 ` // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector.` 698 ` if(!extract_kernel(tmp, eivecs.col(k), eivecs.col(l)))` 699 ` {` 700 ` // Such a case should be impossible, the only reason would be that the three eigenvalues` 701 ` // are very close to each other but not close enough to be detected earlier.` 702 ` eivecs.setIdentity();` 703 ` }` 704 ` }` 705 ` ` 706 ` // Compute eigenvector of index l` 707 ` if(d0<=2*Eigen::NumTraits::epsilon()*(eivals(2)-eivals(0)))` 708 ` {` 709 ` // If d0 is too small, then the two other eigenvalues are numerically the same,` 710 ` // and thus we only have to normalize the orthogonal vector we saved above.` 711 ` eivecs.col(l).normalize();` 640 ` }` 712 ` }` 641 ` else` 713 ` else` 642 ` {` 714 ` {` 643 ` n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();` 715 ` tmp = scaledMat;` 644 716 ` tmp.diagonal().array () -= eivals(l);` 645 ` if(n>safeNorm2)` 717 ` ` 718 ` VectorType dummy;` 719 ` if(!extract_kernel(tmp, eivecs.col(l), dummy))` 646 ` {` 720 ` {` 647 ` eivecs.col(k) = cross / sqrt(n);` 721 ` // Such a case should have been detected earlier.` 648 ` }` 722 ` // 'tmp' appears to be rank one, so let's peak an arbitrary orthogonal unit vector.` 649 ` else` 723 ` eivecs.col(l).normalize();` 650 ` {` 651 ` n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();` 652 653 ` if(n>safeNorm2)` 654 ` {` 655 ` eivecs.col(k) = cross / sqrt(n);` 656 ` }` 657 ` else` 658 ` {` 659 ` // the input matrix and/or the eigenvaues probably contains some inf/NaN,` 660 ` // => exit` 661 ` // scale back to the original size.` 662 ` eivals *= scale;` 663 664 ` solver.m_info = NumericalIssue;` 665 ` solver.m_isInitialized = true;` 666 ` solver.m_eigenvectorsOk = computeEigenvectors;` 667 ` return;` 668 ` }` 669 ` }` 724 ` }` 670 ` }` 725 ` }` 671 726 ` ` 672 ` tmp = scaledMat;` 727 ` // Compute last eigenvector from the other two` 673 ` tmp.diagonal().array() -= eivals(1);` 728 ` eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();` 674 675 ` if(d0<=Eigen::NumTraits::epsilon())` 676 ` {` 677 ` eivecs.col(1) = eivecs.col(k).unitOrthogonal();` 678 ` }` 679 ` else` 680 ` {` 681 ` n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm();` 682 ` if(n>safeNorm2)` 683 ` {` 684 ` eivecs.col(1) = cross / sqrt(n);` 685 ` }` 686 ` else` 687 ` {` 688 ` n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();` 689 ` if(n>safeNorm2)` 690 ` eivecs.col(1) = cross / sqrt(n);` 691 ` else` 692 ` {` 693 ` n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();` 694 ` if(n>safeNorm2)` 695 ` eivecs.col(1) = cross / sqrt(n);` 696 ` else` 697 ` {` 698 ` // we should never reach this point,` 699 ` // if so the last two eigenvalues are likely to be very close to each other` 700 ` eivecs.col(1) = eivecs.col(k).unitOrthogonal();` 701 ` }` 702 ` }` 703 ` }` 704 705 ` // make sure that eivecs[1] is orthogonal to eivecs[2]` 706 ` // FIXME: this step should not be needed` 707 ` Scalar d = eivecs.col(1).dot(eivecs.col(k));` 708 ` eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();` 709 ` }` 710 711 ` eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();` 712 ` }` 729 ` }` 713 ` }` 730 ` }` 731 ` ` 714 ` // Rescale back to the original size.` 732 ` // Rescale back to the original size.` 715 ` eivals *= scale;` 733 ` eivals *= scale;` 734 ` eivals.array() += shift;` 716 ` ` 735 ` ` 717 ` solver.m_info = Success;` 736 ` solver.m_info = Success;` 718 ` solver.m_isInitialized = true;` 737 ` solver.m_isInitialized = true;` 719 ` solver.m_eigenvectorsOk = computeEigenvectors;` 738 ` solver.m_eigenvectorsOk = computeEigenvectors;` 720 ` }` 739 ` }` 721 `};` 740 `};` 722 741 723 `// 2x2 direct eigenvalues decomposition, code from Hauke Heibel` 742 `// 2x2 direct eigenvalues decomposition, code from Hauke Heibel`