synch with main devel branch

46 files changed, 1219 insertions, 944 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 87e7e9a01..2bec68f9c 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -88,6 +88,8 @@
 #include <cstring>
 #include <string>
 #include <limits>
+// for min/max:
+#include <algorithm>
 
 #if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(EIGEN_NO_EXCEPTIONS)
   #define EIGEN_EXCEPTIONS
@@ -162,6 +164,7 @@ namespace Eigen {
 #include "src/Core/MapBase.h"
 #include "src/Core/Map.h"
 #include "src/Core/Block.h"
+#include "src/Core/VectorBlock.h"
 #include "src/Core/Minor.h"
 #include "src/Core/Transpose.h"
 #include "src/Core/DiagonalMatrix.h"
@@ -182,7 +185,6 @@ namespace Eigen {
 #include "src/Core/SolveTriangular.h"
 #include "src/Core/products/SelfadjointRank2Update.h"
 #include "src/Core/products/TriangularMatrixVector.h"
-#include "src/Core/BandMatrix.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/Geometry b/Eigen/Geometry
index 7860d8fca..9cae3459c 100644
--- a/Eigen/Geometry
+++ b/Eigen/Geometry
@@ -6,6 +6,7 @@
 #include "src/Core/util/DisableMSVCWarnings.h"
 
 #include "Array"
+#include "SVD"
 #include <limits>
 
 #ifndef M_PI
diff --git a/Eigen/Sparse b/Eigen/Sparse
index 364fd50c9..a8888daa3 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -84,6 +84,7 @@ namespace Eigen {
 
 #include "src/Sparse/SparseUtil.h"
 #include "src/Sparse/SparseMatrixBase.h"
+#include "src/Sparse/SparseNestByValue.h"
 #include "src/Sparse/CompressedStorage.h"
 #include "src/Sparse/AmbiVector.h"
 #include "src/Sparse/RandomSetter.h"
diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h
index df3afbbdb..02f9c0601 100644
--- a/Eigen/src/Array/Replicate.h
+++ b/Eigen/src/Array/Replicate.h
@@ -140,21 +140,4 @@ VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const
           (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
 }
 
-/** \nonstableyet
-  * \return an expression of the replication of each column (or row) of \c *this
-  *
-  * Example: \include DirectionWise_replicate.cpp
-  * Output: \verbinclude DirectionWise_replicate.out
-  *
-  * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate
-  */
-template<typename ExpressionType, int Direction>
-template<int Factor>
-const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)>
-VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const
-{
-  return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1>
-          (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
-}
-
 #endif // EIGEN_REPLICATE_H
diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h
index 50302bba4..3684d7cd1 100644
--- a/Eigen/src/Array/VectorwiseOp.h
+++ b/Eigen/src/Array/VectorwiseOp.h
@@ -179,6 +179,11 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
                               > Type;
     };
 
+    enum {
+      IsVertical   = (Direction==Vertical) ? 1 : 0,
+      IsHorizontal = (Direction==Horizontal) ? 1 : 0
+    };
+
   protected:
 
     /** \internal
@@ -222,9 +227,17 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
     /** \internal */
     inline const ExpressionType& _expression() const { return m_matrix; }
 
+    /** \returns a row or column vector expression of \c *this reduxed by \a func
+      *
+      * The template parameter \a BinaryOp is the type of the functor
+      * of the custom redux operator. Note that func must be an associative operator.
+      *
+      * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise()
+      */
     template<typename BinaryOp>
     const typename ReduxReturnType<BinaryOp>::Type
-    redux(const BinaryOp& func = BinaryOp()) const;
+    redux(const BinaryOp& func = BinaryOp()) const
+    { return typename ReduxReturnType<BinaryOp>::Type(_expression(), func); }
 
     /** \returns a row (or column) vector expression of the smallest coefficient
       * of each column (or row) of the referenced expression.
@@ -319,16 +332,26 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
       *
       * \sa MatrixBase::reverse() */
     const Reverse<ExpressionType, Direction> reverse() const
-    {
-      return Reverse<ExpressionType, Direction>( _expression() );
-    }
+    { return Reverse<ExpressionType, Direction>( _expression() ); }
 
     const Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1>
     replicate(int factor) const;
 
-    template<int Factor>
-    const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)>
-    replicate(int factor = Factor) const;
+    /** \nonstableyet
+      * \return an expression of the replication of each column (or row) of \c *this
+      *
+      * Example: \include DirectionWise_replicate.cpp
+      * Output: \verbinclude DirectionWise_replicate.out
+      *
+      * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate
+      */
+    // NOTE implemented here because of sunstudio's compilation errors
+    template<int Factor> const Replicate<ExpressionType,(IsVertical?Factor:1),(IsHorizontal?Factor:1)>
+    replicate(int factor = Factor) const
+    {
+      return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1>
+          (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
+    }
 
 /////////// Artithmetic operators ///////////
 
@@ -466,19 +489,4 @@ MatrixBase<Derived>::rowwise()
   return derived();
 }
 
-/** \returns a row or column vector expression of \c *this reduxed by \a func
-  *
-  * The template parameter \a BinaryOp is the type of the functor
-  * of the custom redux operator. Note that func must be an associative operator.
-  *
-  * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise()
-  */
-template<typename ExpressionType, int Direction>
-template<typename BinaryOp>
-const typename VectorwiseOp<ExpressionType,Direction>::template ReduxReturnType<BinaryOp>::Type
-VectorwiseOp<ExpressionType,Direction>::redux(const BinaryOp& func) const
-{
-  return typename ReduxReturnType<BinaryOp>::Type(_expression(), func);
-}
-
 #endif // EIGEN_PARTIAL_REDUX_H
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index de4268344..382518696 100644
--- a/Eigen/src/Core/Block.h
+++ b/Eigen/src/Core/Block.h
@@ -271,13 +271,19 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
 
     inline int stride(void) const { return m_matrix.stride(); }
 
+  #ifndef __SUNPRO_CC
+  // FIXME sunstudio is not friendly with the above friend...
   protected:
+  #endif
 
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal used by allowAligned() */
     inline Block(const MatrixType& matrix, const Scalar* data, int blockRows, int blockCols)
       : Base(data, blockRows, blockCols), m_matrix(matrix)
     {}
+    #endif
 
+  protected:
     const typename MatrixType::Nested m_matrix;
 };
 
@@ -314,249 +320,6 @@ inline const typename BlockReturnType<Derived>::Type MatrixBase<Derived>
   return typename BlockReturnType<Derived>::Type(derived(), startRow, startCol, blockRows, blockCols);
 }
 
-/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this.
-  *
-  * \only_for_vectors
-  *
-  * \addexample SegmentIntInt \label How to reference a sub-vector (dynamic size)
-  *
-  * \param start the first coefficient in the segment
-  * \param size the number of coefficients in the segment
-  *
-  * Example: \include MatrixBase_segment_int_int.cpp
-  * Output: \verbinclude MatrixBase_segment_int_int.out
-  *
-  * \note Even though the returned expression has dynamic size, in the case
-  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
-  * which means that evaluating it does not cause a dynamic memory allocation.
-  *
-  * \sa class Block, segment(int)
-  */
-template<typename Derived>
-inline typename BlockReturnType<Derived>::SubVectorType MatrixBase<Derived>
-  ::segment(int start, int size)
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start,
-                                   ColsAtCompileTime == 1 ? 0 : start,
-                                   RowsAtCompileTime == 1 ? 1 : size,
-                                   ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** This is the const version of segment(int,int).*/
-template<typename Derived>
-inline const typename BlockReturnType<Derived>::SubVectorType
-MatrixBase<Derived>::segment(int start, int size) const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start,
-                                   ColsAtCompileTime == 1 ? 0 : start,
-                                   RowsAtCompileTime == 1 ? 1 : size,
-                                   ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** \returns a dynamic-size expression of the first coefficients of *this.
-  *
-  * \only_for_vectors
-  *
-  * \param size the number of coefficients in the block
-  *
-  * \addexample BlockInt \label How to reference a sub-vector (fixed-size)
-  *
-  * Example: \include MatrixBase_start_int.cpp
-  * Output: \verbinclude MatrixBase_start_int.out
-  *
-  * \note Even though the returned expression has dynamic size, in the case
-  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
-  * which means that evaluating it does not cause a dynamic memory allocation.
-  *
-  * \sa class Block, block(int,int)
-  */
-template<typename Derived>
-inline typename BlockReturnType<Derived,Dynamic>::SubVectorType
-MatrixBase<Derived>::start(int size)
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,
-               RowsAtCompileTime == 1 ? 1 : Dynamic,
-               ColsAtCompileTime == 1 ? 1 : Dynamic>
-              (derived(), 0, 0,
-               RowsAtCompileTime == 1 ? 1 : size,
-               ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** This is the const version of start(int).*/
-template<typename Derived>
-inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType
-MatrixBase<Derived>::start(int size) const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,
-               RowsAtCompileTime == 1 ? 1 : Dynamic,
-               ColsAtCompileTime == 1 ? 1 : Dynamic>
-              (derived(), 0, 0,
-               RowsAtCompileTime == 1 ? 1 : size,
-               ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** \returns a dynamic-size expression of the last coefficients of *this.
-  *
-  * \only_for_vectors
-  *
-  * \param size the number of coefficients in the block
-  *
-  * \addexample BlockEnd \label How to reference the end of a vector (fixed-size)
-  *
-  * Example: \include MatrixBase_end_int.cpp
-  * Output: \verbinclude MatrixBase_end_int.out
-  *
-  * \note Even though the returned expression has dynamic size, in the case
-  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
-  * which means that evaluating it does not cause a dynamic memory allocation.
-  *
-  * \sa class Block, block(int,int)
-  */
-template<typename Derived>
-inline typename BlockReturnType<Derived,Dynamic>::SubVectorType
-MatrixBase<Derived>::end(int size)
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,
-               RowsAtCompileTime == 1 ? 1 : Dynamic,
-               ColsAtCompileTime == 1 ? 1 : Dynamic>
-              (derived(),
-               RowsAtCompileTime == 1 ? 0 : rows() - size,
-               ColsAtCompileTime == 1 ? 0 : cols() - size,
-               RowsAtCompileTime == 1 ? 1 : size,
-               ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** This is the const version of end(int).*/
-template<typename Derived>
-inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType
-MatrixBase<Derived>::end(int size) const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,
-               RowsAtCompileTime == 1 ? 1 : Dynamic,
-               ColsAtCompileTime == 1 ? 1 : Dynamic>
-              (derived(),
-               RowsAtCompileTime == 1 ? 0 : rows() - size,
-               ColsAtCompileTime == 1 ? 0 : cols() - size,
-               RowsAtCompileTime == 1 ? 1 : size,
-               ColsAtCompileTime == 1 ? 1 : size);
-}
-
-/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this
-  *
-  * \only_for_vectors
-  *
-  * The template parameter \a Size is the number of coefficients in the block
-  *
-  * \param start the index of the first element of the sub-vector
-  *
-  * Example: \include MatrixBase_template_int_segment.cpp
-  * Output: \verbinclude MatrixBase_template_int_segment.out
-  *
-  * \sa class Block
-  */
-template<typename Derived>
-template<int Size>
-inline typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::segment(int start)
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,  (RowsAtCompileTime == 1 ? 1 : Size),
-                         (ColsAtCompileTime == 1 ? 1 : Size)>
-              (derived(), RowsAtCompileTime == 1 ? 0 : start,
-                          ColsAtCompileTime == 1 ? 0 : start);
-}
-
-/** This is the const version of segment<int>(int).*/
-template<typename Derived>
-template<int Size>
-inline const typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::segment(int start) const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived,  (RowsAtCompileTime == 1 ? 1 : Size),
-                         (ColsAtCompileTime == 1 ? 1 : Size)>
-              (derived(), RowsAtCompileTime == 1 ? 0 : start,
-                          ColsAtCompileTime == 1 ? 0 : start);
-}
-
-/** \returns a fixed-size expression of the first coefficients of *this.
-  *
-  * \only_for_vectors
-  *
-  * The template parameter \a Size is the number of coefficients in the block
-  *
-  * \addexample BlockStart \label How to reference the start of a vector (fixed-size)
-  *
-  * Example: \include MatrixBase_template_int_start.cpp
-  * Output: \verbinclude MatrixBase_template_int_start.out
-  *
-  * \sa class Block
-  */
-template<typename Derived>
-template<int Size>
-inline typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::start()
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
-                        (ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0);
-}
-
-/** This is the const version of start<int>().*/
-template<typename Derived>
-template<int Size>
-inline const typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::start() const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
-                        (ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0);
-}
-
-/** \returns a fixed-size expression of the last coefficients of *this.
-  *
-  * \only_for_vectors
-  *
-  * The template parameter \a Size is the number of coefficients in the block
-  *
-  * Example: \include MatrixBase_template_int_end.cpp
-  * Output: \verbinclude MatrixBase_template_int_end.out
-  *
-  * \sa class Block
-  */
-template<typename Derived>
-template<int Size>
-inline typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::end()
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
-                        ColsAtCompileTime == 1 ? 1 : Size>
-           (derived(),
-            RowsAtCompileTime == 1 ? 0 : rows() - Size,
-            ColsAtCompileTime == 1 ? 0 : cols() - Size);
-}
-
-/** This is the const version of end<int>.*/
-template<typename Derived>
-template<int Size>
-inline const typename BlockReturnType<Derived,Size>::SubVectorType
-MatrixBase<Derived>::end() const
-{
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
-                        ColsAtCompileTime == 1 ? 1 : Size>
-           (derived(),
-            RowsAtCompileTime == 1 ? 0 : rows() - Size,
-            ColsAtCompileTime == 1 ? 0 : cols() - Size);
-}
-
 /** \returns a dynamic-size expression of a corner of *this.
   *
   * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight,
diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h
index c8fd98ea1..4785f01f4 100644
--- a/Eigen/src/Core/CwiseUnaryView.h
+++ b/Eigen/src/Core/CwiseUnaryView.h
@@ -42,13 +42,13 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> >
  : ei_traits<MatrixType>
 {
   typedef typename ei_result_of<
-                     ViewOp(typename MatrixType::Scalar)
+                     ViewOp(typename ei_traits<MatrixType>::Scalar)
                    >::type Scalar;
   typedef typename MatrixType::Nested MatrixTypeNested;
-  typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
+  typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested;
   enum {
-    Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)),
-    CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<ViewOp>::Cost
+    Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)),
+    CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits<ViewOp>::Cost
   };
 };
 
@@ -62,7 +62,7 @@ class CwiseUnaryView : ei_no_assignment_operator,
 
     inline CwiseUnaryView(const MatrixType& mat, const ViewOp& func = ViewOp())
       : m_matrix(mat), m_functor(func) {}
-    
+
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryView)
 
     EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); }
@@ -77,7 +77,7 @@ class CwiseUnaryView : ei_no_assignment_operator,
     {
       return m_functor(m_matrix.coeff(index));
     }
-    
+
     EIGEN_STRONG_INLINE Scalar& coeffRef(int row, int col)
     {
       return m_functor(m_matrix.const_cast_derived().coeffRef(row, col));
@@ -89,7 +89,8 @@ class CwiseUnaryView : ei_no_assignment_operator,
     }
 
   protected:
-    const typename MatrixType::Nested m_matrix;
+    // FIXME changed from MatrixType::Nested because of a weird compilation error with sun CC
+    const typename ei_nested<MatrixType>::type m_matrix;
     const ViewOp m_functor;
 };
 
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index 3861984da..4f185ea5b 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -295,7 +295,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<
 /** \returns the \em l2 norm of \c *this using a numerically more stable
   * algorithm.
   *
-  * \sa norm(), dot(), squaredNorm()
+  * \sa norm(), dot(), squaredNorm(), blueNorm()
   */
 template<typename Derived>
 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
