diff --git a/Data/Eigen/Internal.hs b/Data/Eigen/Internal.hs
--- a/Data/Eigen/Internal.hs
+++ b/Data/Eigen/Internal.hs
@@ -10,10 +10,12 @@
 import System.IO.Unsafe
 import Data.Complex
 
-class (Num a, Ord a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where
+class (Num a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where
 
 instance Elem Float CFloat where
 instance Elem Double CDouble where
+instance Elem (Complex Float) (CComplex CFloat) where
+instance Elem (Complex Double) (CComplex CDouble) where
 
 class Cast a b where
     cast :: a -> b
@@ -57,6 +59,7 @@
 foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt
 
 foreign import ccall "eigen-proxy.h eigen_random"      c_random      :: CInt -> Ptr a -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_identity"    c_identity    :: CInt -> Ptr a -> CInt -> CInt -> IO CString
 foreign import ccall "eigen-proxy.h eigen_add"         c_add         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
 foreign import ccall "eigen-proxy.h eigen_sub"         c_sub         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
 foreign import ccall "eigen-proxy.h eigen_mul"         c_mul         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
@@ -90,6 +93,7 @@
 instance Code (CComplex CDouble) where; code _ = 3
 
 random      :: forall a . Code a => Ptr a -> CInt -> CInt -> IO CString
+identity    :: forall a . Code a => Ptr a -> CInt -> CInt -> IO CString
 add         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
 sub         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
 mul         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
@@ -116,6 +120,7 @@
 relativeError :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString
 
 random      = c_random      (code (undefined :: a))
+identity    = c_identity    (code (undefined :: a))
 add         = c_add         (code (undefined :: a))
 sub         = c_sub         (code (undefined :: a))
 mul         = c_mul         (code (undefined :: a))
diff --git a/Data/Eigen/Matrix.hs b/Data/Eigen/Matrix.hs
--- a/Data/Eigen/Matrix.hs
+++ b/Data/Eigen/Matrix.hs
@@ -159,13 +159,12 @@
 ones :: I.Elem a b => Int -> Int -> Matrix a b
 ones rows cols = constant rows cols 1
 
--- | Square matrix with 1 on main diagonal and 0 elsewhere
-identity :: forall a b . I.Elem a b => Int -> Matrix a b
-identity size = Matrix size size $ VS.create $ do
-    vm <- VSM.replicate (size * size) (I.cast (0::a))
-    forM_ [0..pred size] $ \n ->
-        VSM.write vm (n * size + n) (I.cast (1::a))
-    return vm
+-- | The identity matrix (not necessarily square).
+identity :: I.Elem a b => Int -> Int -> Matrix a b
+identity rows cols = I.performIO $ do
+    m <- M.new rows cols
+    I.call $ M.unsafeWith m I.identity
+    unsafeFreeze m
 
 -- | The random matrix of a given size
 random :: I.Elem a b => Int -> Int -> IO (Matrix a b)
@@ -224,11 +223,11 @@
 valid (Matrix rows cols vals) = rows >= 0 && cols >= 0 && VS.length vals == rows * cols
 
 -- | The maximum coefficient of the matrix
-maxCoeff :: I.Elem a b => Matrix a b -> a
+maxCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a
 maxCoeff = fold1' max
 
 -- | The minimum coefficient of the matrix
-minCoeff :: I.Elem a b => Matrix a b -> a
+minCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a
 minCoeff = fold1' min
 
 -- | Top @N@ rows of matrix
@@ -248,12 +247,12 @@
 rightCols n m@(Matrix rows cols _) = block 0 (cols - n) rows n m
 
 -- | Construct matrix from a list of rows, column count is detected as maximum row length. Missing values are filled with 0
-fromList :: forall a b . I.Elem a b => [[a]] -> Matrix a b
+fromList :: I.Elem a b => [[a]] -> Matrix a b
 fromList list = Matrix rows cols vals where
     rows = length list
     cols = maximum $ P.map length list
     vals = VS.create $ do
-        vm <- VSM.replicate (rows * cols) (I.cast (0::a))
+        vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head (head list))))
         forM_ (zip [0..] list) $ \(row, vals) ->
             forM_ (zip [0..] vals) $ \(col, val) ->
                 VSM.write vm (col * rows + row) (I.cast val)
diff --git a/cbits/eigen-proxy.cpp b/cbits/eigen-proxy.cpp
--- a/cbits/eigen-proxy.cpp
+++ b/cbits/eigen-proxy.cpp
@@ -120,6 +120,18 @@
     GUARD_END
 }
 
+extern "C" const char* eigen_identity(int code, void* p, int r, int c)
+{
+    GUARD_START
+    switch (code) {
+        case 0: matrix<T0>(p,r,c) = MatrixXf::Identity(r,c); break;
+        case 1: matrix<T1>(p,r,c) = MatrixXd::Identity(r,c); break;
+        case 2: matrix<T2>(p,r,c) = MatrixXcf::Identity(r,c); break;
+        case 3: matrix<T3>(p,r,c) = MatrixXcd::Identity(r,c); break;
+    }
+    GUARD_END
+}
+
 template <class T>
 const char* rank(Decomposition d, int* v, const void* p, int r, int c) {
     typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;
diff --git a/cbits/eigen-proxy.h b/cbits/eigen-proxy.h
--- a/cbits/eigen-proxy.h
+++ b/cbits/eigen-proxy.h
@@ -28,6 +28,7 @@
 const char* eigen_transpose(int, void*, int, int, const void*, int, int);
 const char* eigen_normalize(int, void*, int, int);
 const char* eigen_random(int, void*, int, int);
+const char* eigen_identity(int, void*, int, int);
 const char* eigen_inverse(int, void*, int, int, const void*, int, int);
 const char* eigen_conjugate(int, void*, int, int, const void*, int, int);
 const char* eigen_adjoint(int, void*, int, int, const void*, int, int);
diff --git a/eigen.cabal b/eigen.cabal
--- a/eigen.cabal
+++ b/eigen.cabal
@@ -1,5 +1,5 @@
 name:           eigen
-version:        2.0.0
+version:        2.0.1
 homepage:       https://github.com/osidorkin/haskell-eigen
 synopsis:       Eigen C++ library (linear algebra: matrices, vectors, numerical solvers).
 description:    This module provides Haskell binding for Eigen C++ library.