Index: src/LinearAlgebra/Makefile.am
===================================================================
--- src/LinearAlgebra/Makefile.am	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/Makefile.am	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -19,4 +19,5 @@
   linearsystemofequations.cpp \
   Matrix.cpp \
+  MatrixContent.cpp \
   Plane.cpp \
   Space.cpp \
@@ -34,8 +35,10 @@
   linearsystemofequations.hpp \
   Matrix.hpp \
+  MatrixContent.hpp \
   Plane.hpp \
   Space.hpp \
   vector_ops.hpp \
   Vector.hpp \
+  VectorContent.hpp \
   VectorInterface.hpp
 
Index: src/LinearAlgebra/Matrix.cpp
===================================================================
--- src/LinearAlgebra/Matrix.cpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/Matrix.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -26,6 +26,6 @@
 #include "Helpers/Assert.hpp"
 #include "LinearAlgebra/Vector.hpp"
-#include "VectorContent.hpp"
-#include "MatrixContent.hpp"
+#include "LinearAlgebra/VectorContent.hpp"
+#include "LinearAlgebra/MatrixContent.hpp"
 
 #include <gsl/gsl_blas.h>
@@ -41,36 +41,29 @@
 Matrix::Matrix()
 {
-  content = new MatrixContent();
-  content->content = gsl_matrix_calloc(NDIM, NDIM);
-  createViews();
-}
-
-Matrix::Matrix(const double* src){
-  content = new MatrixContent();
-  content->content = gsl_matrix_alloc(NDIM, NDIM);
-  set(0,0, src[0]);
-  set(1,0, src[1]);
-  set(2,0, src[2]);
-
-  set(0,1, src[3]);
-  set(1,1, src[4]);
-  set(2,1, src[5]);
-
-  set(0,2, src[6]);
-  set(1,2, src[7]);
-  set(2,2, src[8]);
-  createViews();
-}
-
-Matrix::Matrix(const Matrix &src){
-  content = new MatrixContent();
-  content->content = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_matrix_memcpy(content->content,src.content->content);
-  createViews();
-}
-
-Matrix::Matrix(MatrixContent* _content) :
-  content(_content)
-{
+  content = new MatrixContent(NDIM, NDIM);
+  createViews();
+}
+
+Matrix::Matrix(const double* src)
+{
+  content = new MatrixContent(NDIM, NDIM, src);
+  createViews();
+}
+
+Matrix::Matrix(const Matrix &src)
+{
+  content = new MatrixContent(src.content);
+  createViews();
+}
+
+Matrix::Matrix(const MatrixContent &src)
+{
+  content = new MatrixContent(src);
+  createViews();
+}
+
+Matrix::Matrix(MatrixContent* _content)
+{
+  content = new MatrixContent(_content);
   createViews();
 }
@@ -86,5 +79,4 @@
   }
   delete diagonal_ptr;
-  gsl_matrix_free(content->content);
   delete content;
 }
@@ -107,17 +99,9 @@
 
 void Matrix::setIdentity(){
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      set(i,j,i==j);
-    }
-  }
+  content->setIdentity();
 }
 
 void Matrix::setZero(){
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      set(i,j,0.);
-    }
-  }
+  content->setZero();
 }
 
@@ -135,27 +119,31 @@
 }
 