@@ -304,6 +304,142 @@ MatrixBase<Derived>::stableNorm() const
   return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
 }
 
+/** \internal Computes ibeta^iexp by binary expansion of iexp,
+  * exact if ibeta is the machine base */
+template<typename T> inline T bexp(int ibeta, int iexp)
+{
+  T tbeta = T(ibeta);
+  T res = 1.0;
+  int n = iexp;
+  if (n<0)
+  {
+    n = - n;
+    tbeta = 1.0/tbeta;
+  }
+  for(;;)
+  {
+    if ((n % 2)==0)
+      res = res * tbeta;
+    n = n/2;
+    if (n==0) return res;
+    tbeta = tbeta*tbeta;
+  }
+  return res;
+}
+
+/** \returns the \em l2 norm of \c *this using the Blue's algorithm.
+  * A Portable Fortran Program to Find the Euclidean Norm of a Vector,
+  * ACM TOMS, Vol 4, Issue 1, 1978.
+  *
+  * \sa norm(), dot(), squaredNorm(), stableNorm()
+  */
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+MatrixBase<Derived>::blueNorm() const
+{
+  static int nmax;
+  static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
+  int n;
+  Scalar ax, abig, amed, asml;
+
+  if(nmax <= 0)
+  {
+    int nbig, ibeta, it, iemin, iemax, iexp;
+    Scalar abig, eps;
+    // This program calculates the machine-dependent constants
+    // bl, b2, slm, s2m, relerr overfl, nmax
+    // from the "basic" machine-dependent numbers
+    // nbig, ibeta, it, iemin, iemax, rbig.
+    // The following define the basic machine-dependent constants.
+    // For portability, the PORT subprograms "ilmaeh" and "rlmach"
+    // are used. For any specific computer, each of the assignment
+    // statements can be replaced
+    nbig  = std::numeric_limits<int>::max();            // largest integer
+    ibeta = NumTraits<Scalar>::Base;                    // base for floating-point numbers
+    it    = NumTraits<Scalar>::Mantissa;                // number of base-beta digits in mantissa
+    iemin = std::numeric_limits<Scalar>::min_exponent;  // minimum exponent
+    iemax = std::numeric_limits<Scalar>::max_exponent;  // maximum exponent
+    rbig  = std::numeric_limits<Scalar>::max();         // largest floating-point number
+
+    // Check the basic machine-dependent constants.
+    if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
+      || (it<=4 && ibeta <= 3 ) || it<2)
+    {
+      ei_assert(false && "the algorithm cannot be guaranteed on this computer");
+    }
+    iexp  = -((1-iemin)/2);
+    b1    = bexp<Scalar>(ibeta, iexp);  // lower boundary of midrange
+    iexp  = (iemax + 1 - it)/2;
+    b2    = bexp<Scalar>(ibeta,iexp);   // upper boundary of midrange
+
+    iexp  = (2-iemin)/2;
+    s1m   = bexp<Scalar>(ibeta,iexp);   // scaling factor for lower range
+    iexp  = - ((iemax+it)/2);
+    s2m   = bexp<Scalar>(ibeta,iexp);   // scaling factor for upper range
+
+    overfl  = rbig*s2m;          // overfow boundary for abig
+    eps     = bexp<Scalar>(ibeta, 1-it);
+    relerr  = ei_sqrt(eps);      // tolerance for neglecting asml
+    abig    = 1.0/eps - 1.0;
+    if (Scalar(nbig)>abig)  nmax = abig;  // largest safe n
+    else                    nmax = nbig;
+  }
+  n = size();
+  if(n==0)
+    return 0;
+  asml = Scalar(0);
+  amed = Scalar(0);
+  abig = Scalar(0);
+  for(int j=0; j<n; ++j)
+  {
+    ax = ei_abs(coeff(j));
+    if(ax > b2)      abig += ei_abs2(ax*s2m);
+    else if(ax < b2) asml += ei_abs2(ax*s1m);
+    else             amed += ei_abs2(ax);
+  }
+  if(abig > Scalar(0))
+  {
+    abig = ei_sqrt(abig);
+    if(abig > overfl)
+    {
+      ei_assert(false && "overflow");
+      return rbig;
+    }
+    if(amed > Scalar(0))
+    {
+      abig = abig/s2m;
+      amed = ei_sqrt(amed);
+    }
+    else
+    {
+      return abig/s2m;
+    }
+
+  }
+  else if(asml > Scalar(0))
+  {
+    if (amed > Scalar(0))
+    {
+      abig = ei_sqrt(amed);
+      amed = ei_sqrt(asml) / s1m;
+    }
+    else
+    {
+      return ei_sqrt(asml)/s1m;
+    }
+  }
+  else
+  {
+    return ei_sqrt(amed);
+  }
+  asml = std::min(abig, amed);
+  abig = std::max(abig, amed);
+  if(asml <= abig*relerr)
+    return abig;
+  else
+    return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig));
+}
+
 /** \returns an expression of the quotient of *this by its own norm.
   *
   * \only_for_vectors
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 8cb1f7b40..3527042f9 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -151,7 +151,7 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float
 *** double  ***
 **************/
 
-template<> inline double precision<double>() { return 1e-11; }
+template<> inline double precision<double>() { return 1e-12; }
 template<> inline double machine_epsilon<double>() { return 2.220e-16; }
 
 inline double ei_real(double x)  { return x; }
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 18c8f969a..2b4c4634a 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -124,7 +124,7 @@ class Matrix
 {
   public:
     EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix)
-    
+
     enum { Options = _Options };
     friend class Eigen::Map<Matrix, Unaligned>;
     typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType;
@@ -218,7 +218,7 @@ class Matrix
       *
       * This method is intended for dynamic-size matrices, although it is legal to call it on any
       * matrix as long as fixed dimensions are left unchanged. If you only want to change the number
-      * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t). 
+      * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t).
       *
       * If the current number of coefficients of \c *this exactly matches the
       * product \a rows * \a cols, then no memory allocation is performed and
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index cfa41978e..0a3342758 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -163,10 +163,14 @@ template<typename Derived> class MatrixBase
           * constructed from this one. See the \ref flags "list of flags".
           */
 
-      CoeffReadCost = ei_traits<Derived>::CoeffReadCost
+      CoeffReadCost = ei_traits<Derived>::CoeffReadCost,
         /**< This is a rough measure of how expensive it is to read one coefficient from
           * this expression.
           */
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+      _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC
+#endif
     };
 
     /** Default constructor. Just checks at compile-time for self-consistency of the flags. */
@@ -230,7 +234,7 @@ template<typename Derived> class MatrixBase
 
     /** \internal the return type of coeff()
       */
-    typedef typename ei_meta_if<bool(int(Flags)&DirectAccessBit), const Scalar&, Scalar>::ret CoeffReturnType;
+    typedef typename ei_meta_if<_HasDirectAccess, const Scalar&, Scalar>::ret CoeffReturnType;
 
     /** \internal Represents a matrix with all coefficients equal to one another*/
     typedef CwiseNullaryOp<ei_scalar_constant_op<Scalar>,Derived> ConstantReturnType;
@@ -426,8 +430,9 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     Scalar dot(const MatrixBase<OtherDerived>& other) const;
     RealScalar squaredNorm() const;
-    RealScalar norm()  const;
-    RealScalar stableNorm()  const;
+    RealScalar norm() const;
+    RealScalar stableNorm() const;
+    RealScalar blueNorm() const;
     const PlainMatrixType normalized() const;
     void normalize();
 
@@ -450,14 +455,14 @@ template<typename Derived> class MatrixBase
     const typename BlockReturnType<Derived>::Type
     block(int startRow, int startCol, int blockRows, int blockCols) const;
 
-    typename BlockReturnType<Derived>::SubVectorType segment(int start, int size);
-    const typename BlockReturnType<Derived>::SubVectorType segment(int start, int size) const;
+    VectorBlock<Derived> segment(int start, int size);
+    const VectorBlock<Derived> segment(int start, int size) const;
 
-    typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size);
-    const typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size) const;
+    VectorBlock<Derived> start(int size);
+    const VectorBlock<Derived> start(int size) const;
 
-    typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size);
-    const typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size) const;
+    VectorBlock<Derived> end(int size);
+    const VectorBlock<Derived> end(int size) const;
 
     typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols);
     const typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols) const;
@@ -472,14 +477,14 @@ template<typename Derived> class MatrixBase
     template<int CRows, int CCols>
     const typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type) const;
 
-    template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType start(void);
-    template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType start() const;
+    template<int Size> VectorBlock<Derived,Size> start(void);
+    template<int Size> const VectorBlock<Derived,Size> start() const;
 
-    template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType end();
-    template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType end() const;
+    template<int Size> VectorBlock<Derived,Size> end();
+    template<int Size> const VectorBlock<Derived,Size> end() const;
 
-    template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType segment(int start);
-    template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType segment(int start) const;
+    template<int Size> VectorBlock<Derived,Size> segment(int start);
+    template<int Size> const VectorBlock<Derived,Size> segment(int start) const;
 
     Diagonal<Derived,0> diagonal();
     const Diagonal<Derived,0> diagonal() const;
@@ -696,6 +701,7 @@ template<typename Derived> class MatrixBase
     const PartialLU<PlainMatrixType> partialLu() const;
     const PlainMatrixType inverse() const;
     void computeInverse(PlainMatrixType *result) const;
+    bool computeInverseWithCheck( PlainMatrixType *result ) const;
     Scalar determinant() const;
 
 /////////// Cholesky module ///////////
@@ -705,7 +711,7 @@ template<typename Derived> class MatrixBase
 
 /////////// QR module ///////////
 
-    const QR<PlainMatrixType> qr() const;
+    const HouseholderQR<PlainMatrixType> householderQr() const;
 
     EigenvaluesReturnType eigenvalues() const;
     RealScalar operatorNorm() const;
diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h
index 24afe54c5..dec07a036 100644
--- a/Eigen/src/Core/NumTraits.h
+++ b/Eigen/src/Core/NumTraits.h
@@ -70,7 +70,9 @@ template<> struct NumTraits<float>
     HasFloatingPoint = 1,
     ReadCost = 1,
     AddCost = 1,
-    MulCost = 1
+    MulCost = 1,
+    Base = 2,
+    Mantissa = 23
   };
 };
 
