diff --git a/src/GMatrix.cpp b/src/GMatrix.cpp
index 69a552a..7430af9 100644
--- a/src/GMatrix.cpp
+++ b/src/GMatrix.cpp
@@ -134,6 +134,51 @@ ostream& operator<< (ostream& os, const GMatrix& m)
 	if (row != m.m_rows-1)
 	  os << endl;
   }
+  return os;
+}
+
+
+/***************************************************************************
+ *                       GMatrix * Gvector multiplication                  *
+ ***************************************************************************/
+GVector GMatrix::operator* (const GVector& v) const
+{
+  // Raise an exception if the matrix and vector dimensions are not compatible
+  if (m_cols != v.m_num)
+    throw vec_mat_mismatch("GMatrix*GVector operator", v.m_num, m_rows, m_cols);
+
+  // Perform vector multiplication
+  GVector result(m_rows);
+  for (unsigned row = 0; row < m_rows; ++row) {
+    double sum = 0.0;
+	for (unsigned col = 0; col < m_cols; ++col)
+	  sum += (*this)(row,col) * v(col);
+	result(row) = sum;
+  }
+  return result;
+}
+
+
+/***************************************************************************
+ *                       GMatrix * GMatrix multiplication                  *
+ ***************************************************************************/
+GMatrix GMatrix::operator* (const GMatrix& m) const
+{
+  // Raise an exception if the matrix dimensions are not compatible
+  if (m_cols != m.m_rows)
+    throw GMatrix::dim_mult_mismatch("GMatrix multiplication", m_cols, m.m_rows);
+
+  // Perform matrix multiplication
+  GMatrix result(m_rows,m.m_cols);
+  for (unsigned row = 0; row < m_rows; ++row) {
+    for (unsigned col = 0; col < m.m_cols; ++col) {
+	  double sum = 0.0;
+	  for (unsigned i = 0; i < m_cols; ++i)
+	    sum += (*this)(row,i)*m(i,col);
+      result(row,col) = sum;
+    }
+  }
+  return result;
 }
 
 
diff --git a/src/GMatrix.hpp b/src/GMatrix.hpp
index 29c9f02..078957f 100644
--- a/src/GMatrix.hpp
+++ b/src/GMatrix.hpp
@@ -71,7 +71,7 @@ public:
   GMatrix operator+ (const GMatrix& m) const;
   GMatrix operator- (const GMatrix& m) const;
   GMatrix operator* (const GMatrix& m) const;
-  GVector operator* (const GVector& v) const;
+  virtual GVector operator* (const GVector& v) const;
   GMatrix operator+ (const double& d) const;
   GMatrix operator- (const double& d) const;
   GMatrix operator* (const double& d) const;
@@ -128,34 +128,14 @@ protected:
 
 
 /***************************************************************************
- *                   Generic inline members (to be herited)                *
+ *                  Generic inline members (to be inherited)               *
  ***************************************************************************/
