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

(-)a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h (-82 / +101 lines)
Lines 586-723 template<typename SolverType> struct dir Link Here
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<VectorType> res, Ref<VectorType> 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<Scalar>::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<Scalar>::epsilon();
617
      if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
673
      if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::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<Lower>();
680
        scaledMat = scaledMat.template selfadjointView<Lower>();
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<Scalar>::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<Scalar>::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

Return to bug 1014