@@ -83,7 +85,9 @@ template<> struct NumTraits<double>
     HasFloatingPoint = 1,
     ReadCost = 1,
     AddCost = 1,
-    MulCost = 1
+    MulCost = 1,
+    Base = 2,
+    Mantissa = 52
   };
 };
 
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 44e3f606e..f0c4480de 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -155,13 +155,14 @@ template<typename Lhs, typename Rhs> struct ei_product_mode
   typedef typename ei_blas_traits<Lhs>::ActualXprType ActualLhs;
   typedef typename ei_blas_traits<Rhs>::ActualXprType ActualRhs;
   enum{
-
-    value = Lhs::MaxColsAtCompileTime == Dynamic
-            && ( Lhs::MaxRowsAtCompileTime == Dynamic
-              || Rhs::MaxColsAtCompileTime == Dynamic )
-            && (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit)  && (!(ActualLhs::Flags&DirectAccessBit))))
-            && (!(Lhs::IsVectorAtCompileTime && (!(Rhs::Flags&RowMajorBit)) && (!(ActualRhs::Flags&DirectAccessBit))))
-            && (ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)
+    // workaround sun studio:
+    LhsIsVectorAtCompileTime = ei_traits<Lhs>::ColsAtCompileTime==1 || ei_traits<Rhs>::ColsAtCompileTime==1,
+    value = ei_traits<Lhs>::MaxColsAtCompileTime == Dynamic
+            && ( ei_traits<Lhs>::MaxRowsAtCompileTime == Dynamic
+              || ei_traits<Rhs>::MaxColsAtCompileTime == Dynamic )
+            && (!(Rhs::IsVectorAtCompileTime && (ei_traits<Lhs>::Flags&RowMajorBit)  && (!(ei_traits<Lhs>::Flags&DirectAccessBit))))
+            && (!(LhsIsVectorAtCompileTime && (!(ei_traits<Rhs>::Flags&RowMajorBit)) && (!(ei_traits<Rhs>::Flags&DirectAccessBit))))
+            && (ei_is_same_type<typename ei_traits<Lhs>::Scalar, typename ei_traits<Rhs>::Scalar>::ret)
           ? CacheFriendlyProduct
           : NormalProduct };
 };
@@ -577,7 +578,7 @@ struct ei_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMod
 // Forward declarations
 
 template<typename Scalar, bool ConjugateLhs, bool ConjugateRhs>
-void ei_cache_friendly_product(
+static void ei_cache_friendly_product(
   int _rows, int _cols, int depth,
   bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
   bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index 76e4289de..cb162ca91 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h
@@ -180,7 +180,7 @@ struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
 
 template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
 struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> {
-  static void run(const Lhs& lhs, Rhs& rhs) {}
+  static void run(const Lhs&, Rhs&) {}
 };
 
 template<typename Lhs, typename Rhs, int Mode, int StorageOrder>
diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h
new file mode 100644
index 000000000..7ce5977f6
--- /dev/null
+++ b/Eigen/src/Core/VectorBlock.h
@@ -0,0 +1,311 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2006-2008 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_VECTORBLOCK_H
+#define EIGEN_VECTORBLOCK_H
+
+/** \class VectorBlock
+  *
+  * \brief Expression of a fixed-size or dynamic-size sub-vector
+  *
+  * \param VectorType the type of the object in which we are taking a sub-vector
+  * \param Size size of the sub-vector we are taking at compile time (optional)
+  * \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned.
+  *                      The default is AsRequested. This parameter is internaly used by Eigen
+  *                      in expressions such as \code mat.segment() += other; \endcode and most of
+  *                      the time this is the only way it is used.
+  *
+  * This class represents an expression of either a fixed-size or dynamic-size sub-vector.
+  * It is the return type of MatrixBase::segment(int,int) and MatrixBase::segment<int>(int) and
+  * most of the time this is the only way it is used.
+  *
+  * However, if you want to directly maniputate sub-vector expressions,
+  * for instance if you want to write a function returning such an expression, you
+  * will need to use this class.
+  *
+  * Here is an example illustrating the dynamic case:
+  * \include class_VectorBlock.cpp
+  * Output: \verbinclude class_VectorBlock.out
+  *
+  * \note Even though this expression has dynamic size, in the case where \a VectorType
+  * has fixed size, this expression inherits a fixed maximal size which means that evaluating
+  * it does not cause a dynamic memory allocation.
+  *
+  * Here is an example illustrating the fixed-size case:
+  * \include class_FixedVectorBlock.cpp
+  * Output: \verbinclude class_FixedVectorBlock.out
+  *
+  * \sa class Block, MatrixBase::segment(int,int,int,int), MatrixBase::segment(int,int)
+  */
+template<typename VectorType, int Size, int _PacketAccess>
+struct ei_traits<VectorBlock<VectorType, Size, _PacketAccess> >
+  : public ei_traits<Block<VectorType,
+                 ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
+                 ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
+                 _PacketAccess> >
+{
+};
+
+template<typename VectorType, int Size, int PacketAccess> class VectorBlock
+  : public Block<VectorType,
+                 ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
+                 ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
+                 PacketAccess>
+{
+    typedef Block<VectorType,
+                  ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
+                  ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
+                  PacketAccess> Base;
+    enum {
+      IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1
+    };
+  public:
+
+    using Base::operator=;
+    using Base::operator+=;
+    using Base::operator-=;
+    using Base::operator*=;
+    using Base::operator/=;
+    
+    /** Dynamic-size constructor
+      */
+    inline VectorBlock(const VectorType& vector, int start, int size)
+      : Base(vector,
+             IsColVector ? start : 0, IsColVector ? 0 : start,
+             IsColVector ? size  : 1, IsColVector ? 1 : size)
+    {
+      
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
+    }
+
+    /** Fixed-size constructor
+      */
+    inline VectorBlock(const VectorType& vector, int start)
+      : Base(vector, IsColVector ? start : 0, IsColVector ? 0 : start)
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
+    }
+};
+
+
+/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this.
+  *
+  * \only_for_vectors
+  *
+  * \addexample VectorBlockIntInt \label How to reference a sub-vector (dynamic size)
+  *
+  * \param start the first coefficient in the segment
+  * \param size the number of coefficients in the segment
+  *
+  * Example: \include MatrixBase_segment_int_int.cpp
+  * Output: \verbinclude MatrixBase_segment_int_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, segment(int)
+  */
+template<typename Derived>
+inline VectorBlock<Derived> MatrixBase<Derived>
+  ::segment(int start, int size)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), start, size);
+}
+
+/** This is the const version of segment(int,int).*/
+template<typename Derived>
+inline const VectorBlock<Derived>
+MatrixBase<Derived>::segment(int start, int size) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), start, size);
+}
+
+/** \returns a dynamic-size expression of the first coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * \param size the number of coefficients in the block
+  *
+  * \addexample BlockInt \label How to reference a sub-vector (fixed-size)
+  *
+  * Example: \include MatrixBase_start_int.cpp
+  * Output: \verbinclude MatrixBase_start_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, block(int,int)
+  */
+template<typename Derived>
+inline VectorBlock<Derived>
+MatrixBase<Derived>::start(int size)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), 0, size);
+}
+
+/** This is the const version of start(int).*/
+template<typename Derived>
+inline const VectorBlock<Derived>
+MatrixBase<Derived>::start(int size) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), 0, size);
+}
+
+/** \returns a dynamic-size expression of the last coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * \param size the number of coefficients in the block
+  *
+  * \addexample BlockEnd \label How to reference the end of a vector (fixed-size)
+  *
+  * Example: \include MatrixBase_end_int.cpp
+  * Output: \verbinclude MatrixBase_end_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, block(int,int)
+  */
+template<typename Derived>
+inline VectorBlock<Derived>
+MatrixBase<Derived>::end(int size)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), this->size() - size, size);
+}
+
+/** This is the const version of end(int).*/
+template<typename Derived>
+inline const VectorBlock<Derived>
+MatrixBase<Derived>::end(int size) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived>(derived(), this->size() - size, size);
+}
+
+/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * \param start the index of the first element of the sub-vector
+  *
+  * Example: \include MatrixBase_template_int_segment.cpp
+  * Output: \verbinclude MatrixBase_template_int_segment.out
+  *
+  * \sa class Block
+  */
+template<typename Derived>
+template<int Size>
+inline VectorBlock<Derived,Size>
+MatrixBase<Derived>::segment(int start)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived,Size>(derived(), start);
+}
+
+/** This is the const version of segment<int>(int).*/
+template<typename Derived>
+template<int Size>
+inline const VectorBlock<Derived,Size>
+MatrixBase<Derived>::segment(int start) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived,Size>(derived(), start);
+}
+
+/** \returns a fixed-size expression of the first coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * \addexample BlockStart \label How to reference the start of a vector (fixed-size)
+  *
+  * Example: \include MatrixBase_template_int_start.cpp
+  * Output: \verbinclude MatrixBase_template_int_start.out
+  *
+  * \sa class Block
+  */
+template<typename Derived>
+template<int Size>
+inline VectorBlock<Derived,Size>
+MatrixBase<Derived>::start()
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived,Size>(derived(), 0);
+}
+
+/** This is the const version of start<int>().*/
+template<typename Derived>
+template<int Size>
+inline const VectorBlock<Derived,Size>
+MatrixBase<Derived>::start() const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived,Size>(derived(), 0);
+}
+
+/** \returns a fixed-size expression of the last coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_template_int_end.cpp
+  * Output: \verbinclude MatrixBase_template_int_end.out
+  *
+  * \sa class Block
+  */
+template<typename Derived>
+template<int Size>
+inline VectorBlock<Derived,Size>
+MatrixBase<Derived>::end()
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived, Size>(derived(), size() - Size);
+}
+
+/** This is the const version of end<int>.*/
+template<typename Derived>
+template<int Size>
+inline const VectorBlock<Derived,Size>
+MatrixBase<Derived>::end() const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return VectorBlock<Derived, Size>(derived(), size() - Size);
+}
+
+
+#endif // EIGEN_VECTORBLOCK_H
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 454b44e52..0c251022b 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -185,7 +185,7 @@ const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
                                   | EvalBeforeAssigningBit
                                   | SparseBit;
-    
+
 // Possible values for the Mode parameter of part()
 const unsigned int UpperTriangular = UpperTriangularBit;
 const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit;
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 3e752c84b..363972b60 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -42,6 +42,7 @@ template<typename MatrixType> class Minor;
 template<typename MatrixType, int BlockRows=Dynamic, int BlockCols=Dynamic, int PacketAccess=AsRequested,
          int _DirectAccessStatus = ei_traits<MatrixType>::Flags&DirectAccessBit ? DirectAccessBit
                                  : ei_traits<MatrixType>::Flags&SparseBit> class Block;
+template<typename MatrixType, int Size=Dynamic, int PacketAccess=AsRequested> class VectorBlock;
 template<typename MatrixType> class Transpose;
 template<typename MatrixType> class Conjugate;
 template<typename NullaryOp, typename MatrixType>         class CwiseNullaryOp;
@@ -111,7 +112,7 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse;
 
 template<typename MatrixType> class LU;
 template<typename MatrixType> class PartialLU;
-template<typename MatrixType> class QR;
+template<typename MatrixType> class HouseholderQR;
 template<typename MatrixType> class SVD;
 template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
 template<typename MatrixType> class LDLT;
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 0cc5ae9aa..c9f2f4d40 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -51,7 +51,8 @@
 #define EIGEN_GCC3_OR_OLDER 0
 #endif
 
-#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER
+// FIXME vectorization + alignment is completely disabled with sun studio
+#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC)
   #define EIGEN_ARCH_WANTS_ALIGNMENT 1
 #else
   #define EIGEN_ARCH_WANTS_ALIGNMENT 0
@@ -104,7 +105,7 @@
 /** Allows to disable some optimizations which might affect the accuracy of the result.
   * Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them.
   * They currently include:
-  *   - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled. 
+  *   - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled.
   */
 #ifndef EIGEN_FAST_MATH
 #define EIGEN_FAST_MATH 1
@@ -206,13 +207,16 @@ using Eigen::ei_cos;
  * vectorized and non-vectorized code.
  */
 #if !EIGEN_ALIGN
-#define EIGEN_ALIGN_128
+  #define EIGEN_ALIGN_128
 #elif (defined __GNUC__)
-#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
+  #define EIGEN_ALIGN_128 __attribute__((aligned(16)))
 #elif (defined _MSC_VER)
-#define EIGEN_ALIGN_128 __declspec(align(16))
+  #define EIGEN_ALIGN_128 __declspec(align(16))
+#elif (defined __SUNPRO_CC)
+  // FIXME not sure about this one:
+  #define EIGEN_ALIGN_128 __attribute__((aligned(16)))
 #else
-#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
+  #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
 #endif
 
 #define EIGEN_RESTRICT __restrict
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index cc3aa4fac..f8581eebc 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -27,7 +27,17 @@
 #ifndef EIGEN_MEMORY_H
 #define EIGEN_MEMORY_H
 
