Index: src/Matrix.cpp
===================================================================
--- src/Matrix.cpp	(revision aafd7737af09b4111bf30cab050e67fa581a43c9)
+++ src/Matrix.cpp	(revision fee0794e90e50281a6f8a20468a8904a274cfd92)
@@ -13,4 +13,5 @@
 #include "vector.hpp"
 #include "VectorContent.hpp"
+#include "MatrixContent.hpp"
 
 #include <gsl/gsl_blas.h>
@@ -22,9 +23,11 @@
 Matrix::Matrix()
 {
-  content = gsl_matrix_calloc(NDIM, NDIM);
+  content = new MatrixContent();
+  content->content = gsl_matrix_calloc(NDIM, NDIM);
 }
 
 Matrix::Matrix(const double* src){
-  content = gsl_matrix_alloc(NDIM, NDIM);
+  content = new MatrixContent();
+  content->content = gsl_matrix_alloc(NDIM, NDIM);
   set(0,0, src[0]);
   set(1,0, src[1]);
@@ -41,9 +44,10 @@
 
 Matrix::Matrix(const Matrix &src){
-  content = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_matrix_memcpy(content,src.content);
-}
-
-Matrix::Matrix(gsl_matrix* _content) :
+  content = new MatrixContent();
+  content->content = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_matrix_memcpy(content->content,src.content->content);
+}
+
+Matrix::Matrix(MatrixContent* _content) :
   content(_content)
 {}
@@ -51,10 +55,11 @@
 Matrix::~Matrix()
 {
-  gsl_matrix_free(content);
+  gsl_matrix_free(content->content);
+  delete content;
 }
 
 void Matrix::one(){
-  gsl_matrix_free(content);
-  content = gsl_matrix_calloc(NDIM, NDIM);
+  gsl_matrix_free(content->content);
+  content->content = gsl_matrix_calloc(NDIM, NDIM);
   for(int i = NDIM;i--;){
     set(i,i,1.);
@@ -64,5 +69,5 @@
 Matrix &Matrix::operator=(const Matrix &src){
   if(&src!=this){
-    gsl_matrix_memcpy(content,src.content);
+    gsl_matrix_memcpy(content->content,src.content->content);
   }
   return *this;
@@ -70,10 +75,10 @@
 
 Matrix &Matrix::operator+=(const Matrix &rhs){
-  gsl_matrix_add(content, rhs.content);
+  gsl_matrix_add(content->content, rhs.content->content);
   return *this;
 }
 
 Matrix &Matrix::operator-=(const Matrix &rhs){
-  gsl_matrix_sub(content, rhs.content);
+  gsl_matrix_sub(content->content, rhs.content->content);
   return *this;
 }
@@ -98,6 +103,8 @@
 Matrix Matrix::operator*(const Matrix &rhs) const{
   gsl_matrix *res = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, this->content, rhs.content, 0.0, res);
-  return Matrix(res);
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content->content, rhs.content->content, 0.0, res);
+  MatrixContent *content= new MatrixContent();
+  content->content = res;
+  return Matrix(content);
 }
 
@@ -105,5 +112,5 @@
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  return *gsl_matrix_ptr (content, i, j);
+  return *gsl_matrix_ptr (content->content, i, j);
 }
 
@@ -111,5 +118,5 @@
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  return gsl_matrix_get(content, i, j);
+  return gsl_matrix_get(content->content, i, j);
 }
 
@@ -117,5 +124,5 @@
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  gsl_matrix_set(content,i,j,value);
+  gsl_matrix_set(content->content,i,j,value);
 }
 
@@ -150,5 +157,5 @@
 
 Matrix &Matrix::operator*=(const double factor){
-  gsl_matrix_scale(content, factor);
+  gsl_matrix_scale(content->content, factor);
   return *this;
 }
@@ -197,5 +204,5 @@
 Vector operator*(const Matrix &mat,const Vector &vec){
   gsl_vector *res = gsl_vector_calloc(NDIM);
-  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res);
+  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content->content, vec.content->content, 0.0, res);
   VectorContent *content = new VectorContent();
   content->content = res;
Index: src/Matrix.hpp
===================================================================
--- src/Matrix.hpp	(revision aafd7737af09b4111bf30cab050e67fa581a43c9)
+++ src/Matrix.hpp	(revision fee0794e90e50281a6f8a20468a8904a274cfd92)
@@ -9,5 +9,4 @@
 #define MATRIX_HPP_
 
-#include <gsl/gsl_matrix.h>
 #include <iosfwd>
 #include "defs.hpp"
@@ -18,4 +17,5 @@
 
 class Vector;
+class MatrixContent;
 
 class Matrix
@@ -92,6 +92,6 @@
 
 private:
-  Matrix(gsl_matrix*);
-  gsl_matrix *content;
+  Matrix(MatrixContent*);
+  MatrixContent *content;
 };
 
Index: src/MatrixContent.hpp
===================================================================
--- src/MatrixContent.hpp	(revision fee0794e90e50281a6f8a20468a8904a274cfd92)
+++ src/MatrixContent.hpp	(revision fee0794e90e50281a6f8a20468a8904a274cfd92)
@@ -0,0 +1,25 @@
+/*
+ * MatrixContent.hpp
+ *
+ *  Created on: Jul 2, 2010
+ *      Author: crueger
+ */
+
+#ifndef MATRIXCONTENT_HPP_
+#define MATRIXCONTENT_HPP_
+
+/**
+ * !file
+ * The way GSL works does not allow for forward definitions of the structures.
+ * Because of this the pointer to the gsl_matrix struct is wrapped inside another
+ * (dumb) object that allows for forward definitions. All methods of this class
+ * are empty, to allow for a maximum of flexibility.
+ */
+
+#include <gsl/gsl_matrix.h>
+
+struct MatrixContent{
+  gsl_matrix * content;
+};
+
+#endif /* MATRIXCONTENT_HPP_ */
