Replace local variables by member variables in compute() methods.

This is to avoid dynamic memory allocations in the compute() methods of
ComplexEigenSolver, EigenSolver, and SelfAdjointEigenSolver where possible.
As a result, Tridiagonalization::decomposeInPlace() is no longer used.
Biggest remaining issue is the allocation in HouseholderSequence::evalTo().

4 files changed, 101 insertions, 81 deletions

diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index 8e5f1310a..f6b90d70e 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -3,6 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -99,7 +100,8 @@ template<typename _MatrixType> class ComplexEigenSolver
             : m_eivec(),
               m_eivalues(),
               m_schur(),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matX()
     {}
     
     /** \brief Default Constructor with memory preallocation
@@ -112,7 +114,8 @@ template<typename _MatrixType> class ComplexEigenSolver
             : m_eivec(size, size),
               m_eivalues(size),
               m_schur(size),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matX(size, size)
     {}
 
     /** \brief Constructor; computes eigendecomposition of given matrix. 
@@ -125,7 +128,8 @@ template<typename _MatrixType> class ComplexEigenSolver
             : m_eivec(matrix.rows(),matrix.cols()),
               m_eivalues(matrix.cols()),
               m_schur(matrix.rows()),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matX(matrix.rows(),matrix.cols())
     {
       compute(matrix);
     }
@@ -199,6 +203,7 @@ template<typename _MatrixType> class ComplexEigenSolver
     EigenvalueType m_eivalues;
     ComplexSchur<MatrixType> m_schur;
     bool m_isInitialized;
+    EigenvectorType m_matX;
 };
 
 
@@ -217,16 +222,16 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
 
   // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
   // The matrix X is unit triangular.
-  EigenvectorType X = EigenvectorType::Zero(n, n);
+  m_matX = EigenvectorType::Zero(n, n);
   for(int k=n-1 ; k>=0 ; k--)
   {
-    X.coeffRef(k,k) = ComplexScalar(1.0,0.0);
+    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
     for(int i=k-1 ; i>=0 ; i--)
     {
-      X.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
+      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
       if(k-i-1>0)
-        X.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
+        m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
       if(z==ComplexScalar(0))
       {
@@ -234,12 +239,12 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
 	// Use a small value instead, to prevent division by zero.
         ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
       }
-      X.coeffRef(i,k) = X.coeff(i,k) / z;
+      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
     }
   }
 
   // Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
-  m_eivec = m_schur.matrixU() * X;
+  m_eivec.noalias() = m_schur.matrixU() * m_matX;
   // .. and normalize the eigenvectors
   for(int k=0 ; k<n ; k++)
   {
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index f8b953d9b..7713e04b9 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -120,7 +120,7 @@ template<typename _MatrixType> class EigenSolver
       *
       * \sa compute() for an example.
       */
