rework Identity API: no longer restricted to square matrices

Eigen/src/Core/Identity.h

@@ -30,7 +30,7 @@
   *
   * \brief Expression of the identity matrix of some size.
   *
-  * \sa MatrixBase::identity(int)
+  * \sa MatrixBase::identity(), MatrixBase::identity(int,int), MatrixBase::setIdentity()
   */
 template<typename MatrixType> class Identity : NoOperatorEquals,
   public MatrixBase<typename MatrixType::Scalar, Identity<MatrixType> >
@@ -39,9 +39,12 @@ template<typename MatrixType> class Identity : NoOperatorEquals,
     typedef typename MatrixType::Scalar Scalar;
     friend class MatrixBase<Scalar, Identity<MatrixType> >;
     
-    Identity(int rows) : m_rows(rows)
+    Identity(int rows, int cols) : m_rows(rows), m_cols(cols)
     {
-      assert(rows > 0 && RowsAtCompileTime == ColsAtCompileTime);
+      assert(rows > 0
+          && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
+          && cols > 0
+          && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
     }
     
   private:
@@ -52,7 +55,7 @@ template<typename MatrixType> class Identity : NoOperatorEquals,
     
     const Identity& _ref() const { return *this; }
     int _rows() const { return m_rows; }
-    int _cols() const { return m_rows; }
+    int _cols() const { return m_cols; }
     
     Scalar _coeff(int row, int col) const
     {
@@ -60,65 +63,84 @@ template<typename MatrixType> class Identity : NoOperatorEquals,
     }
     
   protected:
-    const int m_rows;
+    const int m_rows, m_cols;
 };
 
-/** \returns an expression of the identity matrix of given type and size.
+/** \returns an expression of the identity matrix (not necessarily square).
   *
-  * \param rows The number of rows of the identity matrix to return. If *this has
+  * The parameters \a rows and \a cols are the number of rows and of columns of
-  *             fixed size, that size is used as the default argument for \a rows
+  * the returned matrix. Must be compatible with this MatrixBase type.
-  *             and is then the only allowed value. If *this has dynamic size,
-  *             you can use any positive value for \a rows.
   *
-  * \note An identity matrix is a square matrix, so it is required that the type of *this
+  * This variant is meant to be used for dynamic-size matrix types. For fixed-size types,
-  *       allows being a square matrix.
+  * it is redundant to pass \a rows and \a cols as arguments, so identity() should be used
+  * instead.
   *
-  * Example: \include MatrixBase_identity_int.cpp
+  * Example: \include MatrixBase_identity_int_int.cpp
-  * Output: \verbinclude MatrixBase_identity_int.out
+  * Output: \verbinclude MatrixBase_identity_int_int.out
   *
-  * \sa class Identity, isIdentity()
+  * \sa identity(), setIdentity(), isIdentity()
   */
 template<typename Scalar, typename Derived>
-const Identity<Derived> MatrixBase<Scalar, Derived>::identity(int rows)
+const Identity<Derived> MatrixBase<Scalar, Derived>::identity(int rows, int cols)
 {
-  return Identity<Derived>(rows);
+  return Identity<Derived>(rows, cols);
 }
 
-/** \returns true if *this is approximately equal to the identity matrix,
+/** \returns an expression of the identity matrix (not necessarily square).
+  *
+  * This variant is only for fixed-size MatrixBase types. For dynamic-size types, you
+  * need to use the variant taking size arguments.
+  *
+  * Example: \include MatrixBase_identity.cpp
+  * Output: \verbinclude MatrixBase_identity.out
+  *
+  * \sa identity(int,int), setIdentity(), isIdentity()
+  */
+template<typename Scalar, typename Derived>
+const Identity<Derived> MatrixBase<Scalar, Derived>::identity()
+{
+  return Identity<Derived>(Traits::RowsAtCompileTime, Traits::ColsAtCompileTime);
+}
+
+/** \returns true if *this is approximately equal to the identity matrix
+  *          (not necessarily square),
   *          within the precision given by \a prec.
   *
   * Example: \include MatrixBase_isIdentity.cpp
   * Output: \verbinclude MatrixBase_isIdentity.out
   *
-  * \sa class Identity, identity(int)
+  * \sa class Identity, identity(), identity(int,int), setIdentity()
   */
 template<typename Scalar, typename Derived>
 bool MatrixBase<Scalar, Derived>::isIdentity
 (typename NumTraits<Scalar>::Real prec) const
 {
-  if(cols() != rows()) return false;
   for(int j = 0; j < cols(); j++)
   {
-    if(!Eigen::isApprox(coeff(j, j), static_cast<Scalar>(1), prec))
+    for(int i = 0; i < rows(); i++)
-      return false;
+    {
-    for(int i = 0; i < j; i++)
+      if(i == j)
-      if(!Eigen::isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1), prec))
+        if(!Eigen::isApprox(coeff(i, j), static_cast<Scalar>(1), prec))
-        return false;
+          return false;
+      else
+        if(!Eigen::isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1), prec))
+          return false;
+    }
   }
   return true;
 }
 
-/** Writes the identity expression into *this.
+/** Writes the identity expression (not necessarily square) into *this.
   *
   * Example: \include MatrixBase_setIdentity.cpp
   * Output: \verbinclude MatrixBase_setIdentity.out
   *
-  * \sa class Identity, identity()
+  * \sa class Identity, identity(), identity(int,int), isIdentity()
   */
 template<typename Scalar, typename Derived>
 Derived& MatrixBase<Scalar, Derived>::setIdentity()
 {
-  return *this = Identity<Derived>(rows());
+  return *this = Identity<Derived>(rows(), cols());
 }
 
 