-Matrix &Matrix::operator=(const Matrix &src){
-  if(&src!=this){
-    gsl_matrix_memcpy(content->content,src.content->content);
-  }
-  return *this;
-}
-
-const Matrix &Matrix::operator+=(const Matrix &rhs){
-  gsl_matrix_add(content->content, rhs.content->content);
-  return *this;
-}
-
-const Matrix &Matrix::operator-=(const Matrix &rhs){
-  gsl_matrix_sub(content->content, rhs.content->content);
-  return *this;
-}
-
-const Matrix &Matrix::operator*=(const Matrix &rhs){
-  (*this) = (*this)*rhs;
-  return *this;
-}
-
-const Matrix Matrix::operator+(const Matrix &rhs) const{
+Matrix &Matrix::operator=(const Matrix &src)
+{
+  // MatrixContent checks for self-assignment
+  *content = *(src.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator+=(const Matrix &rhs)
+{
+  *content += *(rhs.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator-=(const Matrix &rhs)
+{
+  *content -= *(rhs.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator*=(const Matrix &rhs)
+{
+  (*content) *= (*rhs.content);
+  return *this;
+}
+
+const Matrix Matrix::operator+(const Matrix &rhs) const
+{
   Matrix tmp = *this;
   tmp+=rhs;
@@ -163,5 +151,6 @@
 }
 
-const Matrix Matrix::operator-(const Matrix &rhs) const{
+const Matrix Matrix::operator-(const Matrix &rhs) const
+{
   Matrix tmp = *this;
   tmp-=rhs;
@@ -169,56 +158,58 @@
 }
 
-const Matrix Matrix::operator*(const Matrix &rhs) const{
-  gsl_matrix *res = gsl_matrix_alloc(NDIM, NDIM);
-  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);
-}
-
-double &Matrix::at(size_t i, size_t j){
-  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->content, i, j);
-}
-
-const double Matrix::at(size_t i, size_t j) const{
-  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->content, i, j);
-}
-
-Vector &Matrix::row(size_t i){
+const Matrix Matrix::operator*(const Matrix &rhs) const
+{
+  Matrix tmp(content);
+  tmp *= rhs;
+  return tmp;
+}
+
+double &Matrix::at(size_t i, size_t j)
+{
+  return content->at(i,j);
+}
+
+const double Matrix::at(size_t i, size_t j) const
+{
+  return content->at(i,j);
+}
+
+Vector &Matrix::row(size_t i)
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *rows_ptr[i];
 }
 
-const Vector &Matrix::row(size_t i) const{
+const Vector &Matrix::row(size_t i) const
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *rows_ptr[i];
 }
 
-Vector &Matrix::column(size_t i){
+Vector &Matrix::column(size_t i)
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *columns_ptr[i];
 }
 
-const Vector &Matrix::column(size_t i) const{
+const Vector &Matrix::column(size_t i) const
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *columns_ptr[i];
 }
 
-Vector &Matrix::diagonal(){
+Vector &Matrix::diagonal()
+{
   return *diagonal_ptr;
 }
 
-const Vector &Matrix::diagonal() const{
+const Vector &Matrix::diagonal() const
+{
   return *diagonal_ptr;
 }
 
-void Matrix::set(size_t i, size_t j, const double value){
-  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->content,i,j,value);
+void Matrix::set(size_t i, size_t j, const double value)
+{
+  content->set(i,j, value);
 }
 
@@ -253,9 +244,8 @@
 }
 
-Matrix Matrix::transpose() const{
-  MatrixContent *newContent = new MatrixContent();
-  newContent->content = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_matrix_transpose_memcpy(newContent->content, content->content);
-  Matrix res = Matrix(newContent);
+Matrix Matrix::transpose() const
+{
+  std::cout << "const Matrix::transpose()." << std::endl;
+  Matrix res = Matrix(const_cast<const MatrixContent *>(content)->transpose());
   return res;
 }
@@ -263,24 +253,12 @@
 Matrix &Matrix::transpose()
 {
-  double tmp;
-  for (int i=0;i<NDIM;i++)
-    for (int j=i+1;j<NDIM;j++) {
-      tmp = at(j,i);
-      at(j,i) = at(i,j);
-      at(i,j) = tmp;
-    }
-  return *this;
-}
-
+  std::cout << "Matrix::transpose()." << std::endl;
+  content->transpose();
+  return *this;
+}
 
 Vector Matrix::transformToEigenbasis()
 {
-  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
-  gsl_vector *eval = gsl_vector_alloc(NDIM);
-  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_eigen_symmv(content->content, eval, evec, T);
-  gsl_eigen_symmv_free(T);
-  gsl_matrix_memcpy(content->content, evec);
-  gsl_matrix_free(evec);
+  gsl_vector *eval = content->transformToEigenbasis();
   Vector evalues(gsl_vector_get(eval,0), gsl_vector_get(eval,1), gsl_vector_get(eval,2));
   gsl_vector_free(eval);
@@ -288,28 +266,24 @@
 }
 
-const Matrix &Matrix::operator*=(const double factor){
-  gsl_matrix_scale(content->content, factor);
-  return *this;
-}
-
-const Matrix operator*(const double factor,const Matrix& mat){
+const Matrix &Matrix::operator*=(const double factor)
+    {
+  *content *= factor;
+  return *this;
+}
+
+const Matrix operator*(const double factor,const Matrix& mat)
+{
   Matrix tmp = mat;
-  tmp*=factor;
-  return tmp;
-}
-
-const Matrix operator*(const Matrix &mat,const double factor){
+  return tmp *= factor;
+}
+
+const Matrix operator*(const Matrix &mat,const double factor)
+{
   return factor*mat;
 }
 
-bool Matrix::operator==(const Matrix &rhs) const{
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
-        return false;
-      }
-    }
-  }
-  return true;
+bool Matrix::operator==(const Matrix &rhs) const
+{
+  return (*content == *(rhs.content));
 }
 
@@ -333,5 +307,6 @@
 };
 
-ostream &operator<<(ostream &ost,const Matrix &mat){
+ostream &operator<<(ostream &ost,const Matrix &mat)
+{
   for(int i = 0;i<NDIM;++i){
     ost << "\n";
@@ -345,11 +320,11 @@
 }
 
-Vector operator*(const Matrix &mat,const Vector &vec){
-  Vector res;
-  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content->content, vec.content->content, 0.0, res.content->content);
-  return res;
-}
-
-Vector &operator*=(Vector& lhs,const Matrix &mat){
+Vector operator*(const Matrix &mat,const Vector &vec)
+{
+  return (*mat.content) * vec;
+}
+
+Vector &operator*=(Vector& lhs,const Matrix &mat)
+{
   lhs = mat*lhs;
   return lhs;
Index: src/LinearAlgebra/Matrix.hpp
===================================================================
--- src/LinearAlgebra/Matrix.hpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/Matrix.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -43,4 +43,5 @@
   Matrix(const double*);
   Matrix(const Matrix&);
+  Matrix(const MatrixContent&);
   virtual ~Matrix();
 
Index: src/LinearAlgebra/MatrixContent.cpp
===================================================================
--- src/LinearAlgebra/MatrixContent.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ src/LinearAlgebra/MatrixContent.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -0,0 +1,344 @@
+/*
+ * MatrixContent.cpp
+ *
+ *  Created on: Nov 14, 2010
+ *      Author: heber
+ */
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Helpers/MemDebug.hpp"
+
+#include "LinearAlgebra/Matrix.hpp"
+#include "Helpers/Assert.hpp"
+#include "Exceptions/NotInvertibleException.hpp"
+#include "Helpers/fast_functions.hpp"
+#include "Helpers/Assert.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/VectorContent.hpp"
+#include "LinearAlgebra/MatrixContent.hpp"
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_vector.h>
+#include <cmath>
+#include <iostream>
+
+using namespace std;
+
+
+/** Constructor for class MatrixContent.
+ * \param rows number of rows
+ * \param columns number of columns
+ */
+MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
+    rows(_rows),
+    columns(_columns)
+{
+  content = gsl_matrix_calloc(rows, columns);
+}
+
+/** Constructor for class MatrixContent.
+ * \param rows number of rows
+ * \param columns number of columns
+ * \param *src array with components to initialize matrix with
+ */
+MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
+        rows(_rows),
+        columns(_columns)
+{
+  content = gsl_matrix_calloc(rows, columns);
+  set(0,0, src[0]);
+  set(1,0, src[1]);
+  set(2,0, src[2]);
+
+  set(0,1, src[3]);
+  set(1,1, src[4]);
+  set(2,1, src[5]);
+
+  set(0,2, src[6]);
+  set(1,2, src[7]);
+  set(2,2, src[8]);
+}
+
+/** Constructor for class MatrixContent.
+ * We embed the given gls_matrix pointer within this class and set it to NULL
+ * afterwards.
+ * \param *src source gsl_matrix vector to embed within this class
+ */
+MatrixContent::MatrixContent(gsl_matrix *&src) :
+        rows(src->size1),
+        columns(src->size2)
+{
+  content = gsl_matrix_alloc(src->size1, src->size2);
+  gsl_matrix_memcpy(content,src);
+//  content = src;
+//  src = NULL;
+}
+
+/** Copy constructor for class MatrixContent.
+ * \param &src reference to source MatrixContent
+ */
+MatrixContent::MatrixContent(const MatrixContent &src) :
+        rows(src.rows),
+        columns(src.columns)
+{
+  content = gsl_matrix_alloc(src.rows, src.columns);
+  gsl_matrix_memcpy(content,src.content);
+}
+
+/** Copy constructor for class MatrixContent.
+ * \param *src pointer to source MatrixContent
+ */
+MatrixContent::MatrixContent(const MatrixContent *src) :
+        rows(src->rows),
+        columns(src->columns)
+{
+  ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
+  content = gsl_matrix_alloc(src->rows, src->columns);
+  gsl_matrix_memcpy(content,src->content);
+}
+
+/** Destructor for class MatrixContent.
+ */
+MatrixContent::~MatrixContent()
+{
+  gsl_matrix_free(content);
+}
+
+/** Set matrix to identity.
+ */
+void MatrixContent::setIdentity()
+{
+  for(int i=rows;i--;){
+    for(int j=columns;j--;){
+      set(i,j,i==j);
+    }
+  }
+}
+
+/** Set all matrix components to zero.
+ */
+void MatrixContent::setZero()
+{
+  for(int i=rows;i--;){
+    for(int j=columns;j--;){
+      set(i,j,0.);
+    }
+  }
+}
+
+/** Copy operator for MatrixContent with self-assignment check.
+ * \param &src matrix to compare to
+ * \return reference to this
+ */
+MatrixContent &MatrixContent::operator=(const MatrixContent &src)
+{
+  if(&src!=this){
+    gsl_matrix_memcpy(content,src.content);
+  }
+  return *this;
+}
+
+/** Addition operator.
+ * \param &rhs matrix to add
+ * \return reference to this
+ */
+const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
+{
+  gsl_matrix_add(content, rhs.content);
+  return *this;
+}
+
+/** Subtraction operator.
+ * \param &rhs matrix to subtract
+ * \return reference to this
+ */
+const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
+    {
+  gsl_matrix_sub(content, rhs.content);
+  return *this;
+}
+
+/** Multiplication operator.
+ * Note that here matrix have to have same dimensions.
+ * \param &rhs matrix to multiply with
+ * \return reference to this
+ */
+const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
+{
+  ASSERT(rows == rhs.rows,
+      "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
+  ASSERT(columns == rhs.columns,
+      "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
+  (*this) = (*this)*rhs;
+  return *this;
+}
+
+/** Multiplication with copy operator.
+ * \param &rhs matrix to multiply with
+ * \return reference to newly allocated MatrixContent
+ */
+const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
+{
+  gsl_matrix *res = gsl_matrix_alloc(rows, columns);
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
+  // gsl_matrix is taken over by constructor, hence no free
+  MatrixContent tmp(res);
+  gsl_matrix_free(res);
+  return tmp;
+}
+
+/** Accessor for manipulating component (i,j).
+ * \param i row number
+ * \param j column number
+ * \return reference to component (i,j)
+ */
+double &MatrixContent::at(size_t i, size_t j)
+{
+  ASSERT((i>=0) && (i<rows),
+      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
+  ASSERT((j>=0) && (j<columns),
+      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  return *gsl_matrix_ptr (content, i, j);
+}
+
+/** Constant accessor for (value of) component (i,j).
+ * \param i row number
+ * \param j column number
+ * \return const component (i,j)
+ */
+const double MatrixContent::at(size_t i, size_t j) const
+{
+  ASSERT((i>=0) && (i<rows),
+      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
+  ASSERT((j>=0) && (j<columns),
+      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  return gsl_matrix_get(content, i, j);
+}
+
+/** Setter for component (i,j).
+ * \param i row numbr
+ * \param j column numnber
+ * \param value value to set componnt (i,j) to
+ */
+void MatrixContent::set(size_t i, size_t j, const double value)
+{
+  ASSERT((i>=0) && (i<rows),
+      "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
+  ASSERT((j>=0) && (j<columns),
+      "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  gsl_matrix_set(content,i,j,value);
+}
+
+/** Return transposed matrix.
+ * \return new matrix that is transposed of this.
+ */
+MatrixContent MatrixContent::transpose() const
+{
+  std::cout << "const MatrixContent::transpose()." << std::endl;
+  gsl_matrix *res = gsl_matrix_alloc(rows, columns);
+  gsl_matrix_transpose_memcpy(res, content);
+  MatrixContent newContent(res);
+  gsl_matrix_free(res);
+  return newContent;
+}
+
+/** Turn this matrix into its transposed.
+ * Note that this is only possible if rows == columns.
+ */
+MatrixContent &MatrixContent::transpose()
+{
+  std::cout << "MatrixContent::transpose()." << std::endl;
+  ASSERT( rows == columns,
+      "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
+  double tmp;
+  for (size_t i=0;i<rows;i++)
+    for (size_t j=i+1;j<rows;j++) {
+      tmp = at(j,i);
+      at(j,i) = at(i,j);
+      at(i,j) = tmp;
+    }
+  return *this;
+}
+
+/** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
+ * \warn return vector has to be freed'd
+ * \return gsl_vector pointer to vector of eigenvalues
+ */
+gsl_vector* MatrixContent::transformToEigenbasis()
+{
+  ASSERT (rows == columns,
+      "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
+  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
+  gsl_vector *eval = gsl_vector_alloc(rows);
+  gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
+  gsl_eigen_symmv(content, eval, evec, T);
+  gsl_eigen_symmv_free(T);
+  gsl_matrix_memcpy(content, evec);
+  gsl_matrix_free(evec);
+  return eval;
+}
+
+/** Scalar multiplication operator.
+ * \param factor factor to scale with
+ */
+const MatrixContent &MatrixContent::operator*=(const double factor)
+{
+  gsl_matrix_scale(content, factor);
+  return *this;
+}
+
+/** Scalar multiplication and copy operator.
+ * \param factor factor to scale with
+ * \param &mat MatrixContent to scale
+ * \return copied and scaled MatrixContent
+ */
+const MatrixContent operator*(const double factor,const MatrixContent& mat)
+{
+  MatrixContent tmp = mat;
+  tmp*=factor;
+  return tmp;
+}
+
+/** Scalar multiplication and copy operator (with operands exchanged).
+ * \param &mat MatrixContent to scale
+ * \param factor factor to scale with
+ * \return copied and scaled MatrixContent
+ */
+const MatrixContent operator*(const MatrixContent &mat,const double factor)
+{
+  return factor*mat;
+}
+
+/** Equality operator.
+ * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
+ * \param &rhs MatrixContent to checks against
+ */
+bool MatrixContent::operator==(const MatrixContent &rhs) const
+    {
+  if ((rows == rhs.rows) && (columns == rhs.columns)) {
+    for(int i=rows;i--;){
+      for(int j=columns;j--;){
+        if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
+          return false;
+        }
+      }
+    }
+    return true;
+  }
+  return false;
+}
+
+Vector operator*(const MatrixContent &mat,const Vector &vec)
+{
+  Vector res;
+  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
+  return res;
+}
Index: src/LinearAlgebra/MatrixContent.hpp
===================================================================
--- src/LinearAlgebra/MatrixContent.hpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/MatrixContent.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -9,17 +9,77 @@
 #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.
+/** MatrixContent is a wrapper for gsl_matrix.
+ *
+ * The GNU Scientific Library contaisn some very well written routines for
+ * linear algebra problems. However, it's syntax is C and hence it does not
+ * lend itself to readable written code, i.e. C = A * B, where A, B, and C
+ * are all matrices. Writing code this way is very convenient, concise and
+ * also in the same manner as one would type in MatLab.
+ * However, with C++ and its feature to overload functions and its operator
+ * functions such syntax is easy to create.
+ *
+ * Hence, this class is a C++ wrapper to gsl_matrix. There already some other
+ * libraries, however they fall short for the following reasons:
+ * -# gslwrap: not very well written and uses floats instead of doubles
+ * -# gslcpp: last change is from 2007 and only very few commits
+ * -# o2scl: basically, the same functionality as gsl only written in C++,
+ *    however it uses GPLv3 license which is inconvenient for MoleCuilder.
+ *
+ * <h1>Howto use MatrixContent</h1>
+ *
+ * One should not use MatrixContent directly but either have it as a member
+ * variable in a specialized class or inherit from it.
+ *
  */
 
 #include <gsl/gsl_matrix.h>
 
-struct MatrixContent{
+class Vector;
+
+class MatrixContent
+{
+  friend Vector operator*(const MatrixContent &mat,const Vector &vec);
+  friend class Matrix;
+public:
+  MatrixContent(size_t rows, size_t columns);
+  MatrixContent(size_t rows, size_t columns, const double *src);
+  MatrixContent(gsl_matrix *&src);
+  MatrixContent(const MatrixContent &src);
+  MatrixContent(const MatrixContent *src);
+  ~MatrixContent();
+
+  // Accessing
+  double &at(size_t i, size_t j);
+  const double at(size_t i, size_t j) const;
+  void set(size_t i, size_t j, const double value);
+
+  // Transformations
+  MatrixContent transpose() const;
+  MatrixContent &transpose();
+  gsl_vector* transformToEigenbasis();
+
+  // Initializing
+  void setIdentity();
+  void setZero();
+
+  // Operators
+  MatrixContent &operator=(const MatrixContent &src);
+  const MatrixContent &operator+=(const MatrixContent &rhs);
+  const MatrixContent &operator-=(const MatrixContent &rhs);
+  const MatrixContent &operator*=(const MatrixContent &rhs);
+  const MatrixContent operator*(const MatrixContent &rhs) const;
+  const MatrixContent &operator*=(const double factor);
+  bool operator==(const MatrixContent &rhs) const;
+
+
+private:
   gsl_matrix * content;
+  const size_t rows;      // store for internal purposes
+  const size_t columns;  // store for internal purposes
 };
 
+const MatrixContent operator*(const double factor,const MatrixContent& mat);
+const MatrixContent operator*(const MatrixContent &mat,const double factor);
+Vector operator*(const MatrixContent &mat,const Vector &vec);
+
 #endif /* MATRIXCONTENT_HPP_ */
Index: src/LinearAlgebra/Vector.hpp
===================================================================
--- src/LinearAlgebra/Vector.hpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/Vector.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -23,4 +23,5 @@
 class Vector;
 class Matrix;
+class MatrixContent;
 struct VectorContent;
 
@@ -31,5 +32,5 @@
  */
 class Vector : public Space{
-  friend Vector operator*(const Matrix&,const Vector&);
+  friend Vector operator*(const MatrixContent&,const Vector&);
   friend class Matrix;
 public:
Index: src/LinearAlgebra/gslmatrix.cpp
===================================================================
--- src/LinearAlgebra/gslmatrix.cpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/gslmatrix.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -27,5 +27,10 @@
 
 #include <cassert>
+#include <iostream>
 #include <gsl/gsl_linalg.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
 
 /** Constructor of class GSLMatrix.
@@ -46,4 +51,19 @@
   matrix = gsl_matrix_alloc(rows, columns);
   gsl_matrix_memcpy (matrix, src->matrix);
+};
+
+/** Copy constructor of class GSLMatrix.
+ * We take over pointer to gsl_matrix and set the parameter pointer to NULL
+ * afterwards.
+ * \param m row count
+ * \param n column count
+ * \param *&src source gsl_matrix
+ */
+GSLMatrix::GSLMatrix(size_t m, size_t n, gsl_matrix *&src) :
+    rows(m),
+    columns(n)
+{
+  matrix = src;
+  src = NULL;
 };
 
@@ -90,5 +110,5 @@
  * \return (m,n)-th element of matrix
  */
-double GSLMatrix::Get(size_t m, size_t n)
+double GSLMatrix::Get(size_t m, size_t n) const
 {
   return gsl_matrix_get (matrix, m, n);
@@ -121,8 +141,24 @@
  * \return const pointer to \a m-th element
  */
-const double *GSLMatrix::const_Pointer(size_t m, size_t n)
+const double *GSLMatrix::const_Pointer(size_t m, size_t n) const
 {
   return gsl_matrix_const_ptr (matrix, m, n);
 };
+
+/** Returns the number of rows of the matrix.
+ * \return number of rows
+ */
+size_t GSLMatrix::getRowCount() const
+{
+  return rows;
+}
+
+/** Returns the number of columns of the matrix.
+ * \return number of columns
+ */
+size_t GSLMatrix::getColumnCount() const
+{
+  return columns;
+}
 
 /* ========================== Initializing =============================== */
@@ -204,5 +240,5 @@
  * \return true - is null, false - else
  */
-bool GSLMatrix::IsNull()
+bool GSLMatrix::IsNull() const
 {
   return gsl_matrix_isnull (matrix);
@@ -212,5 +248,5 @@
  * \return true - is positive, false - else
  */
-bool GSLMatrix::IsPositive()
+bool GSLMatrix::IsPositive() const
 {
   return gsl_matrix_ispos (matrix);
@@ -220,5 +256,5 @@
  * \return true - is negative, false - else
  */
-bool GSLMatrix::IsNegative()
+bool GSLMatrix::IsNegative() const
 {
   return gsl_matrix_isneg (matrix);
@@ -228,5 +264,5 @@
  * \return true - is non-negative, false - else
  */
-bool GSLMatrix::IsNonNegative()
+bool GSLMatrix::IsNonNegative() const
 {
   return gsl_matrix_isnonneg (matrix);
@@ -237,5 +273,5 @@
  * \return true - matrix is positive-definite, false - else
  */
-bool GSLMatrix::IsPositiveDefinite()
+bool GSLMatrix::IsPositiveDefinite() const
 {
   if (rows != columns)  // only possible for square matrices.
@@ -249,5 +285,5 @@
  * if matrix is square, uses LU decomposition.
  */
-double GSLMatrix::Determinant()
+double GSLMatrix::Determinant() const
 {
   int signum = 0;
@@ -259,2 +295,65 @@
 };
 
+
+/* ============================ Properties ============================== */
+
+const GSLMatrix &GSLMatrix::operator+=(const GSLMatrix &rhs)
+{
+  gsl_matrix_add(matrix, rhs.matrix);
+  return *this;
+}
+
+const GSLMatrix &GSLMatrix::operator-=(const GSLMatrix &rhs)
+{
+  gsl_matrix_sub(matrix, rhs.matrix);
+  return *this;
+}
+
+const GSLMatrix &GSLMatrix::operator*=(const GSLMatrix &rhs)
+{
+  (*this) = (*this)*rhs;
+  return *this;
+}
+
+const GSLMatrix GSLMatrix::operator+(const GSLMatrix &rhs) const
+{
+  GSLMatrix tmp = *this;
+  tmp+=rhs;
+  return tmp;
+}
+
+const GSLMatrix GSLMatrix::operator-(const GSLMatrix &rhs) const
+{
+  GSLMatrix tmp = *this;
+  tmp-=rhs;
+  return tmp;
+}
+
+const GSLMatrix GSLMatrix::operator*(const GSLMatrix &rhs) const
+{
+  gsl_matrix *res = gsl_matrix_alloc(rhs.rows, rhs.columns);
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix, rhs.matrix, 0.0, res);
+  return GSLMatrix(rhs.rows, rhs.columns, res);
+}
+
+/** Print matrix in a matlab style manner.
+ * \param &ost reference to output stream
+ * \param &mat reference to matrix to print
+ * \return reference to obtained output stream for concatenation
+ */
+ostream &operator<<(ostream &ost,const GSLMatrix &mat)
+{
+  ost << "[";
+  for(size_t i = 0;i<mat.rows;++i){
+    for(size_t j = 0; j<mat.columns;++j){
+      ost << mat.Get(i,j);
+      if(j!=mat.columns-1)
+        ost << " ";
+    }
+    if(i!=mat.rows-1)
+      ost << "; ";
+  }
+  ost << "]";
+  return ost;
+}
+
Index: src/LinearAlgebra/gslmatrix.hpp
===================================================================
--- src/LinearAlgebra/gslmatrix.hpp	(revision 9eb7580cd9f1b192263dddcbef45a27b39fc5eca)
+++ src/LinearAlgebra/gslmatrix.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
@@ -9,6 +9,4 @@
 #define GSLMATRIX_HPP_
 
-using namespace std;
-
 /*********************************************** includes ***********************************/
 
@@ -18,5 +16,8 @@
 #endif
 
+#include <iosfwd>
+
 #include <gsl/gsl_matrix.h>
+
 
 /****************************************** forward declarations *****************************/
@@ -28,16 +29,20 @@
 class GSLMatrix {
   friend class LinearSystemOfEquations;
+  friend std::ostream &operator<<(std::ostream &ost,const GSLMatrix &mat);
 
 public:
   GSLMatrix(size_t m, size_t n);
   GSLMatrix(const GSLMatrix * const src);
+  GSLMatrix(size_t m, size_t n, gsl_matrix *&src);
   ~GSLMatrix();
 
   // Accessing
   void SetFromDoubleArray(double * x);
-  double Get(size_t m, size_t n);
+  double Get(size_t m, size_t n) const;
   void Set(size_t m, size_t n, double x);
   double *Pointer(size_t m, size_t n);
-  const double *const_Pointer(size_t m, size_t n);
+  const double *const_Pointer(size_t m, size_t n) const;
+  size_t getRowCount() const;
+  size_t getColumnCount() const;
 
   // Initializing
@@ -50,13 +55,22 @@
   bool SwapColumns(size_t i, size_t j);
   bool SwapRowColumn(size_t i, size_t j);
+
+  // Properties
   bool Transpose();
-  bool IsNull();
-  bool IsPositive();
-  bool IsNegative();
-  bool IsNonNegative();
-  bool IsPositiveDefinite();
-  double Determinant();
+  bool IsNull() const;
+  bool IsPositive() const;
+  bool IsNegative() const;
+  bool IsNonNegative() const;
+  bool IsPositiveDefinite() const;
+  double Determinant() const;
 
-  GSLMatrix& operator=(const GSLMatrix& rhs);
+  // Operators
+  GSLMatrix &operator=(const GSLMatrix &src);
+  const GSLMatrix &operator+=(const GSLMatrix &rhs);
+  const GSLMatrix &operator-=(const GSLMatrix &rhs);
+  const GSLMatrix &operator*=(const GSLMatrix &rhs);
+  const GSLMatrix operator+(const GSLMatrix &rhs) const;
+  const GSLMatrix operator-(const GSLMatrix &rhs) const;
+  const GSLMatrix operator*(const GSLMatrix &rhs) const;
 
 private:
@@ -67,4 +81,5 @@
 };
 
+std::ostream &operator<<(std::ostream &ost,const GSLMatrix &mat);
 
 #endif /* GSLMATRIX_HPP_ */