-    EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
+ EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
 
     /** \brief Default Constructor with memory preallocation
       *
@@ -131,7 +131,11 @@ template<typename _MatrixType> class EigenSolver
     EigenSolver(int size)
       : m_eivec(size, size),
         m_eivalues(size),
-        m_isInitialized(false) {}
+        m_isInitialized(false),
+        m_realSchur(size),
+        m_matT(size, size), 
+        m_tmp(size)
+    {}
 
     /** \brief Constructor; computes eigendecomposition of given matrix. 
       * 
@@ -148,7 +152,10 @@ template<typename _MatrixType> class EigenSolver
     EigenSolver(const MatrixType& matrix)
       : m_eivec(matrix.rows(), matrix.cols()),
         m_eivalues(matrix.cols()),
-        m_isInitialized(false)
+        m_isInitialized(false),
+        m_realSchur(matrix.cols()),
+        m_matT(matrix.rows(), matrix.cols()), 
+        m_tmp(matrix.cols())
     {
       compute(matrix);
     }
@@ -261,12 +268,17 @@ template<typename _MatrixType> class EigenSolver
     EigenSolver& compute(const MatrixType& matrix);
 
   private:
-    void computeEigenvectors(MatrixType& matH);
+    void computeEigenvectors();
 
   protected:
     MatrixType m_eivec;
     EigenvalueType m_eivalues;
     bool m_isInitialized;
+    RealSchur<MatrixType> m_realSchur;
+    MatrixType m_matT;
+
+    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
+    ColumnVectorType m_tmp;
 };
 
 template<typename MatrixType>
@@ -324,32 +336,32 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
   assert(matrix.cols() == matrix.rows());
 
   // Reduce to real Schur form.
-  RealSchur<MatrixType> rs(matrix);
-  MatrixType matT = rs.matrixT();
-  m_eivec = rs.matrixU();
+  m_realSchur.compute(matrix);
+  m_matT = m_realSchur.matrixT();
+  m_eivec = m_realSchur.matrixU();
 
   // Compute eigenvalues from matT
   m_eivalues.resize(matrix.cols());
   int i = 0;
   while (i < matrix.cols()) 
   {
-    if (i == matrix.cols() - 1 || matT.coeff(i+1, i) == Scalar(0)) 
+    if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
     {
-      m_eivalues.coeffRef(i) = matT.coeff(i, i);
+      m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
       ++i;
     }
     else
     {
-      Scalar p = Scalar(0.5) * (matT.coeff(i, i) - matT.coeff(i+1, i+1));
-      Scalar z = ei_sqrt(ei_abs(p * p + matT.coeff(i+1, i) * matT.coeff(i, i+1)));
-      m_eivalues.coeffRef(i)   = ComplexScalar(matT.coeff(i+1, i+1) + p, z);
-      m_eivalues.coeffRef(i+1) = ComplexScalar(matT.coeff(i+1, i+1) + p, -z);
+      Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
+      Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
+      m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
+      m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
       i += 2;
     }
   }
   
   // Compute eigenvectors.
-  computeEigenvectors(matT);
+  computeEigenvectors();
 
   m_isInitialized = true;
   return *this;
@@ -376,7 +388,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
 
 
 template<typename MatrixType>
-void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
+void EigenSolver<MatrixType>::computeEigenvectors()
 {
   const int size = m_eivec.cols();
   const Scalar eps = NumTraits<Scalar>::epsilon();
@@ -385,7 +397,7 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
   Scalar norm = 0.0;
   for (int j = 0; j < size; ++j)
   {
-    norm += matH.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum();
+    norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum();
   }
   
   // Backsubstitute to find vectors of upper triangular form
@@ -405,11 +417,11 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
       Scalar lastr=0, lastw=0;
       int l = n;
 
-      matH.coeffRef(n,n) = 1.0;
+      m_matT.coeffRef(n,n) = 1.0;
       for (int i = n-1; i >= 0; i--)
       {
-        Scalar w = matH.coeff(i,i) - p;
-        Scalar r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1));
+        Scalar w = m_matT.coeff(i,i) - p;
+        Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
 
         if (m_eivalues.coeff(i).imag() < 0.0)
         {
@@ -422,27 +434,27 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
           if (m_eivalues.coeff(i).imag() == 0.0)
           {
             if (w != 0.0)
-              matH.coeffRef(i,n) = -r / w;
+              m_matT.coeffRef(i,n) = -r / w;
             else
-              matH.coeffRef(i,n) = -r / (eps * norm);
+              m_matT.coeffRef(i,n) = -r / (eps * norm);
           }
           else // Solve real equations
           {
-            Scalar x = matH.coeff(i,i+1);
-            Scalar y = matH.coeff(i+1,i);
+            Scalar x = m_matT.coeff(i,i+1);
+            Scalar y = m_matT.coeff(i+1,i);
             Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
             Scalar t = (x * lastr - lastw * r) / denom;
-            matH.coeffRef(i,n) = t;
+            m_matT.coeffRef(i,n) = t;
             if (ei_abs(x) > ei_abs(lastw))
-              matH.coeffRef(i+1,n) = (-r - w * t) / x;
+              m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
             else
-              matH.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
+              m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
           }
 
           // Overflow control
-          Scalar t = ei_abs(matH.coeff(i,n));
+          Scalar t = ei_abs(m_matT.coeff(i,n));
           if ((eps * t) * t > 1)
-            matH.col(n).tail(size-i) /= t;
+            m_matT.col(n).tail(size-i) /= t;
         }
       }
     }
@@ -452,24 +464,24 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
       int l = n-1;
 
       // Last vector component imaginary so matrix is triangular
-      if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n)))
+      if (ei_abs(m_matT.coeff(n,n-1)) > ei_abs(m_matT.coeff(n-1,n)))
       {
-        matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1);
-        matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1);
+        m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
+        m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
       }
       else
       {
-	std::complex<Scalar> cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q);
-        matH.coeffRef(n-1,n-1) = ei_real(cc);
-        matH.coeffRef(n-1,n) = ei_imag(cc);
+	std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
+        m_matT.coeffRef(n-1,n-1) = ei_real(cc);
+        m_matT.coeffRef(n-1,n) = ei_imag(cc);
       }
-      matH.coeffRef(n,n-1) = 0.0;
-      matH.coeffRef(n,n) = 1.0;
+      m_matT.coeffRef(n,n-1) = 0.0;
+      m_matT.coeffRef(n,n) = 1.0;
       for (int i = n-2; i >= 0; i--)
       {
-        Scalar ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1));
-        Scalar sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1));
-        Scalar w = matH.coeff(i,i) - p;
+        Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
+        Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
+        Scalar w = m_matT.coeff(i,i) - p;
 
         if (m_eivalues.coeff(i).imag() < 0.0)
         {
@@ -483,39 +495,39 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
           if (m_eivalues.coeff(i).imag() == 0)
           {
             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
-            matH.coeffRef(i,n-1) = ei_real(cc);
-            matH.coeffRef(i,n) = ei_imag(cc);
+            m_matT.coeffRef(i,n-1) = ei_real(cc);
+            m_matT.coeffRef(i,n) = ei_imag(cc);
           }
           else
           {
             // Solve complex equations
-            Scalar x = matH.coeff(i,i+1);
-            Scalar y = matH.coeff(i+1,i);
+            Scalar x = m_matT.coeff(i,i+1);
+            Scalar y = m_matT.coeff(i+1,i);
             Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
             Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
             if ((vr == 0.0) && (vi == 0.0))
               vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(lastw));
 
 	    std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
-            matH.coeffRef(i,n-1) = ei_real(cc);
-            matH.coeffRef(i,n) = ei_imag(cc);
+            m_matT.coeffRef(i,n-1) = ei_real(cc);
+            m_matT.coeffRef(i,n) = ei_imag(cc);
             if (ei_abs(x) > (ei_abs(lastw) + ei_abs(q)))
             {
-              matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x;
-              matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x;
+              m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
+              m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
             }
             else
             {
-              cc = cdiv(-lastra-y*matH.coeff(i,n-1),-lastsa-y*matH.coeff(i,n),lastw,q);
-              matH.coeffRef(i+1,n-1) = ei_real(cc);
-              matH.coeffRef(i+1,n) = ei_imag(cc);
+              cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
+              m_matT.coeffRef(i+1,n-1) = ei_real(cc);
+              m_matT.coeffRef(i+1,n) = ei_imag(cc);
             }
           }
 
           // Overflow control
-          Scalar t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n)));
+          Scalar t = std::max(ei_abs(m_matT.coeff(i,n-1)),ei_abs(m_matT.coeff(i,n)));
           if ((eps * t) * t > 1)
-            matH.block(i, n-1, size-i, 2) /= t;
+            m_matT.block(i, n-1, size-i, 2) /= t;
 
         }
       }
@@ -525,7 +537,8 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
   // Back transformation to get eigenvectors of original matrix
   for (int j = size-1; j >= 0; j--)
   {
-    m_eivec.col(j).segment(0, size) = m_eivec.leftCols(j+1) * matH.col(j).segment(0, j+1);
+    m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
+    m_eivec.col(j) = m_tmp;
   }
 }
 
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 00a9d368c..25b18dd8d 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -110,9 +111,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
       */
     SelfAdjointEigenSolver()
