diff --git a/src/cppmat/fix_diagonal_matrix.h b/src/cppmat/fix_diagonal_matrix.h
index 06bce76..124a036 100644
--- a/src/cppmat/fix_diagonal_matrix.h
+++ b/src/cppmat/fix_diagonal_matrix.h
@@ -211,6 +211,18 @@ class matrix
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
 // =================================================================================================
 // external arithmetic operators (cppmat::tiny::diagonal)
 // =================================================================================================
diff --git a/src/cppmat/fix_diagonal_matrix.hpp b/src/cppmat/fix_diagonal_matrix.hpp
index bea860e..3e09434 100644
--- a/src/cppmat/fix_diagonal_matrix.hpp
+++ b/src/cppmat/fix_diagonal_matrix.hpp
@@ -1235,6 +1235,34 @@ std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+inline
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+inline
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================
diff --git a/src/cppmat/fix_regular_array.h b/src/cppmat/fix_regular_array.h
index 3d47d79..f9bfe7f 100644
--- a/src/cppmat/fix_regular_array.h
+++ b/src/cppmat/fix_regular_array.h
@@ -317,6 +317,18 @@ class array
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+bool operator!= (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+bool operator== (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B);
+
 // =================================================================================================
 // external arithmetic operators
 // =================================================================================================
diff --git a/src/cppmat/fix_regular_array.hpp b/src/cppmat/fix_regular_array.hpp
index c87dbbd..c12dc67 100644
--- a/src/cppmat/fix_regular_array.hpp
+++ b/src/cppmat/fix_regular_array.hpp
@@ -2458,6 +2458,34 @@ std::ostream& operator<<(std::ostream& out, const array<X,RANK,I,J,K,L,M,N>& src
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+inline
+bool operator!= (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+inline
+bool operator== (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================
diff --git a/src/cppmat/fix_symmetric_matrix.h b/src/cppmat/fix_symmetric_matrix.h
index 9aad37e..eb74939 100644
--- a/src/cppmat/fix_symmetric_matrix.h
+++ b/src/cppmat/fix_symmetric_matrix.h
@@ -215,6 +215,18 @@ class matrix
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
 // =================================================================================================
 // external arithmetic operators (cppmat::tiny::symmetric)
 // =================================================================================================
diff --git a/src/cppmat/fix_symmetric_matrix.hpp b/src/cppmat/fix_symmetric_matrix.hpp
index 22b645b..da431d1 100644
--- a/src/cppmat/fix_symmetric_matrix.hpp
+++ b/src/cppmat/fix_symmetric_matrix.hpp
@@ -1299,6 +1299,34 @@ std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+inline
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+inline
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================
diff --git a/src/cppmat/map_diagonal_matrix.h b/src/cppmat/map_diagonal_matrix.h
index 877953f..5c820c5 100644
--- a/src/cppmat/map_diagonal_matrix.h
+++ b/src/cppmat/map_diagonal_matrix.h
@@ -135,7 +135,22 @@ class matrix
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// =================================================================================================
 // print operator
+// =================================================================================================
+
 template<class X, size_t M, size_t N>
 std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src);
 
diff --git a/src/cppmat/map_diagonal_matrix.hpp b/src/cppmat/map_diagonal_matrix.hpp
index 0b9f6fb..f17ba90 100644
--- a/src/cppmat/map_diagonal_matrix.hpp
+++ b/src/cppmat/map_diagonal_matrix.hpp
@@ -546,6 +546,34 @@ std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+inline
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+inline
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 
 }}} // namespace ...
diff --git a/src/cppmat/map_regular_array.h b/src/cppmat/map_regular_array.h
index 6229d83..9790e66 100644
--- a/src/cppmat/map_regular_array.h
+++ b/src/cppmat/map_regular_array.h
@@ -198,7 +198,22 @@ class array
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+bool operator!= (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+bool operator== (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B);
+
+// =================================================================================================
 // print operator
+// =================================================================================================
+
 template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
 std::ostream& operator<<(std::ostream& out, const array<X,RANK,I,J,K,L,M,N>& src);
 
diff --git a/src/cppmat/map_regular_array.hpp b/src/cppmat/map_regular_array.hpp
index dfc94d1..ca4c227 100644
--- a/src/cppmat/map_regular_array.hpp
+++ b/src/cppmat/map_regular_array.hpp
@@ -1290,6 +1290,34 @@ std::ostream& operator<<(std::ostream& out, const array<X,RANK,I,J,K,L,M,N>& src
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+inline
+bool operator!= (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t RANK, size_t I, size_t J, size_t K, size_t L, size_t M, size_t N>
+inline
+bool operator== (const array<X,RANK,I,J,K,L,M,N> &A, const array<X,RANK,I,J,K,L,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 
 }} // namespace ...
diff --git a/src/cppmat/map_symmetric_matrix.h b/src/cppmat/map_symmetric_matrix.h
index 931b035..b3c5816 100644
--- a/src/cppmat/map_symmetric_matrix.h
+++ b/src/cppmat/map_symmetric_matrix.h
@@ -134,7 +134,22 @@ class matrix
 
 };
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B);
+
+// =================================================================================================
 // print operator
+// =================================================================================================
+
 template<class X, size_t M, size_t N>
 std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src);
 
diff --git a/src/cppmat/map_symmetric_matrix.hpp b/src/cppmat/map_symmetric_matrix.hpp
index 8786dc4..cb583d6 100644
--- a/src/cppmat/map_symmetric_matrix.hpp
+++ b/src/cppmat/map_symmetric_matrix.hpp
@@ -556,6 +556,34 @@ std::ostream& operator<<(std::ostream& out, const matrix<X,M,N>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X, size_t M, size_t N>
+inline
+bool operator!= (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X, size_t M, size_t N>
+inline
+bool operator== (const matrix<X,M,N> &A, const matrix<X,M,N> &B)
+{
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 
 }}} // namespace ...