-#if defined(__APPLE__) || defined(_WIN64) || defined (__FreeBSD__)
+// FreeBSD 6 seems to have 16-byte aligned malloc
+// See http://svn.freebsd.org/viewvc/base/stable/6/lib/libc/stdlib/malloc.c?view=markup
+// FreeBSD 7 seems to have 16-byte aligned malloc except on ARM and MIPS architectures
+// See http://svn.freebsd.org/viewvc/base/stable/7/lib/libc/stdlib/malloc.c?view=markup
+#if defined(__FreeBSD__) && !defined(__arm__) && !defined(__mips__)
+#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 1
+#else
+#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 0
+#endif
+
+#if defined(__APPLE__) || defined(_WIN64) || EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED
   #define EIGEN_MALLOC_ALREADY_ALIGNED 1
 #else
   #define EIGEN_MALLOC_ALREADY_ALIGNED 0
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index 7b139b0c1..a4ffb8a20 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -217,7 +217,7 @@ template<typename Derived,typename Scalar,typename OtherScalar,
          bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret >
 struct ei_special_scalar_op_base
 {
-  // dummy operator* so that the 
+  // dummy operator* so that the
   // "using ei_special_scalar_op_base::operator*" compiles
   void operator*() const;
 };
@@ -237,8 +237,6 @@ struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>
   * TODO: could be a good idea to define a big ReturnType struct ??
   */
 template<typename ExpressionType, int RowsOrSize=Dynamic, int Cols=Dynamic> struct BlockReturnType {
-  typedef Block<ExpressionType, (ei_traits<ExpressionType>::RowsAtCompileTime == 1 ? 1 : RowsOrSize),
-                                (ei_traits<ExpressionType>::ColsAtCompileTime == 1 ? 1 : RowsOrSize)> SubVectorType;
   typedef Block<ExpressionType, RowsOrSize, Cols> Type;
 };
 
@@ -251,7 +249,7 @@ template<typename ExpressionType> struct HNormalizedReturnType {
   typedef Block<ExpressionType,
                 ei_traits<ExpressionType>::ColsAtCompileTime==1 ? SizeMinusOne : 1,
                 ei_traits<ExpressionType>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne;
-  typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>, 
+  typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>,
               NestByValue<StartMinusOne> > Type;
 };
 
diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h
index 66580cf1b..9385f259d 100644
--- a/Eigen/src/Geometry/Quaternion.h
+++ b/Eigen/src/Geometry/Quaternion.h
@@ -303,7 +303,9 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const MatrixBase<Derive
   return *this;
 }
 
-/** Convert the quaternion to a 3x3 rotation matrix */
+/** Convert the quaternion to a 3x3 rotation matrix. The quaternion is required to
+  * be normalized, otherwise the result is undefined.
+  */
 template<typename Scalar>
 inline typename Quaternion<Scalar>::Matrix3
 Quaternion<Scalar>::toRotationMatrix(void) const
@@ -340,11 +342,15 @@ Quaternion<Scalar>::toRotationMatrix(void) const
   return res;
 }
 
-/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b.
+/** Sets \c *this to be a quaternion representing a rotation between
+  * the two arbitrary vectors \a a and \a b. In other words, the built
+  * rotation represent a rotation sending the line of direction \a a
+  * to the line of direction \a b, both lines passing through the origin.
   *
-  * \returns a reference to *this.
+  * \returns a reference to \c *this.
   *
-  * Note that the two input vectors do \b not have to be normalized.
+  * Note that the two input vectors do \b not have to be normalized, and
+  * do not need to have the same norm.
   */
 template<typename Scalar>
 template<typename Derived1, typename Derived2>
@@ -354,21 +360,26 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
   Vector3 v1 = b.normalized();
   Scalar c = v0.dot(v1);
 
-  // if dot == 1, vectors are the same
-  if (ei_isApprox(c,Scalar(1)))
+  // if dot == -1, vectors are nearly opposites
+  // => accuraletly compute the rotation axis by computing the
+  //    intersection of the two planes. This is done by solving:
+  //       x^T v0 = 0
+  //       x^T v1 = 0
+  //    under the constraint:
+  //       ||x|| = 1
+  //    which yields a singular value problem
+  if (c < Scalar(-1)+precision<Scalar>())
   {
-    // set to identity
-    this->w() = 1; this->vec().setZero();
+    c = std::max<Scalar>(c,-1);
+    Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
+    SVD<Matrix<Scalar,2,3> > svd(m);
+    Vector3 axis = svd.matrixV().col(2);
+
+    Scalar w2 = (Scalar(1)+c)*Scalar(0.5);
+    this->w() = ei_sqrt(w2);
+    this->vec() = axis * ei_sqrt(Scalar(1) - w2);
     return *this;
   }
-  // if dot == -1, vectors are opposites
-  if (ei_isApprox(c,Scalar(-1)))
-  {
-    this->vec() = v0.unitOrthogonal();
-    this->w() = 0;
-    return *this;
-  }
-
   Vector3 axis = v0.cross(v1);
   Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2));
   Scalar invs = Scalar(1)/s;
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index b29efb60c..5af14813d 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -29,41 +29,45 @@
 *** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
 ********************************************************************/
 
-template<typename MatrixType>
-void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
+template<typename XprType, typename MatrixType>
+inline void ei_compute_inverse_size2_helper(
+    const XprType& matrix, const typename MatrixType::Scalar& invdet,
+    MatrixType* result)
 {
-  typedef typename MatrixType::Scalar Scalar;
-  const Scalar invdet = Scalar(1) / matrix.determinant();
   result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
   result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
   result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
   result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
 }
 
+template<typename MatrixType>
+inline void ei_compute_inverse_size2(const MatrixType& matrix, MatrixType* result)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  const Scalar invdet = Scalar(1) / matrix.determinant();
+  ei_compute_inverse_size2_helper( matrix, invdet, result );
+}
+
 template<typename XprType, typename MatrixType>
-bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result)
+bool ei_compute_inverse_size2_with_check(const XprType& matrix, MatrixType* result)
 {
   typedef typename MatrixType::Scalar Scalar;
   const Scalar det = matrix.determinant();
   if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
   const Scalar invdet = Scalar(1) / det;
-  result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
-  result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
-  result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
-  result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
+  ei_compute_inverse_size2_helper( matrix, invdet, result );
   return true;
 }
 
-template<typename MatrixType>
-void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
+template<typename XprType, typename MatrixType>
+void ei_compute_inverse_size3_helper(
+    const XprType& matrix,
+    const typename MatrixType::Scalar& invdet,
+    const typename MatrixType::Scalar& det_minor00,
+    const typename MatrixType::Scalar& det_minor10,
+    const typename MatrixType::Scalar& det_minor20,
+    MatrixType* result)
 {
-  typedef typename MatrixType::Scalar Scalar;
-  const Scalar det_minor00 = matrix.minor(0,0).determinant();
-  const Scalar det_minor10 = matrix.minor(1,0).determinant();
-  const Scalar det_minor20 = matrix.minor(2,0).determinant();
-  const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0)
-                                    - det_minor10 * matrix.coeff(1,0)
-                                    + det_minor20 * matrix.coeff(2,0) );
   result->coeffRef(0, 0) = det_minor00 * invdet;
   result->coeffRef(0, 1) = -det_minor10 * invdet;
   result->coeffRef(0, 2) = det_minor20 * invdet;
@@ -75,8 +79,24 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu
   result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
 }
 
+template<bool Check, typename XprType, typename MatrixType>
+bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  const Scalar det_minor00 = matrix.minor(0,0).determinant();
+  const Scalar det_minor10 = matrix.minor(1,0).determinant();
+  const Scalar det_minor20 = matrix.minor(2,0).determinant();
+  const Scalar det = ( det_minor00 * matrix.coeff(0,0)
+      - det_minor10 * matrix.coeff(1,0)
+      + det_minor20 * matrix.coeff(2,0) );
+  if(Check) if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
+  const Scalar invdet = Scalar(1) / det;
+  ei_compute_inverse_size3_helper( matrix, invdet, det_minor00, det_minor10, det_minor20, result );
+  return true;
+}
+
 template<typename MatrixType>
-bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
+bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result)
 {
   /* Let's split M into four 2x2 blocks:
     * (P Q)
@@ -94,7 +114,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
   typedef Block<MatrixType,2,2> XprBlock22;
   typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
   Block22 P_inverse;
-  if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
+  if(ei_compute_inverse_size2_with_check(matrix.template block<2,2>(0,0), &P_inverse))
   {
     const Block22 Q = matrix.template block<2,2>(0,2);
     const Block22 P_inverse_times_Q = P_inverse * Q;
@@ -104,7 +124,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
     const XprBlock22 S = matrix.template block<2,2>(2,2);
     const Block22 X = S - R_times_P_inverse_times_Q;
     Block22 Y;
-    ei_compute_inverse_in_size2_case(X, &Y);
+    ei_compute_inverse_size2(X, &Y);
     result->template block<2,2>(2,2) = Y;
     result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
     const Block22 Z = P_inverse_times_Q * Y;
@@ -118,13 +138,13 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
   }
 }
 
-template<typename MatrixType>
-void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
+template<typename XprType, typename MatrixType>
+bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* result)
 {
-  if(ei_compute_inverse_in_size4_case_helper(matrix, result))
+  if(ei_compute_inverse_size4_helper(matrix, result))
   {
     // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
-    return;
+    return true;
   }
   else
   {
@@ -134,17 +154,17 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
     MatrixType m(matrix);
     m.row(0).swap(m.row(2));
     m.row(1).swap(m.row(3));
-    if(ei_compute_inverse_in_size4_case_helper(m, result))
+    if(ei_compute_inverse_size4_helper(m, result))
     {
       // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
       // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
       // the corresponding columns.
       result->col(0).swap(result->col(2));
       result->col(1).swap(result->col(3));
+    return true;
     }
     else
     {
-      // last possible case. Since matrix is assumed to be invertible, this last case has to work.
       // first, undo the swaps previously made
       m.row(0).swap(m.row(2));
       m.row(1).swap(m.row(3));
@@ -154,13 +174,23 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
       // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
       int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
       m.row(1).swap(m.row(swap1with));
-      ei_compute_inverse_in_size4_case_helper(m, result);
-      result->col(1).swap(result->col(swap1with));
-      result->col(0).swap(result->col(swap0with));
+      if( ei_compute_inverse_size4_helper(m, result) )
+      {
+        result->col(1).swap(result->col(swap1with));
+        result->col(0).swap(result->col(swap0with));
+        return true;
+      }
+      else
+      {
+        // non-invertible matrix
+        return false;
+      }
     }
   }
 }
 
+
+
 /***********************************************
 *** Part 2 : selector and MatrixBase methods ***
 ***********************************************/
@@ -189,7 +219,7 @@ struct ei_compute_inverse<MatrixType, 2>
 {
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
-    ei_compute_inverse_in_size2_case(matrix, result);
+    ei_compute_inverse_size2(matrix, result);
   }
 };
 
@@ -198,7 +228,7 @@ struct ei_compute_inverse<MatrixType, 3>
 {
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
-    ei_compute_inverse_in_size3_case(matrix, result);
+    ei_compute_inverse_size3<false, MatrixType, MatrixType>(matrix, result);
   }
 };
 
@@ -207,7 +237,7 @@ struct ei_compute_inverse<MatrixType, 4>
 {
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
-    ei_compute_inverse_in_size4_case(matrix, result);
+    ei_compute_inverse_size4_with_check(matrix, result);
   }
 };
 
@@ -215,14 +245,15 @@ struct ei_compute_inverse<MatrixType, 4>
   *
   * Computes the matrix inverse of this matrix.
   *
-  * \note This matrix must be invertible, otherwise the result is undefined.
+  * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
+  * computeInverseWithCheck().
   *
   * \param result Pointer to the matrix in which to store the result.
   *
   * Example: \include MatrixBase_computeInverse.cpp
   * Output: \verbinclude MatrixBase_computeInverse.out
   *
-  * \sa inverse()
+  * \sa inverse(), computeInverseWithCheck()
   */
 template<typename Derived>
 inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
@@ -236,7 +267,8 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
   *
   * \returns the matrix inverse of this matrix.
   *
-  * \note This matrix must be invertible, otherwise the result is undefined.
+  * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
+  * computeInverseWithCheck().
   *
   * \note This method returns a matrix by value, which can be inefficient. To avoid that overhead,
   * use computeInverse() instead.
@@ -244,7 +276,7 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
   * Example: \include MatrixBase_inverse.cpp
   * Output: \verbinclude MatrixBase_inverse.out
   *
-  * \sa computeInverse()
+  * \sa computeInverse(), computeInverseWithCheck()
   */
 template<typename Derived>
 inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const
@@ -254,4 +286,81 @@ inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::
   return result;
 }
 