Eigen/src/Core/MatrixBase.h

@@ -182,7 +182,8 @@ template<typename Scalar, typename Derived> class MatrixBase
     static const Ones<Derived> ones(int rows, int cols);
     static const Ones<Derived> ones(int size);
     static const Ones<Derived> ones();
-    static const Identity<Derived> identity(int rows = Derived::RowsAtCompileTime);
+    static const Identity<Derived> identity();
+    static const Identity<Derived> identity(int rows, int cols);
     
     bool isZero(RealScalar prec = precision<Scalar>()) const;
     bool isOnes(RealScalar prec = precision<Scalar>()) const;

Eigen/src/Core/Ones.h

@@ -30,7 +30,8 @@
   *
   * \brief Expression of a matrix where all coefficients equal one.
   *
-  * \sa MatrixBase::ones(), MatrixBase::ones(int), MatrixBase::ones(int,int)
+  * \sa MatrixBase::ones(), MatrixBase::ones(int), MatrixBase::ones(int,int),
+  *     MatrixBase::setOnes(), MatrixBase::isOnes()
   */
 template<typename MatrixType> class Ones : NoOperatorEquals,
   public MatrixBase<typename MatrixType::Scalar, Ones<MatrixType> >

Eigen/src/Core/Random.h

@@ -30,7 +30,8 @@
   *
   * \brief Expression of a random matrix or vector.
   *
-  * \sa MatrixBase::random(), MatrixBase::random(int), MatrixBase::random(int,int)
+  * \sa MatrixBase::random(), MatrixBase::random(int), MatrixBase::random(int,int),
+  *     MatrixBase::setRandom()
   */
 template<typename MatrixType> class Random : NoOperatorEquals,
   public MatrixBase<typename MatrixType::Scalar, Random<MatrixType> >

Eigen/src/Core/Zero.h

@@ -30,7 +30,8 @@
   *
   * \brief Expression of a zero matrix or vector.
   *
-  * \sa MatrixBase::zero(), MatrixBase::zero(int), MatrixBase::zero(int,int)
+  * \sa MatrixBase::zero(), MatrixBase::zero(int), MatrixBase::zero(int,int),
+  *     MatrixBase::setZero(), MatrixBase::isZero()
   */
 template<typename MatrixType> class Zero : NoOperatorEquals,
   public MatrixBase<typename MatrixType::Scalar, Zero<MatrixType> >

doc/snippets/MatrixBase_identity.cpp

@@ -0,0 +1 @@
+cout << Matrix<double, 3, 4>::identity() << endl;

doc/snippets/MatrixBase_identity_int.cpp

@@ -1,2 +0,0 @@
-cout << Matrix2d::identity() << endl;
-cout << MatrixXd::identity(3) << endl;

doc/snippets/MatrixBase_identity_int_int.cpp

@@ -0,0 +1 @@
+cout << MatrixXd::identity(4, 3) << endl;

test/adjoint.cpp

@@ -43,7 +43,7 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
              m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
-                              ::identity(rows),
+                              ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
                               ::random(rows, rows);
   VectorType v1 = VectorType::random(rows),

test/basicstuff.cpp

@@ -42,7 +42,7 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
              m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
-                              ::identity(rows),
+                              ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
                               ::random(rows, rows);
   VectorType v1 = VectorType::random(rows),

test/linearstructure.cpp

@@ -46,7 +46,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
              m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
-                              ::identity(rows),
+                              ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
                               ::random(rows, rows);
   VectorType v1 = VectorType::random(rows),

test/miscmatrices.cpp

@@ -51,7 +51,7 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m)
   else VERIFY_IS_MUCH_SMALLER_THAN(square(r,r2), static_cast<Scalar>(1));
   square = MatrixType::zero(rows, rows);
   square.diagonal() = VectorType::ones(rows);
-  VERIFY_IS_APPROX(square, MatrixType::identity(rows));
+  VERIFY_IS_APPROX(square, MatrixType::identity(rows, rows));
 }
 
 void EigenTest::testMiscMatrices()

test/product.cpp

@@ -46,7 +46,7 @@ template<typename MatrixType> void product(const MatrixType& m)
              m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
-                              ::identity(rows),
+                              ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
                               ::random(rows, rows);
   VectorType v1 = VectorType::random(rows),
@@ -83,7 +83,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   VERIFY_IS_APPROX(m1,                      identity*m1);
   VERIFY_IS_APPROX(v1,                      identity*v1);
   // again, test operator() to check const-qualification
-  VERIFY_IS_APPROX(MatrixType::identity(std::max(rows,cols))(r,c), static_cast<Scalar>(r==c));
+  VERIFY_IS_APPROX(MatrixType::identity(rows, cols)(r,c), static_cast<Scalar>(r==c));
 }
 
 void EigenTest::testProduct()

test/submatrices.cpp

@@ -44,7 +44,7 @@ template<typename MatrixType> void submatrices(const MatrixType& m)
              m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
-                              ::identity(rows),
+                              ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::Traits::RowsAtCompileTime, MatrixType::Traits::RowsAtCompileTime>
                               ::random(rows, rows);
   VectorType v1 = VectorType::random(rows),