merge

Eigen/Jacobi

@@ -0,0 +1,24 @@
+#ifndef EIGEN_JACOBI_MODULE_H
+#define EIGEN_JACOBI_MODULE_H
+
+#include "Core"
+
+#include "src/Core/util/DisableMSVCWarnings.h"
+
+namespace Eigen {
+
+/** \defgroup Jacobi_Module Jacobi module
+  * This module provides Jacobi rotations.
+  *
+  * \code
+  * #include <Eigen/Jacobi>
+  * \endcode
+  */
+
+#include "src/Jacobi/Jacobi.h"
+
+} // namespace Eigen
+
+#include "src/Core/util/EnableMSVCWarnings.h"
+
+#endif // EIGEN_JACOBI_MODULE_H

Eigen/SVD

@@ -2,6 +2,8 @@
 #define EIGEN_SVD_MODULE_H
 
 #include "Core"
+#include "Householder"
+#include "Jacobi"
 
 #include "src/Core/util/DisableMSVCWarnings.h"
 
@@ -20,7 +22,9 @@ namespace Eigen {
   * \endcode
   */
 
+#include "src/SVD/Bidiagonalization.h"
 #include "src/SVD/SVD.h"
+#include "src/SVD/JacobiSquareSVD.h"
 
 } // namespace Eigen
 

Eigen/src/Core/Assign.h

@@ -57,7 +57,10 @@ private:
                   && ((int(Derived::Flags)&RowMajorBit)==(int(OtherDerived::Flags)&RowMajorBit)),
     MayInnerVectorize  = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0
                        && int(DstIsAligned) && int(SrcIsAligned),
-    MayLinearVectorize = MightVectorize && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit),
+    MayLinearVectorize = MightVectorize && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit)
+                                        && (DstIsAligned || InnerMaxSize == Dynamic),/* If the destination isn't aligned,
+      we have to do runtime checks and we don't unroll, so it's only good for large enough sizes. See remark below
+      about InnerMaxSize. */
     MaySliceVectorize  = MightVectorize && int(InnerMaxSize)>=3*PacketSize /* slice vectorization can be slow, so we only
       want it if the slices are big, which is indicated by InnerMaxSize rather than InnerSize, think of the case
       of a dynamic block in a fixed-size matrix */
@@ -90,6 +93,25 @@ public:
               ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) )
               : int(NoUnrolling)
   };