+
+/********************************************
+ * Compute inverse with invertibility check *
+ *******************************************/
+
+template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    LU<MatrixType> lu( matrix );
+    if( !lu.isInvertible() ) return false;
+    lu.computeInverse(result);
+    return true;
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse_with_check<MatrixType, 1>
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    if( 0 == result->coeffRef(0,0) ) return false;
+  
+    typedef typename MatrixType::Scalar Scalar;
+    result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
+    return true;
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse_with_check<MatrixType, 2>
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_size2_with_check(matrix, result);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse_with_check<MatrixType, 3>
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_size3<true, MatrixType, MatrixType>(matrix, result);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse_with_check<MatrixType, 4>
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_size4_with_check(matrix, result);
+  }
+};
+
+/** \lu_module
+  *
+  * Computation of matrix inverse, with invertibility check.
+  *
+  * \returns true if the matrix is invertible, false otherwise.
+  *
+  * \param result Pointer to the matrix in which to store the result.
+  *
+  * \sa inverse(), computeInverse()
+  */
+template<typename Derived>
+inline bool MatrixBase<Derived>::computeInverseWithCheck(PlainMatrixType *result) const
+{
+  ei_assert(rows() == cols());
+  EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
+  return ei_compute_inverse_with_check<PlainMatrixType>::run(eval(), result);
+}
+
+
 #endif // EIGEN_INVERSE_H
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 5883eed65..a23a3a035 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -28,18 +28,20 @@
 /** \ingroup QR_Module
   * \nonstableyet
   *
-  * \class QR
+  * \class HouseholderQR
   *
-  * \brief QR decomposition of a matrix
+  * \brief Householder QR decomposition of a matrix
   *
   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
   *
   * This class performs a QR decomposition using Householder transformations. The result is
   * stored in a compact way compatible with LAPACK.
   *
+  * Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
+  *
   * \sa MatrixBase::qr()
   */
-template<typename MatrixType> class QR
+template<typename MatrixType> class HouseholderQR
 {
   public:
 
@@ -53,88 +55,23 @@ template<typename MatrixType> class QR
     * \brief Default Constructor.
     *
     * The default constructor is useful in cases in which the user intends to
-    * perform decompositions via QR::compute(const MatrixType&).
+    * perform decompositions via HouseholderQR::compute(const MatrixType&).
     */
-    QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+    HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
 
-    QR(const MatrixType& matrix)
+    HouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs(matrix.cols()),
         m_isInitialized(false)
     {
       compute(matrix);
     }
-    
-    /** \deprecated use isInjective()
-      * \returns whether or not the matrix is of full rank
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    EIGEN_DEPRECATED bool isFullRank() const 
-    { 
-      ei_assert(m_isInitialized && "QR is not initialized.");
-      return rank() == m_qr.cols(); 
-    }
-    
-    /** \returns the rank of the matrix of which *this is the QR decomposition.
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    int rank() const;
-    
-    /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    inline int dimensionOfKernel() const
-    {
-      ei_assert(m_isInitialized && "QR is not initialized.");
-      return m_qr.cols() - rank();
-    }
-    
-    /** \returns true if the matrix of which *this is the QR decomposition represents an injective
-      *          linear map, i.e. has trivial kernel; false otherwise.
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    inline bool isInjective() const
-    {
-      ei_assert(m_isInitialized && "QR is not initialized.");
-      return rank() == m_qr.cols();
-    }
-    
-    /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
-      *          linear map; false otherwise.
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    inline bool isSurjective() const
-    {
-      ei_assert(m_isInitialized && "QR is not initialized.");
-      return rank() == m_qr.rows();
-    }
-
-    /** \returns true if the matrix of which *this is the QR decomposition is invertible.
-      *
-      * \note Since the rank is computed only once, i.e. the first time it is needed, this
-      *       method almost does not perform any further computation.
-      */
-    inline bool isInvertible() const
-    {
-      ei_assert(m_isInitialized && "QR is not initialized.");
-      return isInjective() && isSurjective();
-    }
-    
+        
     /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
     const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular>
     matrixR(void) const
     {
-      ei_assert(m_isInitialized && "QR is not initialized.");
+      ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
       int cols = m_qr.cols();
       return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
     }
@@ -148,58 +85,35 @@ template<typename MatrixType> class QR
       *          Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
       *          If no solution exists, *result is left with undefined coefficients.
       *
-      * \returns true if any solution exists, false if no solution exists.
-      *
-      * \note If there exist more than one solution, this method will arbitrarily choose one.
-      *       If you need a complete analysis of the space of solutions, take the one solution obtained
-      *       by this method and add to it elements of the kernel, as determined by kernel().
-      *
       * \note The case where b is a matrix is not yet implemented. Also, this
       *       code is space inefficient.
       *
-      * Example: \include QR_solve.cpp
-      * Output: \verbinclude QR_solve.out
-      *
-      * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
+      * Example: \include HouseholderQR_solve.cpp
+      * Output: \verbinclude HouseholderQR_solve.out
       */
     template<typename OtherDerived, typename ResultType>
-    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+    void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
     MatrixType matrixQ(void) const;
+    
+    /** \returns a reference to the matrix where the Householder QR decomposition is stored
+      * in a LAPACK-compatible way.
+      */
+    const MatrixType& matrixQR() const { return m_qr; }
 
     void compute(const MatrixType& matrix);
 
   protected:
     MatrixType m_qr;
     VectorType m_hCoeffs;
-    mutable int m_rank;
-    mutable bool m_rankIsUptodate;
     bool m_isInitialized;
 };
 
-/** \returns the rank of the matrix of which *this is the QR decomposition. */
-template<typename MatrixType>
-int QR<MatrixType>::rank() const
-{
-  ei_assert(m_isInitialized && "QR is not initialized.");
-  if (!m_rankIsUptodate)
-  {
-    RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
-    int n = m_qr.cols();
-    m_rank = 0;
-    while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
-      ++m_rank;
-    m_rankIsUptodate = true;
-  }
-  return m_rank;
-}
-
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
 template<typename MatrixType>
-void QR<MatrixType>::compute(const MatrixType& matrix)
+void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 { 
-  m_rankIsUptodate = false;
   m_qr = matrix;
   m_hCoeffs.resize(matrix.cols());
 
@@ -262,12 +176,12 @@ void QR<MatrixType>::compute(const MatrixType& matrix)
 
 template<typename MatrixType>
 template<typename OtherDerived, typename ResultType>
-bool QR<MatrixType>::solve(
+void HouseholderQR<MatrixType>::solve(
   const MatrixBase<OtherDerived>& b,
   ResultType *result
 ) const
 {
-  ei_assert(m_isInitialized && "QR is not initialized.");
+  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
   const int rows = m_qr.rows();
   ei_assert(b.rows() == rows);
   result->resize(rows, b.cols());
@@ -276,27 +190,17 @@ bool QR<MatrixType>::solve(
   // Q^T without explicitly forming matrixQ(). Investigate.
   *result = matrixQ().transpose()*b;
   
-  if(!isSurjective())
-  {
-    // is result is in the image of R ?
-    RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff();
-    for(int col = 0; col < result->cols(); ++col)
-      for(int row = m_rank; row < result->rows(); ++row)
-        if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res))
-          return false;
-  }
-  m_qr.corner(TopLeft, m_rank, m_rank)
+  const int rank = std::min(result->rows(), result->cols());
+  m_qr.corner(TopLeft, rank, rank)
       .template marked<UpperTriangular>()
-      .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
-
-  return true;
+      .solveTriangularInPlace(result->corner(TopLeft, rank, result->cols()));
 }
 
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType QR<MatrixType>::matrixQ() const
+MatrixType HouseholderQR<MatrixType>::matrixQ() const
 {
-  ei_assert(m_isInitialized && "QR is not initialized.");
+  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
   // compute the product Q_0 Q_1 ... Q_n-1,
   // where Q_k is the k-th Householder transformation I - h_k v_k v_k'
   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
@@ -319,15 +223,15 @@ MatrixType QR<MatrixType>::matrixQ() const
 
 #endif // EIGEN_HIDE_HEAVY_CODE
 
-/** \return the QR decomposition of \c *this.
+/** \return the Householder QR decomposition of \c *this.
   *
-  * \sa class QR
+  * \sa class HouseholderQR
   */
 template<typename Derived>
-const QR<typename MatrixBase<Derived>::PlainMatrixType>
-MatrixBase<Derived>::qr() const
+const HouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
+MatrixBase<Derived>::householderQr() const
 {
-  return QR<PlainMatrixType>(eval());
+  return HouseholderQR<PlainMatrixType>(eval());
 }
 
 
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 9733fbd21..f9f9feb89 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -34,9 +34,7 @@
   *
   * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
   *
-  * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N
-  * with \c M \>= \c N.
-  *
+  * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N.
   *
   * \sa MatrixBase::SVD()
   */
@@ -55,13 +53,13 @@ template<typename MatrixType> class SVD
     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector;
 
-    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType;
+    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixUType;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType;
-    typedef Matrix<Scalar, MinSize, 1> SingularValuesType;
+    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> SingularValuesType;
 
   public:
 
-    /** 
+    /**
     * \brief Default Constructor.
     *
     * The default constructor is useful in cases in which the user intends to
@@ -70,9 +68,9 @@ template<typename MatrixType> class SVD
     SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {}
 
     SVD(const MatrixType& matrix)
-      : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
+      : m_matU(matrix.rows(), matrix.rows()),
         m_matV(matrix.cols(),matrix.cols()),
-        m_sigma(std::min(matrix.rows(),matrix.cols())),
+        m_sigma(matrix.cols()),
         m_isInitialized(false)
     {
       compute(matrix);
@@ -81,22 +79,22 @@ template<typename MatrixType> class SVD
     template<typename OtherDerived, typename ResultType>
     bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
 
-    const MatrixUType& matrixU() const 
-    { 
+    const MatrixUType& matrixU() const
+    {
       ei_assert(m_isInitialized && "SVD is not initialized.");
-      return m_matU; 
+      return m_matU;
     }
 
-    const SingularValuesType& singularValues() const 
+    const SingularValuesType& singularValues() const
     {
       ei_assert(m_isInitialized && "SVD is not initialized.");
-      return m_sigma; 
+      return m_sigma;
     }
 
-    const MatrixVType& matrixV() const 
+    const MatrixVType& matrixV() const
     {
       ei_assert(m_isInitialized && "SVD is not initialized.");
-      return m_matV; 
+      return m_matV;
     }
 
     void compute(const MatrixType& matrix);
@@ -112,6 +110,23 @@ template<typename MatrixType> class SVD
     void computeScalingRotation(ScalingType *positive, RotationType *unitary) const;
 
   protected:
+    // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow.
+    inline static Scalar pythag(Scalar a, Scalar b)
+    {
+      Scalar abs_a = ei_abs(a);
+      Scalar abs_b = ei_abs(b);
+      if (abs_a > abs_b)
+        return abs_a*ei_sqrt(Scalar(1.0)+ei_abs2(abs_b/abs_a));
+      else
+        return (abs_b == Scalar(0.0) ? Scalar(0.0) : abs_b*ei_sqrt(Scalar(1.0)+ei_abs2(abs_a/abs_b)));
+    }
+
+    inline static Scalar sign(Scalar a, Scalar b)
+    {
+      return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a));
+    }
+    
+  protected:
     /** \internal */
     MatrixUType m_matU;
     /** \internal */
@@ -123,380 +138,271 @@ template<typename MatrixType> class SVD
 
 /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
   *
-  * \note this code has been adapted from JAMA (public domain)
+  * \note this code has been adapted from Numerical Recipes, third edition.
   */
 template<typename MatrixType>
 void SVD<MatrixType>::compute(const MatrixType& matrix)
 {
   const int m = matrix.rows();
   const int n = matrix.cols();
-  const int nu = std::min(m,n);
 
-  m_matU.resize(m, nu);
+  m_matU.resize(m, m);
   m_matU.setZero();
-  m_sigma.resize(std::min(m,n));
+  m_sigma.resize(n);
   m_matV.resize(n,n);
 
-  RowVector e(n);
-  ColVector work(m);
-  MatrixType matA(matrix);
-  const bool wantu = true;
-  const bool wantv = true;
-  int i=0, j=0, k=0;
-
-  // Reduce A to bidiagonal form, storing the diagonal elements
-  // in s and the super-diagonal elements in e.
-  int nct = std::min(m-1,n);
-  int nrt = std::max(0,std::min(n-2,m));
-  for (k = 0; k < std::max(nct,nrt); ++k)
-  {
-    if (k < nct)
-    {
-      // Compute the transformation for the k-th column and
-      // place the k-th diagonal in m_sigma[k].
-      m_sigma[k] = matA.col(k).end(m-k).norm();
-      if (m_sigma[k] != 0.0) // FIXME
-      {
-        if (matA(k,k) < 0.0)
-          m_sigma[k] = -m_sigma[k];
-        matA.col(k).end(m-k) /= m_sigma[k];
-        matA(k,k) += 1.0;
-      }
-      m_sigma[k] = -m_sigma[k];
-    }
-
-    for (j = k+1; j < n; ++j)
-    {
-      if ((k < nct) && (m_sigma[k] != 0.0))
-      {
-        // Apply the transformation.
-        Scalar t = matA.col(k).end(m-k).dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ??
-        t = -t/matA(k,k);
-        matA.col(j).end(m-k) += t * matA.col(k).end(m-k);
-      }
+  int max_iters = 30;
 
-      // Place the k-th row of A into e for the
-      // subsequent calculation of the row transformation.
-      e[j] = matA(k,j);
-    }
+  MatrixVType& V = m_matV;
+  MatrixType A = matrix;
+  SingularValuesType& W = m_sigma;
 
-    // Place the transformation in U for subsequent back multiplication.
-    if (wantu & (k < nct))
-      m_matU.col(k).end(m-k) = matA.col(k).end(m-k);
+  bool flag;
+  int i,its,j,jj,k,l,nm;
+  Scalar anorm, c, f, g, h, s, scale, x, y, z;
+  bool convergence = true;
+  Scalar eps = precision<Scalar>();
 
-    if (k < nrt)
+  Matrix<Scalar,Dynamic,1> rv1(n);
+  g = scale = anorm = 0;
+  // Householder reduction to bidiagonal form.
+  for (i=0; i<n; i++)
+  {
+    l = i+2;
+    rv1[i] = scale*g;
+    g = s = scale = 0.0;
+    if (i < m)
     {
-      // Compute the k-th row transformation and place the
-      // k-th super-diagonal in e[k].
-      e[k] = e.end(n-k-1).norm();
-      if (e[k] != 0.0)
-      {
-          if (e[k+1] < 0.0)
-            e[k] = -e[k];
-          e.end(n-k-1) /= e[k];
-          e[k+1] += 1.0;
-      }
-      e[k] = -e[k];
-      if ((k+1 < m) & (e[k] != 0.0))
+      scale = A.col(i).end(m-i).cwise().abs().sum();
+      if (scale != Scalar(0))
       {
-        // Apply the transformation.
-        work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1);
-        for (j = k+1; j < n; ++j)
-          matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1);
+        for (k=i; k<m; k++)
+        {
+          A(k, i) /= scale;
+          s += A(k, i)*A(k, i);
+        }
+        f = A(i, i);
+        g = -sign( ei_sqrt(s), f );
+        h = f*g - s;
+        A(i, i)=f-g;
+        for (j=l-1; j<n; j++)
+        {
+          s = A.col(i).end(m-i).dot(A.col(j).end(m-i));
+          f = s/h;
+          A.col(j).end(m-i) += f*A.col(i).end(m-i);
+        }
+        A.col(i).end(m-i) *= scale;
       }
-
-      // Place the transformation in V for subsequent back multiplication.
-      if (wantv)
-        m_matV.col(k).end(n-k-1) = e.end(n-k-1);
-    }
-  }
-
-
-  // Set up the final bidiagonal matrix or order p.
-  int p = std::min(n,m+1);
-  if (nct < n)
-    m_sigma[nct] = matA(nct,nct);
-  if (m < p)
-    m_sigma[p-1] = 0.0;
-  if (nrt+1 < p)
-    e[nrt] = matA(nrt,p-1);
-  e[p-1] = 0.0;
-
-  // If required, generate U.
-  if (wantu)
-  {
-    for (j = nct; j < nu; ++j)
-    {
-      m_matU.col(j).setZero();
-      m_matU(j,j) = 1.0;
     }
-    for (k = nct-1; k >= 0; k--)
+    W[i] = scale * g;
+    g = s = scale = 0.0;
+    if (i+1 <= m && i+1 != n)
     {
-      if (m_sigma[k] != 0.0)
+      scale = A.row(i).end(n-l+1).cwise().abs().sum();
+      if (scale != Scalar(0))
       {
-        for (j = k+1; j < nu; ++j)
+        for (k=l-1; k<n; k++)
         {
-          Scalar t = m_matU.col(k).end(m-k).dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ?
-          t = -t/m_matU(k,k);
-          m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k);
+          A(i, k) /= scale;
+          s += A(i, k)*A(i, k);
         }
-        m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k);
-        m_matU(k,k) = Scalar(1) + m_matU(k,k);
-        if (k-1>0)
-          m_matU.col(k).start(k-1).setZero();
-      }
-      else
-      {
-        m_matU.col(k).setZero();
-        m_matU(k,k) = 1.0;
+        f = A(i,l-1);
+        g = -sign(ei_sqrt(s),f);
+        h = f*g - s;
+        A(i,l-1) = f-g;
+        rv1.end(n-l+1) = A.row(i).end(n-l+1)/h;
+        for (j=l-1; j<m; j++)
+        {
+          s = A.row(j).end(n-l+1).dot(A.row(i).end(n-l+1));
+          A.row(j).end(n-l+1) += s*rv1.end(n-l+1).transpose();
+        }
+        A.row(i).end(n-l+1) *= scale;
       }
     }
+    anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) );
   }
