diff --git a/src/Core/Cast.h b/src/Core/Cast.h
index 50fa72697..5dc33092a 100644
--- a/src/Core/Cast.h
+++ b/src/Core/Cast.h
@@ -58,7 +58,7 @@ template<typename NewScalar, typename MatrixType> class Cast : NoDefaultOperator
 
 template<typename Scalar, typename Derived>
 template<typename NewScalar>
-Cast<NewScalar, Derived>
+const Cast<NewScalar, Derived>
 MatrixBase<Scalar, Derived>::cast() const
 {
   return Cast<NewScalar, Derived>(static_cast<const Derived*>(this)->ref());
diff --git a/src/Core/Conjugate.h b/src/Core/Conjugate.h
index a529ce630..04ba143c0 100644
--- a/src/Core/Conjugate.h
+++ b/src/Core/Conjugate.h
@@ -57,7 +57,7 @@ template<typename MatrixType> class Conjugate : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived>
-Conjugate<Derived>
+const Conjugate<Derived>
 MatrixBase<Scalar, Derived>::conjugate() const
 {
   return Conjugate<Derived>(static_cast<const Derived*>(this)->ref());
diff --git a/src/Core/Difference.h b/src/Core/Difference.h
index 8bebfb9df..94412f79e 100644
--- a/src/Core/Difference.h
+++ b/src/Core/Difference.h
@@ -63,7 +63,7 @@ template<typename Lhs, typename Rhs> class Difference : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived1, typename Derived2>
-Difference<Derived1, Derived2>
+const Difference<Derived1, Derived2>
 operator-(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
 {
   return Difference<Derived1, Derived2>(mat1.ref(), mat2.ref());
diff --git a/src/Core/Identity.h b/src/Core/Identity.h
index a5bebbca9..067626f22 100644
--- a/src/Core/Identity.h
+++ b/src/Core/Identity.h
@@ -58,7 +58,7 @@ template<typename MatrixType> class Identity : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived>
-Identity<Derived> MatrixBase<Scalar, Derived>::identity(int rows)
+const Identity<Derived> MatrixBase<Scalar, Derived>::identity(int rows)
 {
   return Identity<Derived>(rows);
 }
diff --git a/src/Core/MatrixBase.h b/src/Core/MatrixBase.h
index bdf1a1c0b..12be45420 100644
--- a/src/Core/MatrixBase.h
+++ b/src/Core/MatrixBase.h
@@ -57,7 +57,7 @@ template<typename Scalar, typename Derived> class MatrixBase
       return this->operator=<Derived>(other);
     }
     
-    template<typename NewScalar> Cast<NewScalar, Derived> cast() const;
+    template<typename NewScalar> const Cast<NewScalar, Derived> cast() const;
     
     Row<Derived> row(int i) const;
     Column<Derived> col(int i) const;
@@ -69,8 +69,9 @@ template<typename Scalar, typename Derived> class MatrixBase
     Block<Derived, BlockRows, BlockCols> block(int startRow, int startCol) const;
     
     Transpose<Derived> transpose() const;
-    Conjugate<Derived> conjugate() const;
-    Transpose<Conjugate<Derived> > adjoint() const { return conjugate().transpose(); }
+    const Conjugate<Derived> conjugate() const;
+    const Transpose<Conjugate<Derived> > adjoint() const
+    { return conjugate().transpose(); }
     Scalar trace() const;
     
     template<typename OtherDerived>
@@ -81,9 +82,9 @@ template<typename Scalar, typename Derived> class MatrixBase
     
     static Eval<Random<Derived> >
     random(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime);
-    static Zero<Derived>
+    static const Zero<Derived>
     zero(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime);
-    static Identity<Derived>
+    static const Identity<Derived>
     identity(int rows = RowsAtCompileTime);
     static FromArray<Derived>
     fromArray(const Scalar* array, int rows = RowsAtCompileTime, int cols = ColsAtCompileTime);
@@ -107,7 +108,7 @@ template<typename Scalar, typename Derived> class MatrixBase
     const Product<Derived, OtherDerived>
     lazyProduct(const MatrixBase<Scalar, OtherDerived>& other) const EIGEN_ALWAYS_INLINE;
     
-    Opposite<Derived> operator-() const;
+    const Opposite<Derived> operator-() const;
     
     template<typename OtherDerived>
     Derived& operator+=(const MatrixBase<Scalar, OtherDerived>& other);
diff --git a/src/Core/Opposite.h b/src/Core/Opposite.h
index 9b6cb3093..93ea645f8 100644
--- a/src/Core/Opposite.h
+++ b/src/Core/Opposite.h
@@ -57,7 +57,7 @@ template<typename MatrixType> class Opposite : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived>
-Opposite<Derived>
+const Opposite<Derived>
 MatrixBase<Scalar, Derived>::operator-() const
 {
   return Opposite<Derived>(static_cast<const Derived*>(this)->ref());
diff --git a/src/Core/ScalarMultiple.h b/src/Core/ScalarMultiple.h
index e644769c2..268b5fa41 100644
--- a/src/Core/ScalarMultiple.h
+++ b/src/Core/ScalarMultiple.h
@@ -58,42 +58,42 @@ template<typename MatrixType> class ScalarMultiple : NoDefaultOperatorEquals,
     const Scalar m_scalar;
 };
 
-#define EIGEN_MAKE_SCALAR_OPS(OtherScalar)                                \
+#define EIGEN_MAKE_SCALAR_OPS(OtherScalar)                             \
 template<typename Scalar, typename Derived>                            \
-ScalarMultiple<Derived>                                               \
-operator*(const MatrixBase<Scalar, Derived>& matrix,                     \
+const ScalarMultiple<Derived>                                          \
+operator*(const MatrixBase<Scalar, Derived>& matrix,                   \
           OtherScalar scalar)                                          \
 {                                                                      \
-  return ScalarMultiple<Derived>(matrix.ref(), scalar);          \
+  return ScalarMultiple<Derived>(matrix.ref(), scalar);                \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
-ScalarMultiple<Derived>                                               \
+const ScalarMultiple<Derived>                                          \
 operator*(OtherScalar scalar,                                          \
-          const MatrixBase<Scalar, Derived>& matrix)                     \
+          const MatrixBase<Scalar, Derived>& matrix)                   \
 {                                                                      \
-  return ScalarMultiple<Derived>(matrix.ref(), scalar);          \
+  return ScalarMultiple<Derived>(matrix.ref(), scalar);                \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
-ScalarMultiple<Derived>                                               \
-operator/(const MatrixBase<Scalar, Derived>& matrix,                     \
+const ScalarMultiple<Derived>                                          \
+operator/(const MatrixBase<Scalar, Derived>& matrix,                   \
           OtherScalar scalar)                                          \
 {                                                                      \
-  assert(NumTraits<Scalar>::HasFloatingPoint);                       \
+  assert(NumTraits<Scalar>::HasFloatingPoint);                         \
   return matrix * (static_cast<Scalar>(1) / scalar);                   \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
 Derived &                                                              \
-MatrixBase<Scalar, Derived>::operator*=(const OtherScalar &other)        \
+MatrixBase<Scalar, Derived>::operator*=(const OtherScalar &other)      \
 {                                                                      \
   return *this = *this * other;                                        \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
 Derived &                                                              \
-MatrixBase<Scalar, Derived>::operator/=(const OtherScalar &other)        \
+MatrixBase<Scalar, Derived>::operator/=(const OtherScalar &other)      \
 {                                                                      \
   return *this = *this / other;                                        \
 }
diff --git a/src/Core/Sum.h b/src/Core/Sum.h
index e110b3771..fa2b206da 100644
--- a/src/Core/Sum.h
+++ b/src/Core/Sum.h
@@ -62,7 +62,7 @@ template<typename Lhs, typename Rhs> class Sum : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived1, typename Derived2>
-Sum<Derived1, Derived2>
+const Sum<Derived1, Derived2>
 operator+(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
 {
   return Sum<Derived1, Derived2>(mat1.ref(), mat2.ref());
diff --git a/src/Core/Zero.h b/src/Core/Zero.h
index 8cb98ad82..9342a66c8 100644
--- a/src/Core/Zero.h
+++ b/src/Core/Zero.h
@@ -56,7 +56,7 @@ template<typename MatrixType> class Zero : NoDefaultOperatorEquals,
 };
 
 template<typename Scalar, typename Derived>
-Zero<Derived> MatrixBase<Scalar, Derived>::zero(int rows, int cols)
+const Zero<Derived> MatrixBase<Scalar, Derived>::zero(int rows, int cols)
 {
   return Zero<Derived>(rows, cols);
 }
diff --git a/test/adjoint.cpp b/test/adjoint.cpp
index 1d98d6a9b..7a33f8b73 100644
--- a/test/adjoint.cpp
+++ b/test/adjoint.cpp
@@ -87,6 +87,13 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   
   // check compatibility of dot and adjoint
   VERIFY_IS_APPROX(v1.dot(square * v2),              (square.adjoint() * v1).dot(v2));
+  
+  // like in testBasicStuff, test operator() to check const-qualification
+  int r = random<int>(0, rows-1),
+      c = random<int>(0, cols-1);
+  VERIFY_IS_APPROX(m1.conjugate()(r,c), conj(m1(r,c)));
+  VERIFY_IS_APPROX(m1.adjoint()(c,r), conj(m1(r,c)));
+  
 }
 
 void EigenTest::testAdjoint()
diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp
index f4bea3e46..dceeee564 100644
--- a/test/basicstuff.cpp
+++ b/test/basicstuff.cpp
@@ -61,6 +61,9 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   Scalar s1 = random<Scalar>(),
          s2 = random<Scalar>();
   
+  int r = random<int>(0, rows-1),
+      c = random<int>(0, cols-1);
+  
   // test Fuzzy.h and Zero.h.
   VERIFY_IS_APPROX(               v1,    v1);
   VERIFY_IS_NOT_APPROX(           v1,    2*v1);
@@ -75,6 +78,14 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   VERIFY_IS_NOT_MUCH_SMALLER_THAN(m1,    m1);
   VERIFY_IS_APPROX(               mzero, m1-m1);
   
+  // always test operator() on each read-only expression class,
+  // in order to check const-qualifiers.
+  // indeed, if an expression class (here Zero) is meant to be read-only,
+  // hence has no _write() method, the corresponding MatrixBase method (here zero())
+  // should return a const-qualified object so that it is the const-qualified
+  // operator() that gets called, which in turn calls _read().
+  VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::zero()(r,c), static_cast<Scalar>(1));
+  
   // test the linear structure, i.e. the following files:
   // Sum.h Difference.h Opposite.h ScalarMultiple.h
   VERIFY_IS_APPROX(-(-m1),                  m1);
@@ -100,6 +111,15 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
     VERIFY_IS_APPROX(m3,                    m2/s1);
   }
   
+  // again, test operator() to check const-qualification
+  VERIFY_IS_APPROX((-m1)(r,c), -(m1(r,c)));
+  VERIFY_IS_APPROX((m1-m2)(r,c), (m1(r,c))-(m2(r,c)));
+  VERIFY_IS_APPROX((m1+m2)(r,c), (m1(r,c))+(m2(r,c)));
+  VERIFY_IS_APPROX((s1*m1)(r,c), s1*(m1(r,c)));
+  VERIFY_IS_APPROX((m1*s1)(r,c), (m1(r,c))*s1);
+  if(NumTraits<Scalar>::HasFloatingPoint)
+    VERIFY_IS_APPROX((m1/s1)(r,c), (m1(r,c))/s1);
+  
   // begin testing Product.h: only associativity for now
   // (we use Transpose.h but this doesn't count as a test for it)
   VERIFY_IS_APPROX((m1*m1.transpose())*m2,  m1*(m1.transpose()*m2));
@@ -118,10 +138,14 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   
   // continue testing Product.h: lazyProduct
   VERIFY_IS_APPROX(square.lazyProduct(m1),  square*m1);
+  // again, test operator() to check const-qualification
+  s1 += square.lazyProduct(m1)(r,c);
   
-  // test Product.h together with Identity.h. This does test Identity.h.
+  // test Product.h together with Identity.h
   VERIFY_IS_APPROX(m1,                      identity*m1);
   VERIFY_IS_APPROX(v1,                      identity*v1);
+  // again, test operator() to check const-qualification
+  VERIFY_IS_APPROX(MatrixType::identity()(r,c), static_cast<Scalar>(r==c));
   
   // test FromArray.h
   Scalar* array1 = new Scalar[rows];