Index: src/LinearAlgebra/Makefile.am
===================================================================
--- src/LinearAlgebra/Makefile.am	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ src/LinearAlgebra/Makefile.am	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -18,7 +18,7 @@
   LineSegmentSet.cpp \
   linearsystemofequations.cpp \
-  Matrix.cpp \
   MatrixContent.cpp \
   Plane.cpp \
+  RealSpaceMatrix.cpp \
   Space.cpp \
   vector_ops.cpp \
@@ -34,6 +34,6 @@
   LineSegmentSet.hpp \
   linearsystemofequations.hpp \
-  Matrix.hpp \
   MatrixContent.hpp \
+  RealSpaceMatrix.hpp \
   Plane.hpp \
   Space.hpp \
Index: src/LinearAlgebra/Matrix.cpp
===================================================================
--- src/LinearAlgebra/Matrix.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ 	(revision )
@@ -1,332 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010 University of Bonn. All rights reserved.
- * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
- */
-
-/*
- * Matrix.cpp
- *
- *  Created on: Jun 25, 2010
- *      Author: crueger
- */
-
-// 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;
-
-Matrix::Matrix()
-{
-  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();
-}
-
-Matrix::~Matrix()
-{
-  // delete all views
-  for(int i=NDIM;i--;){
-    delete rows_ptr[i];
-  }
-  for(int i=NDIM;i--;){
-    delete columns_ptr[i];
-  }
-  delete diagonal_ptr;
-  delete content;
-}
-
-void Matrix::createViews(){
-  // create row views
-  for(int i=NDIM;i--;){
-    VectorContent *rowContent = new VectorViewContent(gsl_matrix_row(content->content,i));
-    rows_ptr[i] = new Vector(rowContent);
-  }
-  // create column views
-  for(int i=NDIM;i--;){
-    VectorContent *columnContent = new VectorViewContent(gsl_matrix_column(content->content,i));
-    columns_ptr[i] = new Vector(columnContent);
-  }
-  // create diagonal view
-  VectorContent *diagonalContent = new VectorViewContent(gsl_matrix_diagonal(content->content));
-  diagonal_ptr = new Vector(diagonalContent);
-}
-
-void Matrix::setIdentity(){
-  content->setIdentity();
-}
-
-void Matrix::setZero(){
-  content->setZero();
-}
-
-void Matrix::setRotation(const double x, const double y, const double z)
-{
-  set(0,0, cos(y)*cos(z));
-  set(0,1, cos(z)*sin(x)*sin(y) - cos(x)*sin(z));
-  set(0,2, cos(x)*cos(z)*sin(y) + sin(x) * sin(z));
-  set(1,0, cos(y)*sin(z));
-  set(1,1, cos(x)*cos(z) + sin(x)*sin(y)*sin(z));
-  set(1,2, -cos(z)*sin(x) + cos(x)*sin(y)*sin(z));
-  set(2,0, -sin(y));
-  set(2,1, cos(y)*sin(x));
-  set(2,2, cos(x)*cos(y));
-}
-
-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;
-  return tmp;
-}
-
-const Matrix Matrix::operator-(const Matrix &rhs) const
-{
-  Matrix tmp = *this;
-  tmp-=rhs;
-  return tmp;
-}
-
-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
-{
-  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
-  return *rows_ptr[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
-{
-  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
-  return *columns_ptr[i];
-}
-
-Vector &Matrix::diagonal()
-{
-  return *diagonal_ptr;
-}
-
-const Vector &Matrix::diagonal() const
-{
-  return *diagonal_ptr;
-}
-
-void Matrix::set(size_t i, size_t j, const double value)
-{
-  content->set(i,j, value);
-}
-
-double Matrix::determinant() const{
-  return at(0,0)*at(1,1)*at(2,2)
-       + at(0,1)*at(1,2)*at(2,0)
-       + at(0,2)*at(1,0)*at(2,1)
-       - at(2,0)*at(1,1)*at(0,2)
-       - at(2,1)*at(1,2)*at(0,0)
-       - at(2,2)*at(1,0)*at(0,1);
-}
-
-
-Matrix Matrix::invert() const{
-  double det = determinant();
-  if(fabs(det)<MYEPSILON){
-    throw NotInvertibleException(__FILE__,__LINE__);
-  }
-
-  double detReci = 1./det;
-  Matrix res;
-  res.set(0,0,  detReci*RDET2(at(1,1),at(2,1),at(1,2),at(2,2)));    // A_11
-  res.set(1,0, -detReci*RDET2(at(1,0),at(2,0),at(1,2),at(2,2)));    // A_21
-  res.set(2,0,  detReci*RDET2(at(1,0),at(2,0),at(1,1),at(2,1)));    // A_31
-  res.set(0,1, -detReci*RDET2(at(0,1),at(2,1),at(0,2),at(2,2)));    // A_12
-  res.set(1,1,  detReci*RDET2(at(0,0),at(2,0),at(0,2),at(2,2)));    // A_22
-  res.set(2,1, -detReci*RDET2(at(0,0),at(2,0),at(0,1),at(2,1)));    // A_32
-  res.set(0,2,  detReci*RDET2(at(0,1),at(1,1),at(0,2),at(1,2)));    // A_13
-  res.set(1,2, -detReci*RDET2(at(0,0),at(1,0),at(0,2),at(1,2)));    // A_23
-  res.set(2,2,  detReci*RDET2(at(0,0),at(1,0),at(0,1),at(1,1)));    // A_33
-  return res;
-}
-
-Matrix Matrix::transpose() const
-{
-  std::cout << "const Matrix::transpose()." << std::endl;
-  Matrix res = Matrix(const_cast<const MatrixContent *>(content)->transpose());
-  return res;
-}
-
-Matrix &Matrix::transpose()
-{
-  std::cout << "Matrix::transpose()." << std::endl;
-  content->transpose();
-  return *this;
-}
-
-Vector Matrix::transformToEigenbasis()
-{
-  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);
-  return evalues;
-}
-
-const Matrix &Matrix::operator*=(const double factor)
-    {
-  *content *= factor;
-  return *this;
-}
-
-const Matrix operator*(const double factor,const Matrix& mat)
-{
-  Matrix tmp = mat;
-  return tmp *= factor;
-}
-
-const Matrix operator*(const Matrix &mat,const double factor)
-{
-  return factor*mat;
-}
-
-bool Matrix::operator==(const Matrix &rhs) const
-{
-  return (*content == *(rhs.content));
-}
-
-/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
- * \param *symm 6-dim array of unique symmetric matrix components
- * \return allocated NDIM*NDIM array with the symmetric matrix
- */
-Matrix ReturnFullMatrixforSymmetric(const double * const symm)
-{
-  Matrix matrix;
-  matrix.set(0,0, symm[0]);
-  matrix.set(1,0, symm[1]);
-  matrix.set(2,0, symm[3]);
-  matrix.set(0,1, symm[1]);
-  matrix.set(1,1, symm[2]);
-  matrix.set(2,1, symm[4]);
-  matrix.set(0,2, symm[3]);
-  matrix.set(1,2, symm[4]);
-  matrix.set(2,2, symm[5]);
-  return matrix;
-};
-
-ostream &operator<<(ostream &ost,const Matrix &mat)
-{
-  for(int i = 0;i<NDIM;++i){
-    ost << "\n";
-    for(int j = 0; j<NDIM;++j){
-      ost << mat.at(i,j);
-      if(j!=NDIM-1)
-        ost << "; ";
-    }
-  }
-  return ost;
-}
-
-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 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ 	(revision )
@@ -1,173 +1,0 @@
-/*
- * Matrix.hpp
- *
- *  Created on: Jun 25, 2010
- *      Author: crueger
- */
-
-#ifndef MATRIX_HPP_
-#define MATRIX_HPP_
-
-#include <iosfwd>
-#include "defs.hpp"
-
-/**
- * Simple class to store matrices with a few extra functions
- */
-
-class Vector;
-class MatrixContent;
-
-class Matrix
-{
-  friend Vector operator*(const Matrix&,const Vector&);
-public:
-  Matrix();
-
-  /**
-   * construct a matrix from a 3x3 double array that contains all Elements.
-   *
-   * The elements are laid out in the following way:
-   * array -> matrix
-   * 0 -> (0,0)
-   * 1 -> (1,0)
-   * 2 -> (2,0)
-   * 3 -> (0,1)
-   * 4 -> (1,1)
-   * 5 -> (2,1)
-   * 6 -> (0,2)
-   * 7 -> (1,2)
-   * 8 -> (2,2)
-   *
-   */
-  Matrix(const double*);
-  Matrix(const Matrix&);
-  Matrix(const MatrixContent&);
-  virtual ~Matrix();
-
-  /**
-   * Set the matrix to a unit matrix.
-   */
-  void setIdentity();
-
-  /**
-   * Set all matrix entries to zero.
-   */
-  void setZero();
-
-  /**
-   * Sets all matrix corresponding to a rotation matrix,
-   * first around the x-, then the y-, finally the z-axis
-   * with given angles.
-   */
-  void setRotation(const double x, const double y, const double z);
-
-  /**
-   * Access the matrix at index (i,j)
-   */
-  double &at(size_t i, size_t j);
-  /**
-   * Access the matrix at index (i,j)
-   */
-  const double at(size_t i, size_t j) const;
-
-  /**
-   * Set the matrix at index (i,j).
-   *
-   * Slightly faster than at(i,j)=x
-   */
-  void set(size_t i, size_t j, const double value);
-
-  /**
-   * get the ith row of the matrix as a vector
-   */
-  Vector &row(size_t);
-  const Vector &row(size_t) const;
-
-  /**
-   * get the ith column of the matrix as a vector
-   */
-  Vector &column(size_t);
-  const Vector &column(size_t) const;
-
-  /**
-   * get the diagonal of the matrix as a vector
-   */
-  Vector &diagonal();
-  const Vector &diagonal() const;
-
-  /**
-   * Calculate the determinant of the matrix
-   */
-  double determinant() const;
-
-  /**
-   * Calculate the inverse of the matrix.
-   *
-   * Rather costly, so use precomputation as often as possible.
-   */
-  Matrix invert() const;
-
-  /**
-   * Diagonalizes a matrix and sets its rows to the resulting eigenvalues.
-   * The eigenvalues are returned as a vector.
-   *
-   * Rather costly, so use precomputation as often as possible.
-   */
-  Vector transformToEigenbasis();
-
-  /**
-   * Calculate the transpose of the matrix.
-   */
-  Matrix transpose() const;
-  Matrix &transpose();
-
-  // operators
-  Matrix &operator=(const Matrix&);
-
-  const Matrix &operator+=(const Matrix&);
-  const Matrix &operator-=(const Matrix&);
-  const Matrix &operator*=(const Matrix&);
-
-  const Matrix &operator*=(const double);
-
-  const Matrix operator+(const Matrix&) const;
-  const Matrix operator-(const Matrix&) const;
-  const Matrix operator*(const Matrix&) const;
-
-  bool operator==(const Matrix&) const;
-
-private:
-  Matrix(MatrixContent*);
-  void createViews();
-  MatrixContent *content;
-  // we keep around some Vector views of the matrix, to return references
-  Vector* rows_ptr[NDIM];
-  Vector* columns_ptr[NDIM];
-  Vector* diagonal_ptr;
-};
-
-const Matrix operator*(const double,const Matrix&);
-const Matrix operator*(const Matrix&,const double);
-
-/**
- * Takes a symmetric matrix that stores the lower diagonal and produces a
- * full matrix.
- *
- * The array is laid out as follows:
- *
- * array -> matrix
- * 0 -> (0,0)
- * 1 -> (1,0);(0,1)
- * 2 -> (1,1)
- * 3 -> (2,0);(0,2)
- * 4 -> (2,1);(1,2)
- * 5 -> (2,2)
- */
-Matrix ReturnFullMatrixforSymmetric(const double * const cell_size);
-
-std::ostream &operator<<(std::ostream&,const Matrix&);
-Vector operator*(const Matrix&,const Vector&);
-Vector& operator*=(Vector&,const Matrix&);
-
-#endif /* MATRIX_HPP_ */
Index: src/LinearAlgebra/MatrixContent.cpp
===================================================================
--- src/LinearAlgebra/MatrixContent.cpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ src/LinearAlgebra/MatrixContent.cpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -14,5 +14,5 @@
 #include "Helpers/MemDebug.hpp"
 
-#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
 #include "Helpers/Assert.hpp"
 #include "Exceptions/NotInvertibleException.hpp"
Index: src/LinearAlgebra/MatrixContent.hpp
===================================================================
--- src/LinearAlgebra/MatrixContent.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ src/LinearAlgebra/MatrixContent.hpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -40,5 +40,5 @@
 {
   friend Vector operator*(const MatrixContent &mat,const Vector &vec);
-  friend class Matrix;
+  friend class RealSpaceMatrix;
 public:
   MatrixContent(size_t rows, size_t columns);
Index: src/LinearAlgebra/RealSpaceMatrix.cpp
===================================================================
--- src/LinearAlgebra/RealSpaceMatrix.cpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
+++ src/LinearAlgebra/RealSpaceMatrix.cpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -0,0 +1,332 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * RealSpaceMatrix.cpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Helpers/MemDebug.hpp"
+
+#include "LinearAlgebra/RealSpaceMatrix.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;
+
+RealSpaceMatrix::RealSpaceMatrix()
+{
+  content = new MatrixContent(NDIM, NDIM);
+  createViews();
+}
+
+RealSpaceMatrix::RealSpaceMatrix(const double* src)
+{
+  content = new MatrixContent(NDIM, NDIM, src);
+  createViews();
+}
+
+RealSpaceMatrix::RealSpaceMatrix(const RealSpaceMatrix &src)
+{
+  content = new MatrixContent(src.content);
+  createViews();
+}
+
+RealSpaceMatrix::RealSpaceMatrix(const MatrixContent &src)
+{
+  content = new MatrixContent(src);
+  createViews();
+}
+
+RealSpaceMatrix::RealSpaceMatrix(MatrixContent* _content)
+{
+  content = new MatrixContent(_content);
+  createViews();
+}
+
+RealSpaceMatrix::~RealSpaceMatrix()
+{
+  // delete all views
+  for(int i=NDIM;i--;){
+    delete rows_ptr[i];
+  }
+  for(int i=NDIM;i--;){
+    delete columns_ptr[i];
+  }
+  delete diagonal_ptr;
+  delete content;
+}
+
+void RealSpaceMatrix::createViews(){
+  // create row views
+  for(int i=NDIM;i--;){
+    VectorContent *rowContent = new VectorViewContent(gsl_matrix_row(content->content,i));
+    rows_ptr[i] = new Vector(rowContent);
+  }
+  // create column views
+  for(int i=NDIM;i--;){
+    VectorContent *columnContent = new VectorViewContent(gsl_matrix_column(content->content,i));
+    columns_ptr[i] = new Vector(columnContent);
+  }
+  // create diagonal view
+  VectorContent *diagonalContent = new VectorViewContent(gsl_matrix_diagonal(content->content));
+  diagonal_ptr = new Vector(diagonalContent);
+}
+
+void RealSpaceMatrix::setIdentity(){
+  content->setIdentity();
+}
+
+void RealSpaceMatrix::setZero(){
+  content->setZero();
+}
+
+void RealSpaceMatrix::setRotation(const double x, const double y, const double z)
+{
+  set(0,0, cos(y)*cos(z));
+  set(0,1, cos(z)*sin(x)*sin(y) - cos(x)*sin(z));
+  set(0,2, cos(x)*cos(z)*sin(y) + sin(x) * sin(z));
+  set(1,0, cos(y)*sin(z));
+  set(1,1, cos(x)*cos(z) + sin(x)*sin(y)*sin(z));
+  set(1,2, -cos(z)*sin(x) + cos(x)*sin(y)*sin(z));
+  set(2,0, -sin(y));
+  set(2,1, cos(y)*sin(x));
+  set(2,2, cos(x)*cos(y));
+}
+
+RealSpaceMatrix &RealSpaceMatrix::operator=(const RealSpaceMatrix &src)
+{
+  // MatrixContent checks for self-assignment
+  *content = *(src.content);
+  return *this;
+}
+
+const RealSpaceMatrix &RealSpaceMatrix::operator+=(const RealSpaceMatrix &rhs)
+{
+  *content += *(rhs.content);
+  return *this;
+}
+
+const RealSpaceMatrix &RealSpaceMatrix::operator-=(const RealSpaceMatrix &rhs)
+{
+  *content -= *(rhs.content);
+  return *this;
+}
+
+const RealSpaceMatrix &RealSpaceMatrix::operator*=(const RealSpaceMatrix &rhs)
+{
+  (*content) *= (*rhs.content);
+  return *this;
+}
+
+const RealSpaceMatrix RealSpaceMatrix::operator+(const RealSpaceMatrix &rhs) const
+{
+  RealSpaceMatrix tmp = *this;
+  tmp+=rhs;
+  return tmp;
+}
+
+const RealSpaceMatrix RealSpaceMatrix::operator-(const RealSpaceMatrix &rhs) const
+{
+  RealSpaceMatrix tmp = *this;
+  tmp-=rhs;
+  return tmp;
+}
+
+const RealSpaceMatrix RealSpaceMatrix::operator*(const RealSpaceMatrix &rhs) const
+{
+  RealSpaceMatrix tmp(content);
+  tmp *= rhs;
+  return tmp;
+}
+
+double &RealSpaceMatrix::at(size_t i, size_t j)
+{
+  return content->at(i,j);
+}
+
+const double RealSpaceMatrix::at(size_t i, size_t j) const
+{
+  return content->at(i,j);
+}
+
+Vector &RealSpaceMatrix::row(size_t i)
+{
+  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
+  return *rows_ptr[i];
+}
+
+const Vector &RealSpaceMatrix::row(size_t i) const
+{
+  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
+  return *rows_ptr[i];
+}
+
+Vector &RealSpaceMatrix::column(size_t i)
+{
+  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
+  return *columns_ptr[i];
+}
+
+const Vector &RealSpaceMatrix::column(size_t i) const
+{
+  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
+  return *columns_ptr[i];
+}
+
+Vector &RealSpaceMatrix::diagonal()
+{
+  return *diagonal_ptr;
+}
+
+const Vector &RealSpaceMatrix::diagonal() const
+{
+  return *diagonal_ptr;
+}
+
+void RealSpaceMatrix::set(size_t i, size_t j, const double value)
+{
+  content->set(i,j, value);
+}
+
+double RealSpaceMatrix::determinant() const{
+  return at(0,0)*at(1,1)*at(2,2)
+       + at(0,1)*at(1,2)*at(2,0)
+       + at(0,2)*at(1,0)*at(2,1)
+       - at(2,0)*at(1,1)*at(0,2)
+       - at(2,1)*at(1,2)*at(0,0)
+       - at(2,2)*at(1,0)*at(0,1);
+}
+
+
+RealSpaceMatrix RealSpaceMatrix::invert() const{
+  double det = determinant();
+  if(fabs(det)<MYEPSILON){
+    throw NotInvertibleException(__FILE__,__LINE__);
+  }
+
+  double detReci = 1./det;
+  RealSpaceMatrix res;
+  res.set(0,0,  detReci*RDET2(at(1,1),at(2,1),at(1,2),at(2,2)));    // A_11
+  res.set(1,0, -detReci*RDET2(at(1,0),at(2,0),at(1,2),at(2,2)));    // A_21
+  res.set(2,0,  detReci*RDET2(at(1,0),at(2,0),at(1,1),at(2,1)));    // A_31
+  res.set(0,1, -detReci*RDET2(at(0,1),at(2,1),at(0,2),at(2,2)));    // A_12
+  res.set(1,1,  detReci*RDET2(at(0,0),at(2,0),at(0,2),at(2,2)));    // A_22
+  res.set(2,1, -detReci*RDET2(at(0,0),at(2,0),at(0,1),at(2,1)));    // A_32
+  res.set(0,2,  detReci*RDET2(at(0,1),at(1,1),at(0,2),at(1,2)));    // A_13
+  res.set(1,2, -detReci*RDET2(at(0,0),at(1,0),at(0,2),at(1,2)));    // A_23
+  res.set(2,2,  detReci*RDET2(at(0,0),at(1,0),at(0,1),at(1,1)));    // A_33
+  return res;
+}
+
+RealSpaceMatrix RealSpaceMatrix::transpose() const
+{
+  std::cout << "const Matrix::transpose()." << std::endl;
+  RealSpaceMatrix res = RealSpaceMatrix(const_cast<const MatrixContent *>(content)->transpose());
+  return res;
+}
+
+RealSpaceMatrix &RealSpaceMatrix::transpose()
+{
+  std::cout << "Matrix::transpose()." << std::endl;
+  content->transpose();
+  return *this;
+}
+
+Vector RealSpaceMatrix::transformToEigenbasis()
+{
+  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);
+  return evalues;
+}
+
+const RealSpaceMatrix &RealSpaceMatrix::operator*=(const double factor)
+    {
+  *content *= factor;
+  return *this;
+}
+
+const RealSpaceMatrix operator*(const double factor,const RealSpaceMatrix& mat)
+{
+  RealSpaceMatrix tmp = mat;
+  return tmp *= factor;
+}
+
+const RealSpaceMatrix operator*(const RealSpaceMatrix &mat,const double factor)
+{
+  return factor*mat;
+}
+
+bool RealSpaceMatrix::operator==(const RealSpaceMatrix &rhs) const
+{
+  return (*content == *(rhs.content));
+}
+
+/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
+ * \param *symm 6-dim array of unique symmetric matrix components
+ * \return allocated NDIM*NDIM array with the symmetric matrix
+ */
+RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const symm)
+{
+  RealSpaceMatrix matrix;
+  matrix.set(0,0, symm[0]);
+  matrix.set(1,0, symm[1]);
+  matrix.set(2,0, symm[3]);
+  matrix.set(0,1, symm[1]);
+  matrix.set(1,1, symm[2]);
+  matrix.set(2,1, symm[4]);
+  matrix.set(0,2, symm[3]);
+  matrix.set(1,2, symm[4]);
+  matrix.set(2,2, symm[5]);
+  return matrix;
+};
+
+ostream &operator<<(ostream &ost,const RealSpaceMatrix &mat)
+{
+  for(int i = 0;i<NDIM;++i){
+    ost << "\n";
+    for(int j = 0; j<NDIM;++j){
+      ost << mat.at(i,j);
+      if(j!=NDIM-1)
+        ost << "; ";
+    }
+  }
+  return ost;
+}
+
+Vector operator*(const RealSpaceMatrix &mat,const Vector &vec)
+{
+  return (*mat.content) * vec;
+}
+
+Vector &operator*=(Vector& lhs,const RealSpaceMatrix &mat)
+{
+  lhs = mat*lhs;
+  return lhs;
+}
+
Index: src/LinearAlgebra/RealSpaceMatrix.hpp
===================================================================
--- src/LinearAlgebra/RealSpaceMatrix.hpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
+++ src/LinearAlgebra/RealSpaceMatrix.hpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -0,0 +1,177 @@
+/*
+ * RealSpaceMatrix.hpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+#ifndef REALSPACEMATRIX_HPP_
+#define REALSPACEMATRIX_HPP_
+
+#include <iosfwd>
+#include "defs.hpp"
+
+/**
+ * Simple class to store matrices with a few extra functions
+ */
+
+class Vector;
+class MatrixContent;
+
+/** 3x3 Matrix class.
+ * This class has some specific features for 3D space such as rotations or
+ * hard-coded determinants.
+ */
+class RealSpaceMatrix
+{
+  friend Vector operator*(const RealSpaceMatrix&,const Vector&);
+public:
+  RealSpaceMatrix();
+
+  /**
+   * construct a matrix from a 3x3 double array that contains all Elements.
+   *
+   * The elements are laid out in the following way:
+   * array -> matrix
+   * 0 -> (0,0)
+   * 1 -> (1,0)
+   * 2 -> (2,0)
+   * 3 -> (0,1)
+   * 4 -> (1,1)
+   * 5 -> (2,1)
+   * 6 -> (0,2)
+   * 7 -> (1,2)
+   * 8 -> (2,2)
+   *
+   */
+  RealSpaceMatrix(const double*);
+  RealSpaceMatrix(const RealSpaceMatrix&);
+  RealSpaceMatrix(const MatrixContent&);
+  virtual ~RealSpaceMatrix();
+
+  /**
+   * Set the matrix to a unit matrix.
+   */
+  void setIdentity();
+
+  /**
+   * Set all matrix entries to zero.
+   */
+  void setZero();
+
+  /**
+   * Sets all matrix corresponding to a rotation matrix,
+   * first around the x-, then the y-, finally the z-axis
+   * with given angles.
+   */
+  void setRotation(const double x, const double y, const double z);
+
+  /**
+   * Access the matrix at index (i,j)
+   */
+  double &at(size_t i, size_t j);
+  /**
+   * Access the matrix at index (i,j)
+   */
+  const double at(size_t i, size_t j) const;
+
+  /**
+   * Set the matrix at index (i,j).
+   *
+   * Slightly faster than at(i,j)=x
+   */
+  void set(size_t i, size_t j, const double value);
+
+  /**
+   * get the ith row of the matrix as a vector
+   */
+  Vector &row(size_t);
+  const Vector &row(size_t) const;
+
+  /**
+   * get the ith column of the matrix as a vector
+   */
+  Vector &column(size_t);
+  const Vector &column(size_t) const;
+
+  /**
+   * get the diagonal of the matrix as a vector
+   */
+  Vector &diagonal();
+  const Vector &diagonal() const;
+
+  /**
+   * Calculate the determinant of the matrix
+   */
+  double determinant() const;
+
+  /**
+   * Calculate the inverse of the matrix.
+   *
+   * Rather costly, so use precomputation as often as possible.
+   */
+  RealSpaceMatrix invert() const;
+
+  /**
+   * Diagonalizes a matrix and sets its rows to the resulting eigenvalues.
+   * The eigenvalues are returned as a vector.
+   *
+   * Rather costly, so use precomputation as often as possible.
+   */
+  Vector transformToEigenbasis();
+
+  /**
+   * Calculate the transpose of the matrix.
+   */
+  RealSpaceMatrix transpose() const;
+  RealSpaceMatrix &transpose();
+
+  // operators
+  RealSpaceMatrix &operator=(const RealSpaceMatrix&);
+
+  const RealSpaceMatrix &operator+=(const RealSpaceMatrix&);
+  const RealSpaceMatrix &operator-=(const RealSpaceMatrix&);
+  const RealSpaceMatrix &operator*=(const RealSpaceMatrix&);
+
+  const RealSpaceMatrix &operator*=(const double);
+
+  const RealSpaceMatrix operator+(const RealSpaceMatrix&) const;
+  const RealSpaceMatrix operator-(const RealSpaceMatrix&) const;
+  const RealSpaceMatrix operator*(const RealSpaceMatrix&) const;
+
+  bool operator==(const RealSpaceMatrix&) const;
+
+private:
+  RealSpaceMatrix(MatrixContent*);
+  void createViews();
+  MatrixContent *content;
+  // we keep around some Vector views of the matrix, to return references
+  Vector* rows_ptr[NDIM];
+  Vector* columns_ptr[NDIM];
+  Vector* diagonal_ptr;
+};
+
+const RealSpaceMatrix operator*(const double,const RealSpaceMatrix&);
+const RealSpaceMatrix operator*(const RealSpaceMatrix&,const double);
+
+/**
+ * Takes a symmetric matrix that stores the lower diagonal and produces a
+ * full matrix.
+ *
+ * The array is laid out as follows:
+ *
+ * array -> matrix
+ * 0 -> (0,0)
+ * 1 -> (1,0);(0,1)
+ * 2 -> (1,1)
+ * 3 -> (2,0);(0,2)
+ * 4 -> (2,1);(1,2)
+ * 5 -> (2,2)
+ */
+RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const cell_size);
+
+std::ostream &operator<<(std::ostream&,const RealSpaceMatrix&);
+Vector operator*(const RealSpaceMatrix&,const Vector&);
+Vector& operator*=(Vector&,const RealSpaceMatrix&);
+
+#endif /* REALSPACEMATRIX_HPP_ */
Index: src/LinearAlgebra/Vector.hpp
===================================================================
--- src/LinearAlgebra/Vector.hpp	(revision 3bc9264606122a0c59b0d7bf20669cc353efbbec)
+++ src/LinearAlgebra/Vector.hpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -22,5 +22,5 @@
 
 class Vector;
-class Matrix;
+class RealSpaceMatrix;
 class MatrixContent;
 struct VectorContent;
@@ -33,5 +33,5 @@
 class Vector : public Space{
   friend Vector operator*(const MatrixContent&,const Vector&);
-  friend class Matrix;
+  friend class RealSpaceMatrix;
 public:
   Vector();
