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;