-        : m_eivec(int(Size), int(Size)),
-          m_eivalues(int(Size)),
-          m_subdiag(int(TridiagonalizationType::SizeMinusOne))
+        : m_eivec(),
+          m_eivalues(),
+          m_tridiag(),
+          m_subdiag()
     {
       ei_assert(Size!=Dynamic);
     }
@@ -133,7 +135,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     SelfAdjointEigenSolver(int size)
         : m_eivec(size, size),
           m_eivalues(size),
-          m_subdiag(TridiagonalizationType::SizeMinusOne)
+          m_tridiag(size),
+          m_subdiag(size > 1 ? size - 1 : 1)
     {}
 
     /** \brief Constructor; computes eigendecomposition of given matrix. 
@@ -157,9 +160,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
       : m_eivec(matrix.rows(), matrix.cols()),
         m_eivalues(matrix.cols()),
-        m_subdiag()
+        m_tridiag(matrix.rows()),
+        m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1)
     {
-      if (matrix.rows() > 1) m_subdiag.resize(matrix.rows() - 1);
       compute(matrix, computeEigenvectors);
     }
 
@@ -187,9 +190,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
       : m_eivec(matA.rows(), matA.cols()),
         m_eivalues(matA.cols()),
-        m_subdiag()
+        m_tridiag(matA.rows()),
+        m_subdiag(matA.rows() > 1 ? matA.rows() - 1 : 1)
     {
-      if (matA.rows() > 1) m_subdiag.resize(matA.rows() - 1);
       compute(matA, matB, computeEigenvectors);
     }
 
@@ -351,6 +354,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
   protected:
     MatrixType m_eivec;
     RealVectorType m_eivalues;
+    TridiagonalizationType m_tridiag;
     typename TridiagonalizationType::SubDiagonalType m_subdiag;
     #ifndef NDEBUG
     bool m_eigenvectorsOk;
@@ -396,14 +400,12 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
     return *this;
   }
 
-  m_eivec = matrix;
-
-  // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?
-  // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
-  // (same for diag and subdiag)
+  m_tridiag.compute(matrix);
   RealVectorType& diag = m_eivalues;
-  m_subdiag.resize(n-1);
-  TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag, computeEigenvectors);
+  diag = m_tridiag.diagonal();
+  m_subdiag = m_tridiag.subDiagonal();
+  if (computeEigenvectors)
+    m_eivec = m_tridiag.matrixQ();
 
   int end = n-1;
   int start = 0;
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index 43509863a..6ea852a6b 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -70,10 +70,10 @@ template<typename _MatrixType> class Tridiagonalization
 
     enum {
       Size = MatrixType::RowsAtCompileTime,
-      SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
+      SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
       Options = MatrixType::Options,
       MaxSize = MatrixType::MaxRowsAtCompileTime,
-      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
+      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
     };
 
     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
@@ -108,7 +108,7 @@ template<typename _MatrixType> class Tridiagonalization
       * \sa compute() for an example.
       */
     Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
-      : m_matrix(size,size), m_hCoeffs(size-1)
+      : m_matrix(size,size), m_hCoeffs(size > 1 ? size-1 : 1)
     {}
 
     /** \brief Constructor; computes tridiagonal decomposition of given matrix. 
@@ -122,7 +122,7 @@ template<typename _MatrixType> class Tridiagonalization
       * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
       */
     Tridiagonalization(const MatrixType& matrix)
-      : m_matrix(matrix), m_hCoeffs(matrix.cols()-1)
+      : m_matrix(matrix), m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1)
     {
       _compute(m_matrix, m_hCoeffs);
     }