diff --git a/src/GMatrix.cpp b/src/GMatrix.cpp
index a092062..af86b75 100644
--- a/src/GMatrix.cpp
+++ b/src/GMatrix.cpp
@@ -495,6 +495,20 @@ GMatrix& GMatrix::operator/= (const double& d)
  ==========================================================================*/
 
 /***************************************************************************
+ *                       Set all matrix elements to zero                   *
+ ***************************************************************************/
+void GMatrix::clear()
+{
+  // Loop over all matrix elements
+  for (int i = 0; i < m_elements; ++i)
+    m_data[i] = 0.0;
+  
+  // Return
+  return;
+}
+
+
+/***************************************************************************
  *                    Get minimum physical matrix element                  *
  ***************************************************************************/
 double GMatrix::min() const
diff --git a/src/GMatrix.hpp b/src/GMatrix.hpp
index e308811..e6c31c8 100644
--- a/src/GMatrix.hpp
+++ b/src/GMatrix.hpp
@@ -76,6 +76,7 @@ public:
   // Matrix functions
   virtual int     rows() const { return m_rows; }
   virtual int     cols() const { return m_cols; }
+  virtual void    clear();
   virtual double  min() const;
   virtual double  max() const;
   virtual double  sum() const;
diff --git a/src/GSparseMatrix.hpp b/src/GSparseMatrix.hpp
index 51ddf50..b5376d4 100644
--- a/src/GSparseMatrix.hpp
+++ b/src/GSparseMatrix.hpp
@@ -84,7 +84,8 @@ public:
   // Matrix functions
   // USE_BASE: virtual int rows() const { return m_rows; }
   // USE_BASE: virtual int cols() const { return m_cols; }
-  virtual double  fill() const { return (double(m_elements)/double(m_rows*m_cols)); }
+  virtual void    clear();
+  virtual double  fill() const;
   virtual double  min() const;
   virtual double  max() const;
   virtual double  sum() const;
@@ -244,6 +245,21 @@ GSparseMatrix operator/ (const GSparseMatrix& a, const double& b)
   return result;
 }
 
+// Matrix clear function
+inline
+void GSparseMatrix::clear()
+{
+  free_elements(0, m_elements);
+  return;
+}
+
+// Matrix fill function
+inline
+double GSparseMatrix::fill() const
+{
+  return (double(m_elements)/double(m_rows*m_cols));
+}
+
 // Matrix transpose function
 inline
 GSparseMatrix transpose(const GSparseMatrix& m)
diff --git a/src/GSymMatrix.hpp b/src/GSymMatrix.hpp
index 8fcd5ad..8448adb 100644
--- a/src/GSymMatrix.hpp
+++ b/src/GSymMatrix.hpp
@@ -77,8 +77,9 @@ public:
   virtual GSymMatrix& operator/= (const double& d);
 
   // Matrix functions
-  // USE_BASE: virtual int rows() const { return m_rows; }
-  // USE_BASE: virtual int cols() const { return m_cols; }
+  // USE_BASE: virtual int      rows() const { return m_rows; }
+  // USE_BASE: virtual int      cols() const { return m_cols; }
+  // USE_BASE: virtual void     clear();
   // USE_BASE: virtual double   min() const;
   // USE_BASE: virtual double   max() const;
   virtual double  sum() const;
@@ -92,12 +93,6 @@ public:
   virtual GVector cholesky_solver(const GVector& v, int compress = 1);
   virtual void    cholesky_invert(int compress = 1);
     
-  // Exception: All matrix elements are zero
-  class all_zero : public GException {
-  public:
-    all_zero(string origin);
-  };
-
 private:
   void set_inx(void);  // Set indices of non-zero rows/columns
   int  m_num_inx;      // Number of non-zero rows/columns
diff --git a/test/test_GSparseMatrix.cpp b/test/test_GSparseMatrix.cpp
index 91e4c8b..05b0879 100644
--- a/test/test_GSparseMatrix.cpp
+++ b/test/test_GSparseMatrix.cpp
@@ -349,6 +349,15 @@ void test_assign_values(const GSparseMatrix& m_test)
 	  throw;
 	}
 	cout << ".";
+	//
+	result.clear();
+    if (!check_matrix(result, 0.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {
+      cout << endl << "TEST ERROR: Corrupt GSparseMatrix.clear()." << endl;
+	  cout << "result:" << endl << result << endl;
+	  cout << "m_test:" << endl << m_test << endl;
+	  throw;
+	}
+	cout << ".";
   }
   catch (exception &e) {
     cout << e.what() << endl;