-
-  // If required, generate V.
-  if (wantv)
+  // Accumulation of right-hand transformations.
+  for (i=n-1; i>=0; i--)
   {
-    for (k = n-1; k >= 0; k--)
+    //Accumulation of right-hand transformations.
+    if (i < n-1)
     {
-      if ((k < nrt) & (e[k] != 0.0))
+      if (g != Scalar(0.0))
       {
-        for (j = k+1; j < nu; ++j)
+        for (j=l; j<n; j++) //Double division to avoid possible underflow.
+          V(j, i) = (A(i, j)/A(i, l))/g;
+        for (j=l; j<n; j++)
         {
-          Scalar t = m_matV.col(k).end(n-k-1).dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ?
-          t = -t/m_matV(k+1,k);
-          m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1);
+          s = A.row(i).end(n-l).dot(V.col(j).end(n-l));
+          V.col(j).end(n-l) += s * V.col(i).end(n-l);
         }
       }
-      m_matV.col(k).setZero();
-      m_matV(k,k) = 1.0;
+      V.row(i).end(n-l).setZero();
+      V.col(i).end(n-l).setZero();
     }
+    V(i, i) = 1.0;
+    g = rv1[i];
+    l = i;
   }
-
-  // Main iteration loop for the singular values.
-  int pp = p-1;
-  int iter = 0;
-  Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
-  while (p > 0)
+  // Accumulation of left-hand transformations.
+  for (i=std::min(m,n)-1; i>=0; i--)
   {
-    int k=0;
-    int kase=0;
-
-    // Here is where a test for too many iterations would go.
-
-    // This section of the program inspects for
-    // negligible elements in the s and e arrays.  On
-    // completion the variables kase and k are set as follows.
-
-    // kase = 1     if s(p) and e[k-1] are negligible and k<p
-    // kase = 2     if s(k) is negligible and k<p
-    // kase = 3     if e[k-1] is negligible, k<p, and
-    //              s(k), ..., s(p) are not negligible (qr step).
-    // kase = 4     if e(p-1) is negligible (convergence).
-
-    for (k = p-2; k >= -1; --k)
+    l = i+1;
+    g = W[i];
+    if (n-l>0)
+      A.row(i).end(n-l).setZero();
+    if (g != Scalar(0.0))
     {
-      if (k == -1)
-          break;
-      if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1])))
+      g = Scalar(1.0)/g;
+      for (j=l; j<n; j++)
       {
-          e[k] = 0.0;
-          break;
+        s = A.col(i).end(m-l).dot(A.col(j).end(m-l));
+        f = (s/A(i,i))*g;
+        A.col(j).end(m-i) += f * A.col(i).end(m-i);
       }
-    }
-    if (k == p-2)
-    {
-      kase = 4;
+      A.col(i).end(m-i) *= g;
     }
     else
+      A.col(i).end(m-i).setZero();
+    ++A(i,i);
+  }
+  // Diagonalization of the bidiagonal form: Loop over
+  // singular values, and over allowed iterations.
+  for (k=n-1; k>=0; k--)
+  {
+    for (its=0; its<max_iters; its++)
     {
-      int ks;
-      for (ks = p-1; ks >= k; --ks)
+      flag = true;
+      for (l=k; l>=0; l--)
       {
-        if (ks == k)
-          break;
-        Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0));
-        if (ei_abs(m_sigma[ks]) <= eps*t)
+        // Test for splitting.
+        nm = l-1;
+        // Note that rv1[1] is always zero.
+        //if ((double)(ei_abs(rv1[l])+anorm) == anorm)
+        if (l==0 || ei_abs(rv1[l]) <= eps*anorm)
         {
-          m_sigma[ks] = 0.0;
+          flag = false;
           break;
         }
+        //if ((double)(ei_abs(W[nm])+anorm) == anorm)
+        if (ei_abs(W[nm]) <= eps*anorm)
+          break;
       }
-      if (ks == k)
-      {
-        kase = 3;
-      }
-      else if (ks == p-1)
-      {
-        kase = 1;
-      }
-      else
-      {
-        kase = 2;
-        k = ks;
-      }
-    }
-    ++k;
-
-    // Perform the task indicated by kase.
-    switch (kase)
-    {
-
-      // Deflate negligible s(p).
-      case 1:
+      if (flag)
       {
-        Scalar f(e[p-2]);
-        e[p-2] = 0.0;
-        for (j = p-2; j >= k; --j)
+        c = 0.0;  //Cancellation of rv1[l], if l > 0.
+        s = 1.0;
+        for (i=l ;i<k+1; i++)
         {
-          Scalar t(ei_hypot(m_sigma[j],f));
-          Scalar cs(m_sigma[j]/t);
-          Scalar sn(f/t);
-          m_sigma[j] = t;
-          if (j != k)
-          {
-            f = -sn*e[j-1];
-            e[j-1] = cs*e[j-1];
-          }
-          if (wantv)
+          f = s*rv1[i];
+          rv1[i] = c*rv1[i];
+          //if ((double)(ei_abs(f)+anorm) == anorm)
+          if (ei_abs(f) <= eps*anorm)
+            break;
+          g = W[i];
+          h = pythag(f,g);
+          W[i] = h;
+          h = Scalar(1.0)/h;
+          c = g*h;
+          s = -f*h;
+          for (j=0; j<m; j++)
           {
-            for (i = 0; i < n; ++i)
-            {
-              t = cs*m_matV(i,j) + sn*m_matV(i,p-1);
-              m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1);
-              m_matV(i,j) = t;
-            }
+            y = A(j,nm);
+            z = A(j,i);
+            A(j,nm) = y*c + z*s;
+            A(j,i)  = z*c - y*s;
           }
         }
       }
-      break;
-
-      // Split at negligible s(k).
-      case 2:
+      z = W[k];
+      if (l == k)  //Convergence.
       {
-        Scalar f(e[k-1]);
-        e[k-1] = 0.0;
-        for (j = k; j < p; ++j)
-        {
-          Scalar t(ei_hypot(m_sigma[j],f));
-          Scalar cs( m_sigma[j]/t);
-          Scalar sn(f/t);
-          m_sigma[j] = t;
-          f = -sn*e[j];
-          e[j] = cs*e[j];
-          if (wantu)
-          {
-            for (i = 0; i < m; ++i)
-            {
-              t = cs*m_matU(i,j) + sn*m_matU(i,k-1);
-              m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1);
-              m_matU(i,j) = t;
-            }
-          }
+        if (z < 0.0) { // Singular value is made nonnegative.
+          W[k] = -z;
+          V.col(k) = -V.col(k);
         }
+        break;
       }
-      break;
-
-      // Perform one qr step.
-      case 3:
+      if (its+1 == max_iters)
       {
-        // Calculate the shift.
-        Scalar scale = std::max(std::max(std::max(std::max(
-                        ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])),
-                        ei_abs(m_sigma[k])),ei_abs(e[k]));
-        Scalar sp = m_sigma[p-1]/scale;
-        Scalar spm1 = m_sigma[p-2]/scale;
-        Scalar epm1 = e[p-2]/scale;
-        Scalar sk = m_sigma[k]/scale;
-        Scalar ek = e[k]/scale;
-        Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2);
-        Scalar c = (sp*epm1)*(sp*epm1);
-        Scalar shift = 0.0;
-        if ((b != 0.0) || (c != 0.0))
+        convergence = false;
+      }
+      x = W[l]; // Shift from bottom 2-by-2 minor.
+      nm = k-1;
+      y = W[nm];
+      g = rv1[nm];
+      h = rv1[k];
+      f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y);
+      g = pythag(f,1.0);
+      f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
+      c = s = 1.0;
+      //Next QR transformation:
+      for (j=l; j<=nm; j++)
+      {
+        i = j+1;
+        g = rv1[i];
+        y = W[i];
+        h = s*g;
+        g = c*g;
+        z = pythag(f,h);
+        rv1[j] = z;
+        c = f/z;
+        s = h/z;
+        f = x*c + g*s;
+        g = g*c - x*s;
+        h = y*s;
+        y *= c;
+        for (jj=0; jj<n; jj++)
         {
-          shift = ei_sqrt(b*b + c);
-          if (b < 0.0)
-            shift = -shift;
-          shift = c/(b + shift);
+          x = V(jj,j);
+          z = V(jj,i);
+          V(jj,j) = x*c + z*s;
+          V(jj,i) = z*c - x*s;
         }
-        Scalar f = (sk + sp)*(sk - sp) + shift;
-        Scalar g = sk*ek;
-
-        // Chase zeros.
-
-        for (j = k; j < p-1; ++j)
+        z = pythag(f,h);
+        W[j] = z;
+        // Rotation can be arbitrary if z = 0.
+        if (z!=Scalar(0))
         {
-          Scalar t = ei_hypot(f,g);
-          Scalar cs = f/t;
-          Scalar sn = g/t;
-          if (j != k)
-            e[j-1] = t;
-          f = cs*m_sigma[j] + sn*e[j];
-          e[j] = cs*e[j] - sn*m_sigma[j];
-          g = sn*m_sigma[j+1];
-          m_sigma[j+1] = cs*m_sigma[j+1];
-          if (wantv)
-          {
-            for (i = 0; i < n; ++i)
-            {
-              t = cs*m_matV(i,j) + sn*m_matV(i,j+1);
-              m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1);
-              m_matV(i,j) = t;
-            }
-          }
-          t = ei_hypot(f,g);
-          cs = f/t;
-          sn = g/t;
-          m_sigma[j] = t;
-          f = cs*e[j] + sn*m_sigma[j+1];
-          m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1];
-          g = sn*e[j+1];
-          e[j+1] = cs*e[j+1];
-          if (wantu && (j < m-1))
-          {
-            for (i = 0; i < m; ++i)
-            {
-              t = cs*m_matU(i,j) + sn*m_matU(i,j+1);
-              m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1);
-              m_matU(i,j) = t;
-            }
-          }
+          z = Scalar(1.0)/z;
+          c = f*z;
+          s = h*z;
         }