-// Matrix/Vector multiplication: GMatrix * GVector
-inline
-GVector GMatrix::operator* (const GVector& v) const
-{
-  #if G_RANGE_CHECK
-  if (m_cols != v.m_num)
-    throw vec_mat_mismatch("GMatrix*GVector operator", v.m_num, m_rows, m_cols);
-  #endif
-  GVector result(m_rows);
-  for (unsigned row = 0; row < m_rows; ++row) {
-    double sum = 0.0;
-	for (unsigned col = 0; col < m_cols; ++col)
-	  sum += (*this)(row,col) * v(col);
-	result(row) = sum;
-  }
-  return result;
-}
-
 // Matrix unary addition operator: GMatrix += GMatrix
 inline
 GMatrix& GMatrix::operator+= (const GMatrix& m)
 {
-  #if G_RANGE_CHECK
   if (m_rows != m.m_rows || m_cols != m.m_cols)
     throw dim_add_mismatch("GMatrix += operator", m_rows, m.m_rows, m_cols, m.m_cols);
-  #endif
   for (unsigned i = 0; i < m_elements; ++i)
     m_data[i] += m.m_data[i];
   return *this;
@@ -165,10 +145,8 @@ GMatrix& GMatrix::operator+= (const GMatrix& m)
 inline
 GMatrix& GMatrix::operator-= (const GMatrix& m)
 {
-  #if G_RANGE_CHECK
   if (m_rows != m.m_rows || m_cols != m.m_cols)
     throw dim_add_mismatch("GMatrix -= operator", m_rows, m.m_rows, m_cols, m.m_cols);
-  #endif
   for (unsigned i = 0; i < m_elements; ++i)
     m_data[i] -= m.m_data[i];
   return *this;
@@ -313,7 +291,7 @@ GMatrix& GMatrix::operator= (const double& d)
 inline
 double& GMatrix::operator() (unsigned row, unsigned col)
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (row >= m_rows || col >= m_cols)
     throw out_of_range("GMatrix access", row, col, m_rows, m_cols);
   #endif
@@ -324,7 +302,7 @@ double& GMatrix::operator() (unsigned row, unsigned col)
 inline
 const double& GMatrix::operator() (unsigned row, unsigned col) const
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (row >= m_rows || col >= m_cols)
     throw out_of_range("GMatrix access", row, col, m_rows, m_cols);
   #endif
@@ -349,26 +327,6 @@ GMatrix GMatrix::operator- (const GMatrix& m) const
   return result;
 }
 
-// Matrix multiplication: GMatrix * GMatrix
-inline
-GMatrix GMatrix::operator* (const GMatrix& m) const
-{
-  #if G_RANGE_CHECK
-  if (m_cols != m.m_rows)
-    throw GMatrix::dim_mult_mismatch("GMatrix multiplication", m_cols, m.m_rows);
-  #endif
-  GMatrix result(m_rows,m.m_cols);
-  for (unsigned row = 0; row < m_rows; ++row) {
-    for (unsigned col = 0; col < m.m_cols; ++col) {
-	  double sum = 0.0;
-	  for (unsigned i = 0; i < m_cols; ++i)
-	    sum += (*this)(row,i)*m(i,col);
-      result(row,col) = sum;
-    }
-  }
-  return result;
-}
-
 // Matrix scaling: GMatrix + double
 inline 
 GMatrix GMatrix::operator+ (const double& d) const
diff --git a/src/GSymMatrix.cpp b/src/GSymMatrix.cpp
index 810f140..63b7d90 100644
--- a/src/GSymMatrix.cpp
+++ b/src/GSymMatrix.cpp
@@ -114,6 +114,29 @@ GSymMatrix& GSymMatrix::operator= (const GSymMatrix& m)
 
 
 /***************************************************************************
+ *                   GSymMatrix * GSymMatrix multiplication                *
+ ***************************************************************************/
+GSymMatrix GSymMatrix::operator* (const GSymMatrix& m) const
+{
+  // Raise an exception if the matrix dimensions are not compatible
+  if (m_cols != m.m_rows)
+    throw GMatrix::dim_mult_mismatch("GSymMatrix multiplication", m_cols, m.m_rows);
+
+  // Perform matrix multiplication
+  GSymMatrix result(m_rows,m.m_cols);
+  for (unsigned row = 0; row < m_rows; ++row) {
+    for (unsigned col = 0; col < m.m_cols; ++col) {
+	  double sum = 0.0;
+	  for (unsigned i = 0; i < m_cols; ++i)
+	    sum += (*this)(row,i)*m(i,col);
+      result(row,col) = sum;
+    }
+  }
+  return result;
+}
+
+
+/***************************************************************************
  *                     GSymMatrix Cholesky decomposition                   *
  * ----------------------------------------------------------------------- *
  * Inplace Cholesky decomposition inspired by Numerical Recipes algorithm. *
diff --git a/src/GSymMatrix.hpp b/src/GSymMatrix.hpp
index 6f84760..a47375d 100644
--- a/src/GSymMatrix.hpp
+++ b/src/GSymMatrix.hpp
@@ -111,14 +111,15 @@ private:
 inline
 GSymMatrix& GSymMatrix::operator= (const double& d) 
 { 
-  this->GMatrix::operator=(d); 
+  this->GMatrix::operator=(d);
+  return *this;
 }
 
 // Matrix element access operator: GSymMatrix(row,col)
 inline
 double& GSymMatrix::operator() (unsigned row, unsigned col)
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (row >= m_rows || col >= m_cols)
     throw out_of_range("GSymMatrix access", row, col, m_rows, m_cols);
   #endif
@@ -130,7 +131,7 @@ double& GSymMatrix::operator() (unsigned row, unsigned col)
 inline
 const double& GSymMatrix::operator() (unsigned row, unsigned col) const
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (row >= m_rows || col >= m_cols)
     throw out_of_range("GSymMatrix access", row, col, m_rows, m_cols);
   #endif
@@ -156,26 +157,6 @@ GSymMatrix GSymMatrix::operator- (const GSymMatrix& m) const
   return result;
 }
 
-// Matrix multiplication: GSymMatrix * GSymMatrix
-inline
-GSymMatrix GSymMatrix::operator* (const GSymMatrix& m) const
-{
-  #if G_RANGE_CHECK
-  if (m_cols != m.m_rows)
-    throw GMatrix::dim_mult_mismatch("GSymMatrix multiplication", m_cols, m.m_rows);
-  #endif
-  GSymMatrix result(m_rows,m.m_cols);
-  for (unsigned row = 0; row < m_rows; ++row) {
-    for (unsigned col = 0; col < m.m_cols; ++col) {
-	  double sum = 0.0;
-	  for (unsigned i = 0; i < m_cols; ++i)
-	    sum += (*this)(row,i)*m(i,col);
-      result(row,col) = sum;
-    }
-  }
-  return result;
-}
-
 // Matrix/Vector multiplication: GSymMatrix * GVector
 inline
 GVector GSymMatrix::operator* (const GVector& v) const
diff --git a/src/GVector.cpp b/src/GVector.cpp
index c250ef9..d073b1f 100644
--- a/src/GVector.cpp
+++ b/src/GVector.cpp
@@ -59,6 +59,7 @@ ostream& operator<< (ostream& os, const GVector& v)
 	if (i != v.m_num-1)
 	  os << ", ";
   }
+  return os;
 }
 
 
@@ -67,10 +68,8 @@ ostream& operator<< (ostream& os, const GVector& v)
  ***************************************************************************/
 GVector cross(const GVector &a, const GVector &b)
 {
-  #if G_RANGE_CHECK
   if (a.m_num != b.m_num)
     throw GVector::dim_mismatch("GVector cross product", a.m_num, b.m_num);
-  #endif
   if (a.m_num != 3)
     throw GVector::bad_cross_dim(a.m_num);
   GVector result(3);
diff --git a/src/GVector.hpp b/src/GVector.hpp
index 7cb82be..c770b21 100644
--- a/src/GVector.hpp
+++ b/src/GVector.hpp
@@ -146,7 +146,7 @@ GVector::GVector(const GVector& v)
 inline
 double& GVector::operator() (unsigned inx)
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (inx >= m_num)
     throw out_of_range("GVector access", inx, m_num);
   #endif
@@ -157,7 +157,7 @@ double& GVector::operator() (unsigned inx)
 inline
 const double& GVector::operator() (unsigned inx) const
 {
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   if (inx >= m_num)
     throw out_of_range("const GVector access", inx, m_num);
   #endif
@@ -181,10 +181,8 @@ GVector& GVector::operator= (const GVector& v)
 inline
 GVector& GVector::operator+= (const GVector& v)
 {
-  #if G_RANGE_CHECK
   if (m_num != v.m_num)
     throw dim_mismatch("GVector += operator", m_num, v.m_num);
-  #endif
   for (unsigned i = 0; i < m_num; ++i)
     m_data[i] += v.m_data[i];
   return *this;
@@ -194,10 +192,8 @@ GVector& GVector::operator+= (const GVector& v)
 inline
 GVector& GVector::operator-= (const GVector& v)
 {
-  #if G_RANGE_CHECK
   if (m_num != v.m_num)
     throw dim_mismatch("GVector -= operator", m_num, v.m_num);
-  #endif
   for (unsigned i = 0; i < m_num; ++i)
     m_data[i] -= v.m_data[i];
   return *this;
@@ -327,10 +323,8 @@ GVector operator- (const double& a, const GVector& b)
 inline
 double operator* (const GVector& a, const GVector& b)
 {
-  #if G_RANGE_CHECK
   if (a.m_num != b.m_num)
     throw GVector::dim_mismatch("GVector scalar product", a.m_num, b.m_num);
-  #endif
   double result = 0.0;
   for (unsigned i = 0; i < a.m_num; ++i)
     result += (a.m_data[i] * b.m_data[i]);
diff --git a/test/test_GMatrix.cpp b/test/test_GMatrix.cpp
index 818ad91..24d2dc0 100644
--- a/test/test_GMatrix.cpp
+++ b/test/test_GMatrix.cpp
@@ -304,7 +304,7 @@ int main(void)
     cout << e.what() << endl;
 	throw;
   }
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   try {
     GMatrix test3(3,3);
     test3(3,3) = 1.0;
diff --git a/test/test_GSymMatrix.cpp b/test/test_GSymMatrix.cpp
index 2396943..9974d9a 100644
--- a/test/test_GSymMatrix.cpp
+++ b/test/test_GSymMatrix.cpp
@@ -356,7 +356,7 @@ int main(void)
     cout << e.what() << endl;
 	throw;
   }
-  #if G_RANGE_CHECK
+  #if defined(G_RANGE_CHECK)
   try {
     GSymMatrix test3(3,3);
     test3(3,3) = 1.0;
diff --git a/test/test_GVector.cpp b/test/test_GVector.cpp
index fe560be..c6ce1c4 100644
--- a/test/test_GVector.cpp
+++ b/test/test_GVector.cpp
@@ -18,6 +18,12 @@ using namespace std;
 
 int main(void)
 {
+  // Dump header
+  cout << endl;
+  cout << "*************************" << endl;
+  cout << "* GVector class testing *" << endl;
+  cout << "*************************" << endl;
+
   // Set parameters
   unsigned num = 5;
 
@@ -64,7 +70,7 @@ int main(void)
 	cout << "  Size = " << test3.size() << endl;
     //
 	cout << "  GVector out of range access: ";
-    #if G_RANGE_CHECK
+    #if defined(G_RANGE_CHECK)
 	cout << test3(3) << endl;
 	#endif
   }
@@ -323,4 +329,4 @@ int main(void)
     cout << e.what() << endl;
   }
 
-}
\ No newline at end of file
+}