+  
+  static void debug()
+  {
+    EIGEN_DEBUG_VAR(DstIsAligned)
+    EIGEN_DEBUG_VAR(SrcIsAligned)
+    EIGEN_DEBUG_VAR(SrcAlignment)
+    EIGEN_DEBUG_VAR(InnerSize)
+    EIGEN_DEBUG_VAR(InnerMaxSize)
+    EIGEN_DEBUG_VAR(PacketSize)
+    EIGEN_DEBUG_VAR(MightVectorize)
+    EIGEN_DEBUG_VAR(MayInnerVectorize)
+    EIGEN_DEBUG_VAR(MayLinearVectorize)
+    EIGEN_DEBUG_VAR(MaySliceVectorize)
+    EIGEN_DEBUG_VAR(Vectorization)
+    EIGEN_DEBUG_VAR(UnrollingLimit)
+    EIGEN_DEBUG_VAR(MayUnrollCompletely)
+    EIGEN_DEBUG_VAR(MayUnrollInner)
+    EIGEN_DEBUG_VAR(Unrolling)    
+  }
 };
 
 /***************************************************************************
@@ -405,6 +427,9 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>
     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
   ei_assert(rows() == other.rows() && cols() == other.cols());
   ei_assign_impl<Derived, OtherDerived>::run(derived(),other.derived());
+#ifdef EIGEN_DEBUG_ASSIGN
+  ei_assign_traits<Derived, OtherDerived>::debug();
+#endif
   return derived();
 }
 

Eigen/src/Core/MatrixBase.h

@@ -778,6 +778,12 @@ template<typename Derived> class MatrixBase
     void applyHouseholderOnTheRight(const EssentialPart& essential,
                                     const RealScalar& beta);
 
+///////// Jacobi module /////////
+
+    void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s);
+    void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s);
+    bool makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s);
+    bool makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s);
 
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN

Eigen/src/Householder/Householder.h

@@ -32,6 +32,15 @@ template<int n> struct ei_decrement_size
   };
 };
 
+template<typename EssentialPart>
+void makeTrivialHouseholder(
+  EssentialPart *essential,
+  typename EssentialPart::RealScalar *beta)
+{
+  *beta = typename EssentialPart::RealScalar(0);
+  essential->setZero();
+}
+
 template<typename Derived>
 template<typename EssentialPart>
 void MatrixBase<Derived>::makeHouseholder(

Eigen/src/Jacobi/Jacobi.h

@@ -0,0 +1,91 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_JACOBI_H
+#define EIGEN_JACOBI_H
+
+template<typename Derived>
+void MatrixBase<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s)
+{
+  for(int i = 0; i < cols(); ++i)
+  {
+    Scalar tmp = coeff(p,i);
+    coeffRef(p,i) = c * tmp - s * coeff(q,i);
+    coeffRef(q,i) = s * tmp + c * coeff(q,i);
+  }
+}
+
+template<typename Derived>
+void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s)
+{
+  for(int i = 0; i < rows(); ++i)
+  {
+    Scalar tmp = coeff(i,p);
+    coeffRef(i,p) = c * tmp - s * coeff(i,q);
+    coeffRef(i,q) = s * tmp + c * coeff(i,q);
+  }
+}
+
+template<typename Scalar>
+bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Scalar *s)
+{
+  if(ei_abs(y) < max_coeff * 0.5 * machine_epsilon<Scalar>())
+  {
+    *c = Scalar(1);
+    *s = Scalar(0);
+    return true;
+  }
+  else
+  {
+    Scalar tau = (z - x) / (2 * y);
+    Scalar w = ei_sqrt(1 + ei_abs2(tau));
+    Scalar t;
+    if(tau>0)
+      t = Scalar(1) / (tau + w);
+    else
+      t = Scalar(1) / (tau - w);
+    *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t));
+    *s = *c * t;
+    return false;
+  }
+}
+
+template<typename Derived>
+inline bool MatrixBase<Derived>::makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s)
+{
+  return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), max_coeff, c, s);
+}
+
+template<typename Derived>
+inline bool MatrixBase<Derived>::makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s)
+{
+  return ei_makeJacobi(col(p).squaredNorm(),
+                       col(p).dot(col(q)),
+                       col(q).squaredNorm(),
+                       max_coeff,
+                       c,s);
+}
+
+
+#endif // EIGEN_JACOBI_H

Eigen/src/SVD/JacobiSquareSVD.h

@@ -0,0 +1,170 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_JACOBISQUARESVD_H
+#define EIGEN_JACOBISQUARESVD_H
+
+/** \ingroup SVD_Module
+  * \nonstableyet
+  *
+  * \class JacobiSquareSVD
+  *
+  * \brief Jacobi SVD decomposition of a square matrix
+  *
+  * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
+  * \param ComputeU whether the U matrix should be computed
+  * \param ComputeV whether the V matrix should be computed
+  *
+  * \sa MatrixBase::jacobiSvd()
+  */
+template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSVD
+{
+  private:
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+      Options = MatrixType::Options
+    };
+    
+    typedef Matrix<Scalar, Dynamic, Dynamic, Options> DummyMatrixType;
+    typedef typename ei_meta_if<ComputeU,
+                                Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
+                                       Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime>,
+                                DummyMatrixType>::ret MatrixUType;
+    typedef typename Diagonal<MatrixType,0>::PlainMatrixType SingularValuesType;
+    typedef Matrix<Scalar, 1, RowsAtCompileTime, Options, 1, MaxRowsAtCompileTime> RowType;
+    typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> ColType;
+
+  public:
+
+    JacobiSquareSVD() : m_isInitialized(false) {}
+
+    JacobiSquareSVD(const MatrixType& matrix) : m_isInitialized(false) 
+    {
+      compute(matrix);
+    }
+    
+    void compute(const MatrixType& matrix);
+    
+    const MatrixUType& matrixU() const
+    {
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_matrixU;
+    }
+
+    const SingularValuesType& singularValues() const
+    {
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_singularValues;
+    }
+
+    const MatrixUType& matrixV() const
+    {
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_matrixV;
+    }
+
+  protected:
+    MatrixUType m_matrixU;
+    MatrixUType m_matrixV;
+    SingularValuesType m_singularValues;
+    bool m_isInitialized;
+};
+
+template<typename MatrixType, bool ComputeU, bool ComputeV>
+void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix)
+{
+  MatrixType work_matrix(matrix);
+  int size = matrix.rows();
+  if(ComputeU) m_matrixU = MatrixUType::Identity(size,size);
+  if(ComputeV) m_matrixV = MatrixUType::Identity(size,size);
+  m_singularValues.resize(size);
+  RealScalar max_coeff = work_matrix.cwise().abs().maxCoeff();
+  for(int k = 1; k < 40; ++k) {
+    bool finished = true;
+    for(int p = 1; p < size; ++p)
+    {
+      for(int q = 0; q < p; ++q)
+      {
+        Scalar c, s;
+        finished &= work_matrix.makeJacobiForAtA(p,q,max_coeff,&c,&s);
+        work_matrix.applyJacobiOnTheRight(p,q,c,s);
+        if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s);
+      }
+    }
+    if(finished) break;
+  }
+  
+  for(int i = 0; i < size; ++i)
+  {
+    m_singularValues.coeffRef(i) = work_matrix.col(i).norm();
+  }
+
+  int first_zero = size;
+  RealScalar biggest = m_singularValues.maxCoeff();
+  for(int i = 0; i < size; i++)
+  {
+    int pos;
+    RealScalar biggest_remaining = m_singularValues.end(size-i).maxCoeff(&pos);
+    if(first_zero == size && ei_isMuchSmallerThan(biggest_remaining, biggest)) first_zero = pos + i;
+    if(pos)
+    {
+      pos += i;
+      std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
+      if(ComputeU) work_matrix.col(pos).swap(work_matrix.col(i));
+      if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i));
+    }
+  }
+  
+  if(ComputeU)
+  {
+    for(int i = 0; i < first_zero; ++i)
+    {
+      m_matrixU.col(i) = work_matrix.col(i) / m_singularValues.coeff(i);
+    }
+    if(first_zero < size)
+    {
+      for(int i = first_zero; i < size; ++i)
+      {
+        for(int j = 0; j < size; ++j)
+        {
+          m_matrixU.col(i).setZero();
+          m_matrixU.coeffRef(j,i) = Scalar(1);
+          for(int k = 0; k < first_zero; ++k)
+            m_matrixU.col(i) -= m_matrixU.col(i).dot(m_matrixU.col(k)) * m_matrixU.col(k);
+          RealScalar n = m_matrixU.col(i).norm();
+          if(!ei_isMuchSmallerThan(n, biggest))
+          {
+            m_matrixU.col(i) /= n;
+            break;
+          }
+        }
+      }
+    }     
+  }
+  m_isInitialized = true;
+}
+#endif // EIGEN_JACOBISQUARESVD_H

Eigen/src/SVD/SVD.h

@@ -393,8 +393,9 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
     {
       int k;
       W.end(n-i).minCoeff(&k);
-      if (k != i)
+      if (k != 0)
       {
+        k += i;
         std::swap(W[k],W[i]);
         A.col(i).swap(A.col(k));
         V.col(i).swap(V.col(k));