-        e[p-2] = f;
-        iter = iter + 1;
-      }
-      break;
-
-      // Convergence.
-      case 4:
-      {
-        // Make the singular values positive.
-        if (m_sigma[k] <= 0.0)
+        f = c*g + s*y;
+        x = c*y - s*g;
+        for (jj=0; jj<m; jj++)
         {
-          m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0);
-          if (wantv)
-            m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1);
+          y = A(jj,j);
+          z = A(jj,i);
+          A(jj,j) = y*c + z*s;
+          A(jj,i) = z*c - y*s;
         }
+      }
+      rv1[l] = 0.0;
+      rv1[k] = f;
+      W[k]   = x;
+    }
+  }
 
-        // Order the singular values.
-        while (k < pp)
-        {
-          if (m_sigma[k] >= m_sigma[k+1])
-            break;
-          Scalar t = m_sigma[k];
-          m_sigma[k] = m_sigma[k+1];
-          m_sigma[k+1] = t;
-          if (wantv && (k < n-1))
-            m_matV.col(k).swap(m_matV.col(k+1));
-          if (wantu && (k < m-1))
-            m_matU.col(k).swap(m_matU.col(k+1));
-          ++k;
-        }
-        iter = 0;
-        p--;
+  // sort the singular values:
+  {
+    for (int i=0; i<n; i++)
+    {
+      int k;
+      W.end(n-i).minCoeff(&k);
+      if (k != i)
+      {
+        std::swap(W[k],W[i]);
+        A.col(i).swap(A.col(k));
+        V.col(i).swap(V.col(k));
       }
-      break;
-    } // end big switch
-  } // end iterations
+    }
+  }
+  m_matU.setZero();
+  if (m>=n)
+    m_matU.block(0,0,m,n) = A;
+  else
+    m_matU = A.block(0,0,m,m);
 
   m_isInitialized = true;
 }
@@ -554,6 +460,8 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul
   const int rows = m_matU.rows();
   ei_assert(b.rows() == rows);
 
+  result->resize(m_matV.rows(), b.cols());
+
   Scalar maxVal = m_sigma.cwise().abs().maxCoeff();
   for (int j=0; j<b.cols(); ++j)
   {
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index b751b5cb9..af3b5e5eb 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -443,9 +443,8 @@ class SparseMatrix
         // two passes algorithm:
         //  1 - compute the number of coeffs per dest inner vector
         //  2 - do the actual copy/eval
-        // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed
-        //typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
-        typedef typename ei_eval<OtherDerived>::type OtherCopy;
+        // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed
+        typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
         typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
         OtherCopy otherCopy(other.derived());
 
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index d62d23a78..1af452251 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -102,8 +102,10 @@ template<typename Derived> class SparseMatrixBase
     /** \internal the return type of MatrixBase::imag() */
     typedef SparseCwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
     /** \internal the return type of MatrixBase::adjoint() */
-    typedef SparseTranspose</*NestByValue<*/typename ei_cleantype<ConjugateReturnType>::type> /*>*/
-            AdjointReturnType;
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+                        SparseCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, SparseNestByValue<Eigen::SparseTranspose<Derived> > >,
+                        SparseTranspose<Derived>
+                     >::ret AdjointReturnType;
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** This is the "real scalar" type; if the \a Scalar type is already real numbers
@@ -357,7 +359,7 @@ template<typename Derived> class SparseMatrixBase
     SparseTranspose<Derived> transpose() { return derived(); }
     const SparseTranspose<Derived> transpose() const { return derived(); }
     // void transposeInPlace();
-    const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
+    const AdjointReturnType adjoint() const { return transpose().nestByValue(); }
 
     // sub-vector
     SparseInnerVectorSet<Derived,1> row(int i);
@@ -529,7 +531,7 @@ template<typename Derived> class SparseMatrixBase
       */
 //     inline int stride(void) const { return derived().stride(); }
 
-//     inline const NestByValue<Derived> nestByValue() const;
+    inline const SparseNestByValue<Derived> nestByValue() const;
 
 
     ConjugateReturnType conjugate() const;
diff --git a/Eigen/src/Sparse/SparseNestByValue.h b/Eigen/src/Sparse/SparseNestByValue.h
new file mode 100644
index 000000000..b48277232
--- /dev/null
+++ b/Eigen/src/Sparse/SparseNestByValue.h
@@ -0,0 +1,84 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2006-2008 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_SPARSENESTBYVALUE_H
+#define EIGEN_SPARSENESTBYVALUE_H
+
+/** \class SparseNestByValue
+  *
+  * \brief Expression which must be nested by value
+  *
+  * \param ExpressionType the type of the object of which we are requiring nesting-by-value
+  *
+  * This class is the return type of MatrixBase::nestByValue()
+  * and most of the time this is the only way it is used.
+  *
+  * \sa SparseMatrixBase::nestByValue(), class NestByValue
+  */
+template<typename ExpressionType>
+struct ei_traits<SparseNestByValue<ExpressionType> > : public ei_traits<ExpressionType>
+{};
+
+template<typename ExpressionType> class SparseNestByValue
+  : public SparseMatrixBase<SparseNestByValue<ExpressionType> >
+{
+  public:
+
+    typedef typename ExpressionType::InnerIterator InnerIterator;
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue)
+
+    inline SparseNestByValue(const ExpressionType& matrix) : m_expression(matrix) {}
+
+    EIGEN_STRONG_INLINE int rows() const { return m_expression.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return m_expression.cols(); }
+
+    operator const ExpressionType&() const { return m_expression; }
+
+  protected:
+    const ExpressionType m_expression;
+};
+
+/** \returns an expression of the temporary version of *this.
+  */
+template<typename Derived>
+inline const SparseNestByValue<Derived>
+SparseMatrixBase<Derived>::nestByValue() const
+{
+  return SparseNestByValue<Derived>(derived());
+}
+
+// template<typename MatrixType>
+// class SparseNestByValue<MatrixType>::InnerIterator : public MatrixType::InnerIterator
+// {
+//     typedef typename MatrixType::InnerIterator Base;
+//   public:
+// 
+//     EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer)
+//       : Base(expr.m_expression, outer)
+//     {}
+// };
+
+#endif // EIGEN_SPARSENESTBYVALUE_H
diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h
index 52781aa46..b5fc7c7b7 100644
--- a/Eigen/src/Sparse/SparseUtil.h
+++ b/Eigen/src/Sparse/SparseUtil.h
@@ -106,6 +106,7 @@ template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix;
 template<typename _Scalar, int _Flags = 0> class SparseVector;
 template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
 
+template<typename MatrixType>                            class SparseNestByValue;
 template<typename MatrixType>                            class SparseTranspose;
 template<typename MatrixType, int Size>                  class SparseInnerVectorSet;
 template<typename Derived>                               class SparseCwise;
@@ -146,4 +147,6 @@ template<typename T> class ei_eval<T,IsSparse>
     typedef SparseMatrix<_Scalar, _Flags> type;
 };
 
+template<typename T> struct ei_must_nest_by_value<SparseNestByValue<T> > { enum { ret = true }; };
+
 #endif // EIGEN_SPARSEUTIL_H
diff --git a/cmake/FindEigen2.cmake b/cmake/FindEigen2.cmake
index 0e86aaf29..ee054b0da 100644
--- a/cmake/FindEigen2.cmake
+++ b/cmake/FindEigen2.cmake
@@ -46,7 +46,7 @@ macro(_eigen2_check_version)
   endif(${EIGEN2_VERSION} VERSION_LESS ${Eigen2_FIND_VERSION})
 
   if(NOT EIGEN2_VERSION_OK)
-  
+
     message(STATUS "Eigen2 version ${EIGEN2_VERSION} found in ${EIGEN2_INCLUDE_DIR}, "
                    "but at least version ${Eigen2_FIND_VERSION} is required")
   endif(NOT EIGEN2_VERSION_OK)
diff --git a/doc/D09_StructHavingEigenMembers.dox b/doc/D09_StructHavingEigenMembers.dox
index 4ff2b23aa..07b6570dc 100644
--- a/doc/D09_StructHavingEigenMembers.dox
+++ b/doc/D09_StructHavingEigenMembers.dox
@@ -88,7 +88,7 @@ The solution is to let class Foo have an aligned "operator new", as we showed in
 
 \section movetotop Should I then put all the members of Eigen types at the beginning of my class?
 
-No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like
+That's not required. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So code like this works fine:
 
 \code
 class Foo
@@ -100,25 +100,13 @@ public:
 };
 \endcode
 
-it will work just fine. You do \b not need to rewrite it as
-
-\code
-class Foo
-{
-  Eigen::Vector2d v;
-  double x;
-public:
-  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-};
-\endcode
-
 \section dynamicsize What about dynamic-size matrices and vectors?
 
-Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors.
+Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with \ref FixedSizeVectorizable  "fixed-size vectorizable matrices and vectors".
 
 \section bugineigen So is this a bug in Eigen?
 
-No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data.
+No, it's not our bug. It's more like an inherent problem of the C++98 language specification, and seems to be taken care of in the upcoming language revision: <a href="http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2341.pdf">see this document</a>.
 
 \section conditional What if I want to do this conditionnally (depending on template parameters) ?
 
diff --git a/doc/D11_UnalignedArrayAssert.dox b/doc/D11_UnalignedArrayAssert.dox
index da1b49a08..e95c767b4 100644
--- a/doc/D11_UnalignedArrayAssert.dox
+++ b/doc/D11_UnalignedArrayAssert.dox
@@ -15,12 +15,24 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
 There are 4 known causes for this issue. Please read on to understand them and learn how to fix them.
 
 \b Table \b of \b contents
+ - \ref where
  - \ref c1
  - \ref c2
  - \ref c3
  - \ref c4
  - \ref explanation
 
+\section where Where in my own code is the cause of the problem?
+
+First of all, you need to find out where in your own code this assertion was triggered from. At first glance, the error message doesn't look helpful, as it refers to a file inside Eigen! However, since your program crashed, if you can reproduce the crash, you can get a backtrace using any debugger. For example, if you're using GCC, you can use the GDB debugger as follows:
+\code
+$ gdb ./my_program          # Start GDB on your program
+> run                       # Start running your program
+...                         # Now reproduce the crash!
+> bt                        # Obtain the backtrace
+\endcode
+Now that you know precisely where in your own code the problem is happening, read on to understand what you need to change.
+
 \section c1 Cause 1: Structures having Eigen objects as members
 
 If you have code like this,
diff --git a/doc/eigendoxy_header.html.in b/doc/eigendoxy_header.html.in
index 97cafd41d..572c47158 100644
--- a/doc/eigendoxy_header.html.in
+++ b/doc/eigendoxy_header.html.in
@@ -6,5 +6,5 @@
 </head><body>
 <a name="top"></a>
 <a class="logo" href="http://eigen.tuxfamily.org/">
-<img src="Eigen_Silly_Professor_64x64.png" width=64 height=64 alt="Eiegn's silly professor"
+<img src="Eigen_Silly_Professor_64x64.png" width=64 height=64 alt="Eigen's silly professor"
      style="position:absolute; border:none" /></a>