diff --git a/src/cppmat/var_diagonal_matrix.h b/src/cppmat/var_diagonal_matrix.h
index 7010ae2..5b19e48 100644
--- a/src/cppmat/var_diagonal_matrix.h
+++ b/src/cppmat/var_diagonal_matrix.h
@@ -213,6 +213,10 @@ class matrix
 
 };
 
+// equality operators
+template<class X> bool operator!= (const matrix<X> &A, const matrix<X> &B);
+template<class X> bool operator== (const matrix<X> &A, const matrix<X> &B);
+
 // external arithmetic operators (cppmat::diagonal)
 template<class X> matrix<X> operator* (const matrix<X> &A, const matrix<X> &B);
 template<class X> matrix<X> operator+ (const matrix<X> &A, const matrix<X> &B);
diff --git a/src/cppmat/var_diagonal_matrix.hpp b/src/cppmat/var_diagonal_matrix.hpp
index 5df6e9a..d6241cb 100644
--- a/src/cppmat/var_diagonal_matrix.hpp
+++ b/src/cppmat/var_diagonal_matrix.hpp
@@ -1303,6 +1303,42 @@ std::ostream& operator<<(std::ostream& out, const matrix<X>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X>
+inline
+bool operator!= (const matrix<X> &A, const matrix<X> &B)
+{
+  assert( A.shape() == B.shape() );
+  assert( A.rank () == B.rank () );
+  assert( A.size () == B.size () );
+
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X>
+inline
+bool operator== (const matrix<X> &A, const matrix<X> &B)
+{
+  assert( A.shape() == B.shape() );
+  assert( A.rank () == B.rank () );
+  assert( A.size () == B.size () );
+
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================
diff --git a/src/cppmat/var_regular_array.h b/src/cppmat/var_regular_array.h
index ac4914d..0e906df 100644
--- a/src/cppmat/var_regular_array.h
+++ b/src/cppmat/var_regular_array.h
@@ -349,6 +349,7 @@ class array
 
 // equality operators
 template<class X> bool operator!= (const array<X> &A, const array<X> &B);
+template<class X> bool operator== (const array<X> &A, const array<X> &B);
 
 // external arithmetic operators
 template<class X> array<X> operator* (const array<X> &A, const array<X> &B);
diff --git a/src/cppmat/var_regular_array.hpp b/src/cppmat/var_regular_array.hpp
index 996201f..159ea76 100644
--- a/src/cppmat/var_regular_array.hpp
+++ b/src/cppmat/var_regular_array.hpp
@@ -3082,7 +3082,6 @@ std::ostream& operator<<(std::ostream& out, const array<X>& src)
 // equality operators
 // =================================================================================================
 
-// equality operators
 template<class X>
 inline
 bool operator!= (const array<X> &A, const array<X> &B)
@@ -3098,6 +3097,23 @@ bool operator!= (const array<X> &A, const array<X> &B)
   return false;
 }
 
+// -------------------------------------------------------------------------------------------------
+
+template<class X>
+inline
+bool operator== (const array<X> &A, const array<X> &B)
+{
+  assert( A.shape() == B.shape() );
+  assert( A.rank () == B.rank () );
+  assert( A.size () == B.size () );
+
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================
diff --git a/src/cppmat/var_symmetric_matrix.h b/src/cppmat/var_symmetric_matrix.h
index ccfa8e4..e314f0d 100644
--- a/src/cppmat/var_symmetric_matrix.h
+++ b/src/cppmat/var_symmetric_matrix.h
@@ -217,6 +217,10 @@ class matrix
 
 };
 
+// equality operators
+template<class X> bool operator!= (const matrix<X> &A, const matrix<X> &B);
+template<class X> bool operator== (const matrix<X> &A, const matrix<X> &B);
+
 // external arithmetic operators (cppmat::symmetric::matrix)
 template<class X> matrix<X> operator* (const matrix<X> &A, const matrix<X> &B);
 template<class X> matrix<X> operator/ (const matrix<X> &A, const matrix<X> &B);
diff --git a/src/cppmat/var_symmetric_matrix.hpp b/src/cppmat/var_symmetric_matrix.hpp
index bd53700..da7e5aa 100644
--- a/src/cppmat/var_symmetric_matrix.hpp
+++ b/src/cppmat/var_symmetric_matrix.hpp
@@ -1359,6 +1359,42 @@ std::ostream& operator<<(std::ostream& out, const matrix<X>& src)
   return out;
 }
 
+// =================================================================================================
+// equality operators
+// =================================================================================================
+
+template<class X>
+inline
+bool operator!= (const matrix<X> &A, const matrix<X> &B)
+{
+  assert( A.shape() == B.shape() );
+  assert( A.rank () == B.rank () );
+  assert( A.size () == B.size () );
+
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return true;
+
+  return false;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+template<class X>
+inline
+bool operator== (const matrix<X> &A, const matrix<X> &B)
+{
+  assert( A.shape() == B.shape() );
+  assert( A.rank () == B.rank () );
+  assert( A.size () == B.size () );
+
+  for ( size_t i = 0 ; i < A.size() ; ++i )
+    if ( A[i] != B[i] )
+      return false;
+
+  return true;
+}
+
 // =================================================================================================
 // arithmetic operators: external
 // =================================================================================================