\ No newline at end of file
diff --git a/doc/examples/class_FixedVectorBlock.cpp b/doc/examples/class_FixedVectorBlock.cpp
new file mode 100644
index 000000000..c176be495
--- /dev/null
+++ b/doc/examples/class_FixedVectorBlock.cpp
@@ -0,0 +1,26 @@
+#include <Eigen/Core>
+USING_PART_OF_NAMESPACE_EIGEN
+using namespace std;
+
+template<typename Derived>
+Eigen::VectorBlock<Derived, 2>
+firstTwo(MatrixBase<Derived>& v)
+{
+  return Eigen::VectorBlock<Derived, 2>(v.derived(), 0);
+}
+
+template<typename Derived>
+const Eigen::VectorBlock<Derived, 2>
+firstTwo(const MatrixBase<Derived>& v)
+{
+  return Eigen::VectorBlock<Derived, 2>(v.derived(), 0);
+}
+
+int main(int, char**)
+{
+  Matrix<int,1,6> v; v << 1,2,3,4,5,6;
+  cout << firstTwo(4*v) << endl; // calls the const version
+  firstTwo(v) *= 2;              // calls the non-const version
+  cout << "Now the vector v is:" << endl << v << endl;
+  return 0;
+}
diff --git a/doc/examples/class_VectorBlock.cpp b/doc/examples/class_VectorBlock.cpp
new file mode 100644
index 000000000..d979f973a
--- /dev/null
+++ b/doc/examples/class_VectorBlock.cpp
@@ -0,0 +1,26 @@
+#include <Eigen/Core>
+USING_PART_OF_NAMESPACE_EIGEN
+using namespace std;
+
+template<typename Derived>
+Eigen::VectorBlock<Derived>
+segmentFromRange(MatrixBase<Derived>& v, int start, int end)
+{
+  return Eigen::VectorBlock<Derived>(v.derived(), start, end-start);
+}
+
+template<typename Derived>
+const Eigen::VectorBlock<Derived>
+segmentFromRange(const MatrixBase<Derived>& v, int start, int end)
+{
+  return Eigen::VectorBlock<Derived>(v.derived(), start, end-start);
+}
+
+int main(int, char**)
+{
+  Matrix<int,1,6> v; v << 1,2,3,4,5,6;
+  cout << segmentFromRange(2*v, 2, 4) << endl; // calls the const version
+  segmentFromRange(v, 1, 3) *= 5;              // calls the non-const version
+  cout << "Now the vector v is:" << endl << v << endl;
+  return 0;
+}
diff --git a/doc/snippets/QR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp
index 7e629f851..aa9404951 100644
--- a/doc/snippets/QR_solve.cpp
+++ b/doc/snippets/HouseholderQR_solve.cpp
@@ -4,11 +4,6 @@ Matrix3f y = Matrix3f::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is the matrix y:" << endl << y << endl;
 Matrix3f x;
-if(m.qr().solve(y, &x))
-{
-  assert(y.isApprox(m*x));
-  cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
-}
-else
-  cout << "The equation mx=y does not have any solution." << endl;
-
+m.householderQr().solve(y, &x));
+assert(y.isApprox(m*x));
+cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
diff --git a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
new file mode 100644
index 000000000..19e24c90b
--- /dev/null
+++ b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
@@ -0,0 +1,9 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d inv;
+if(m.computeInverseWithCheck(&inv)) {
+  cout << "It is invertible, and its inverse is:" << endl << inv << endl;
+}
+else {
+  cout << "It is not invertible." << endl;
+}
diff --git a/test/adjoint.cpp b/test/adjoint.cpp
index 965e4d256..403f16a45 100644
--- a/test/adjoint.cpp
+++ b/test/adjoint.cpp
@@ -76,6 +76,7 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   {
     VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(),         static_cast<RealScalar>(1));
     VERIFY_IS_APPROX(v1.norm(),                       v1.stableNorm());
+    VERIFY_IS_APPROX(v1.blueNorm(),                       v1.stableNorm());
   }
 
   // check compatibility of dot and adjoint
diff --git a/test/array.cpp b/test/array.cpp
index 37deeaa4f..a308f7366 100644
--- a/test/array.cpp
+++ b/test/array.cpp
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp
index 2f088988f..52aff93ee 100644
--- a/test/basicstuff.cpp
+++ b/test/basicstuff.cpp
@@ -107,7 +107,7 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   {
     VERIFY_IS_NOT_APPROX(m3, m1);
   }
-  
+
   m3.real() = m1.real();
   VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), static_cast<const MatrixType&>(m1).real());
   VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), m1.real());
@@ -121,16 +121,16 @@ template<typename MatrixType> void basicStuffComplex(const MatrixType& m)
 
   int rows = m.rows();
   int cols = m.cols();
-  
+
   Scalar s1 = ei_random<Scalar>(),
          s2 = ei_random<Scalar>();
-  
+
   VERIFY(ei_real(s1)==ei_real_ref(s1));
   VERIFY(ei_imag(s1)==ei_imag_ref(s1));
   ei_real_ref(s1) = ei_real(s2);
   ei_imag_ref(s1) = ei_imag(s2);
   VERIFY(s1==s2);
-  
+
   RealMatrixType rm1 = RealMatrixType::Random(rows,cols),
                  rm2 = RealMatrixType::Random(rows,cols);
   MatrixType cm(rows,cols);
@@ -162,7 +162,7 @@ void test_basicstuff()
     CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) );
     CALL_SUBTEST( basicStuff(Matrix<float, 100, 100>()) );
     CALL_SUBTEST( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
-    
+
     CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) );
     CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) );
   }
diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp
index c5d2715db..f8281a16b 100644
--- a/test/geo_hyperplane.cpp
+++ b/test/geo_hyperplane.cpp
@@ -64,7 +64,7 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane)
   // transform
   if (!NumTraits<Scalar>::IsComplex)
   {
-    MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ();
+    MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ();
     DiagonalMatrix<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random());
     Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random());
 
diff --git a/test/inverse.cpp b/test/inverse.cpp
index 90b897197..5ac39e35a 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp
@@ -65,6 +65,21 @@ template<typename MatrixType> void inverse(const MatrixType& m)
 
   // since for the general case we implement separately row-major and col-major, test that
   VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose());
+
+  //computeInverseWithCheck tests
+  //First: an invertible matrix
+  bool invertible = m1.computeInverseWithCheck(&m2);
+  VERIFY(invertible);
+  VERIFY_IS_APPROX(identity, m1*m2);
+
+  //Second: a rank one matrix (not invertible, except for 1x1 matrices)
+  VectorType v3 = VectorType::Random(rows);
+  MatrixType m3 = v3*v3.transpose(), m4(rows,cols);
+  invertible = m3.computeInverseWithCheck( &m4 );
+  if( 1 == rows ){
+      VERIFY( invertible ); }
+  else{
+      VERIFY( !invertible ); }
 }
 
 void test_inverse()
@@ -77,7 +92,7 @@ void test_inverse()
     CALL_SUBTEST( inverse(MatrixXf(8,8)) );
     CALL_SUBTEST( inverse(MatrixXcd(7,7)) );
   }
-  
+
   // test some tricky cases for 4x4 matrices
   VERIFY_IS_APPROX((Matrix4f() << 0,0,1,0, 1,0,0,0, 0,1,0,0, 0,0,0,1).finished().inverse(),
                    (Matrix4f() << 0,1,0,0, 0,0,1,0, 1,0,0,0, 0,0,0,1).finished());
diff --git a/test/lu.cpp b/test/lu.cpp
index 60e45ff2b..4ad92bb11 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -57,6 +57,11 @@ template<typename MatrixType> void lu_non_invertible()
   VERIFY_IS_APPROX(m3, m1*m2);
   m3 = MatrixType::Random(rows,cols2);
   VERIFY(!lu.solve(m3, &m2));
+  
+  typedef Matrix<typename MatrixType::Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
+  SquareMatrixType m4(rows, rows), m5(rows, rows);
+  createRandomMatrixOfRank(rows/2, rows, rows, m4);
+  VERIFY(!m4.computeInverseWithCheck(&m5));
 }
 
 template<typename MatrixType> void lu_invertible()
diff --git a/test/main.h b/test/main.h
index 76572d9b3..3947451cc 100644
--- a/test/main.h
+++ b/test/main.h
@@ -242,8 +242,8 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri
   const int diag_size = std::min(d.rows(),d.cols());
   d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
 
-  QR<MatrixType> qra(a);
-  QR<MatrixType> qrb(b);
+  HouseholderQR<MatrixType> qra(a);
+  HouseholderQR<MatrixType> qrb(b);
   m = (qra.matrixQ() * d * qrb.matrixQ()).lazy();
 }
 
diff --git a/test/qr.cpp b/test/qr.cpp
index 6e96f1e97..88a447c4b 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -36,7 +36,7 @@ template<typename MatrixType> void qr(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
   MatrixType a = MatrixType::Random(rows,cols);
-  QR<MatrixType> qrOfA(a);
+  HouseholderQR<MatrixType> qrOfA(a);
   VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR());
   VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR());
 
@@ -55,40 +55,6 @@ template<typename MatrixType> void qr(const MatrixType& m)
   VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
 }
 
-template<typename MatrixType> void qr_non_invertible()
-{
-  /* this test covers the following files: QR.h */
-  int rows = ei_random<int>(20,200), cols = ei_random<int>(20,rows), cols2 = ei_random<int>(20,rows);
-  int rank = ei_random<int>(1, std::min(rows, cols)-1);
-
-  MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
-  createRandomMatrixOfRank(rank, rows, cols, m1);
-
-  QR<MatrixType> lu(m1);
-//   typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel();
-//   typename LU<MatrixType>::ImageResultType m1image = lu.image();
-  std::cerr << rows << "x" << cols << "   " << rank << " " << lu.rank() << "\n";
-  if (rank != lu.rank())
-    std::cerr << lu.matrixR().diagonal().transpose() << "\n";
-  VERIFY(rank == lu.rank());
-  VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
-  VERIFY(!lu.isInjective());
-  VERIFY(!lu.isInvertible());
-  VERIFY(lu.isSurjective() == (lu.rank() == rows));
-//   VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
-//   VERIFY(m1image.lu().rank() == rank);
-//   MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
-//   sidebyside << m1, m1image;
-//   VERIFY(sidebyside.lu().rank() == rank);
-  m2 = MatrixType::Random(cols,cols2);
-  m3 = m1*m2;
-  m2 = MatrixType::Random(cols,cols2);
-  lu.solve(m3, &m2);
-  VERIFY_IS_APPROX(m3, m1*m2);
-  m3 = MatrixType::Random(rows,cols2);
-  VERIFY(!lu.solve(m3, &m2));
-}
-
 template<typename MatrixType> void qr_invertible()
 {
   /* this test covers the following files: QR.h */
@@ -105,33 +71,18 @@ template<typename MatrixType> void qr_invertible()
     m1 += a * a.adjoint();
   }
 
-  QR<MatrixType> lu(m1);
-  VERIFY(0 == lu.dimensionOfKernel());
-  VERIFY(size == lu.rank());
-  VERIFY(lu.isInjective());
-  VERIFY(lu.isSurjective());
-  VERIFY(lu.isInvertible());
-//   VERIFY(lu.image().lu().isInvertible());
+  HouseholderQR<MatrixType> qr(m1);
   m3 = MatrixType::Random(size,size);
-  lu.solve(m3, &m2);
+  qr.solve(m3, &m2);
   //std::cerr << m3 - m1*m2 << "\n\n";
   VERIFY_IS_APPROX(m3, m1*m2);
-//   VERIFY_IS_APPROX(m2, lu.inverse()*m3);
-  m3 = MatrixType::Random(size,size);
-  VERIFY(lu.solve(m3, &m2));
 }
 
 template<typename MatrixType> void qr_verify_assert()
 {
   MatrixType tmp;
 
-  QR<MatrixType> qr;
-  VERIFY_RAISES_ASSERT(qr.isFullRank())
-  VERIFY_RAISES_ASSERT(qr.rank())
-  VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
-  VERIFY_RAISES_ASSERT(qr.isInjective())
-  VERIFY_RAISES_ASSERT(qr.isSurjective())
-  VERIFY_RAISES_ASSERT(qr.isInvertible())
+  HouseholderQR<MatrixType> qr;
   VERIFY_RAISES_ASSERT(qr.matrixR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
   VERIFY_RAISES_ASSERT(qr.matrixQ())
@@ -149,11 +100,6 @@ void test_qr()
   }
 
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST( qr_non_invertible<MatrixXf>() );
-    CALL_SUBTEST( qr_non_invertible<MatrixXd>() );
-    // TODO fix issue with complex
-//     CALL_SUBTEST( qr_non_invertible<MatrixXcf>() );
-//     CALL_SUBTEST( qr_non_invertible<MatrixXcd>() );
     CALL_SUBTEST( qr_invertible<MatrixXf>() );
     CALL_SUBTEST( qr_invertible<MatrixXd>() );
     // TODO fix issue with complex
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 9b7d7c4e6..666eea872 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -299,6 +299,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     initSparse<Scalar>(density, refMat2, m2);
     VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
     VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
+
+    VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
   }
   
   // test prune
diff --git a/test/svd.cpp b/test/svd.cpp
index 2eb7881a0..93e5ea017 100644
--- a/test/svd.cpp
+++ b/test/svd.cpp
@@ -50,7 +50,7 @@ template<typename MatrixType> void svd(const MatrixType& m)
     MatrixType sigma = MatrixType::Zero(rows,cols);
     MatrixType matU  = MatrixType::Zero(rows,rows);
     sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal();
-    matU.block(0,0,rows,cols) = svd.matrixU();
+    matU = svd.matrixU();
     VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose());
   }