Index: src/Actions/WorldAction/RepeatBoxAction.cpp
===================================================================
--- src/Actions/WorldAction/RepeatBoxAction.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/Actions/WorldAction/RepeatBoxAction.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -13,4 +13,5 @@
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
@@ -60,5 +61,7 @@
     (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
     double * const cell_size = World::getInstance().getDomain();
-    double *M = ReturnFullMatrixforSymmetric(cell_size);
+    double *M_double = ReturnFullMatrixforSymmetric(cell_size);
+    Matrix M = Matrix(M_double);
+    delete[] M_double;
     Vector x,y;
     int n[NDIM];
@@ -118,5 +121,4 @@
       }
     }
-    delete(M);
     delete dialog;
     return Action::success;
Index: src/Exceptions/NotInvertibleException.cpp
===================================================================
--- src/Exceptions/NotInvertibleException.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
+++ src/Exceptions/NotInvertibleException.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -0,0 +1,21 @@
+/*
+ * NotInvertibleException.cpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+#include "NotInvertibleException.hpp"
+
+using namespace std;
+
+NotInvertibleException::NotInvertibleException(string _file, int _line) :
+  MathException(_file,_line)
+{}
+
+NotInvertibleException::~NotInvertibleException() throw ()
+{}
+
+const char* NotInvertibleException::what() const throw(){
+  return "Given Matrix is not invertible";
+}
Index: src/Exceptions/NotInvertibleException.hpp
===================================================================
--- src/Exceptions/NotInvertibleException.hpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
+++ src/Exceptions/NotInvertibleException.hpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -0,0 +1,22 @@
+/*
+ * NotInvertibleException.hpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+#ifndef NOTINVERTIBLEEXCEPTION_HPP_
+#define NOTINVERTIBLEEXCEPTION_HPP_
+
+#include "Exceptions/MathException.hpp"
+
+class NotInvertibleException : public MathException
+{
+public:
+  NotInvertibleException(std::string _file, int _line);
+  virtual ~NotInvertibleException() throw ();
+
+  virtual const char* what() const throw();
+};
+
+#endif /* NOTINVERTIBLEEXCEPTION_HPP_ */
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/Makefile.am	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -30,15 +30,36 @@
   atom_trajectoryparticleinfo.hpp
 
+EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
+				  Exceptions/LinearDependenceException.cpp \
+				  Exceptions/MathException.cpp \
+				  Exceptions/NotInvertibleException.cpp \
+				  Exceptions/SkewException.cpp \
+				  Exceptions/ZeroVectorException.cpp
+				  
+EXCEPTIONHEADER = Exceptions/CustomException.hpp \
+				  Exceptions/LinearDependenceException.hpp \
+				  Exceptions/MathException.hpp \
+				  Exceptions/NotInvertibleException.hpp \
+				  Exceptions/SkewException.hpp \
+				  Exceptions/ZeroVectorException.hpp
+
+# TODO: Move Exceptionsource back down, when transition is done
 LINALGSOURCE = \
+  ${EXCEPTIONSOURCE} \
   ${HELPERSOURCE} \
   gslmatrix.cpp \
   gslvector.cpp \
   linearsystemofequations.cpp \
+  Matrix.cpp \
   Space.cpp \
   vector.cpp
+  
 			   
-LINALGHEADER = gslmatrix.hpp \
+LINALGHEADER = \
+  ${EXCEPTIONHEADER} \
+  gslmatrix.hpp \
   gslvector.hpp \
   linearsystemofequations.hpp \
+  Matrix.hpp \
   Space.hpp \
   vector.hpp
@@ -118,17 +139,6 @@
   Descriptors/MoleculeNameDescriptor.hpp \
   Descriptors/MoleculePtrDescriptor.hpp
-				   
-EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
-				  Exceptions/LinearDependenceException.cpp \
-				  Exceptions/MathException.cpp \
-				  Exceptions/SkewException.cpp \
-				  Exceptions/ZeroVectorException.cpp
-				  
-EXCEPTIONHEADER = Exceptions/CustomException.hpp \
-				  Exceptions/LinearDependenceException.hpp \
-				  Exceptions/MathException.hpp \
-				  Exceptions/SkewException.hpp \
-				  Exceptions/ZeroVectorException.hpp
-
+
+#TODO: Exceptionsource should go here, when the transisition makes this possible				   
 SOURCE = \
   ${ANALYSISSOURCE} \
@@ -139,5 +149,4 @@
   ${DESCRIPTORSOURCE} \
   ${HELPERSOURCE} \
-  ${EXCEPTIONSOURCE} \
   bond.cpp \
   bondgraph.cpp \
@@ -185,5 +194,4 @@
   ${PATTERNHEADER} \
   ${DESCRIPTORHEADER} \
-  ${EXCEPTIONHEADER} \
   bond.hpp \
   bondgraph.hpp \
Index: src/Matrix.cpp
===================================================================
--- src/Matrix.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
+++ src/Matrix.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -0,0 +1,136 @@
+/*
+ * Matrix.cpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+#include "Matrix.hpp"
+#include "Helpers/Assert.hpp"
+#include "Exceptions/NotInvertibleException.hpp"
+#include "Helpers/fast_functions.hpp"
+
+#include <gsl/gsl_blas.h>
+#include <cmath>
+
+Matrix::Matrix()
+{
+  content = gsl_matrix_calloc(NDIM, NDIM);
+}
+
+Matrix::Matrix(const double* src){
+  content = gsl_matrix_alloc(NDIM, NDIM);
+  this->at(0,0) = src[0];
+  this->at(1,0) = src[1];
+  this->at(2,0) = src[2];
+
+  this->at(0,1) = src[3];
+  this->at(1,1) = src[4];
+  this->at(2,1) = src[5];
+
+  this->at(0,2) = src[6];
+  this->at(1,2) = src[7];
+  this->at(2,2) = src[8];
+}
+
+Matrix::Matrix(const Matrix &src){
+  content = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_matrix_memcpy(content,src.content);
+}
+
+Matrix::Matrix(gsl_matrix* _content) :
+  content(_content)
+{}
+
+Matrix::~Matrix()
+{
+  gsl_matrix_free(content);
+}
+
+Matrix &Matrix::operator=(const Matrix &src){
+  if(&src!=this){
+    gsl_matrix_memcpy(content,src.content);
+  }
+  return *this;
+}
+
+Matrix &Matrix::operator+=(const Matrix &rhs){
+  gsl_matrix_add(content, rhs.content);
+  return *this;
+}
+
+Matrix &Matrix::operator-=(const Matrix &rhs){
+  gsl_matrix_sub(content, rhs.content);
+  return *this;
+}
+
+Matrix &Matrix::operator*=(const Matrix &rhs){
+  (*this) = (*this)*rhs;
+  return *this;
+}
+
+Matrix Matrix::operator+(const Matrix &rhs) const{
+  Matrix tmp = *this;
+  tmp+=rhs;
+  return tmp;
+}
+
+Matrix Matrix::operator-(const Matrix &rhs) const{
+  Matrix tmp = *this;
+  tmp-=rhs;
+  return tmp;
+}
+
+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);
+}
+
+double &Matrix::at(size_t i, size_t j){
+  return *gsl_matrix_ptr (content, i, j);
+}
+
+double Matrix::determinant(){
+  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(){
+  double det = determinant();
+  if(fabs(det)<MYEPSILON){
+    throw NotInvertibleException(__FILE__,__LINE__);
+  }
+
+  double detReci = 1./det;
+  Matrix res;
+  res.at(0,0) =  detReci*RDET2(at(1,1),at(2,1),at(1,2),at(2,2));    // A_11
+  res.at(1,0) = -detReci*RDET2(at(1,0),at(2,0),at(1,2),at(2,2));    // A_21
+  res.at(2,0) =  detReci*RDET2(at(1,0),at(2,0),at(1,1),at(2,1));    // A_31
+  res.at(0,1) = -detReci*RDET2(at(0,1),at(2,1),at(0,2),at(2,2));    // A_12
+  res.at(1,1) =  detReci*RDET2(at(0,0),at(2,0),at(0,2),at(2,2));    // A_22
+  res.at(2,1) = -detReci*RDET2(at(0,0),at(2,0),at(0,1),at(2,1));    // A_32
+  res.at(0,2) =  detReci*RDET2(at(0,1),at(1,1),at(0,2),at(1,2));    // A_13
+  res.at(1,2) = -detReci*RDET2(at(0,0),at(1,0),at(0,2),at(1,2));    // A_23
+  res.at(2,2) =  detReci*RDET2(at(0,0),at(1,0),at(0,1),at(1,1));    // A_33
+  return res;
+}
+
+Matrix &Matrix::operator*=(const double factor){
+  gsl_matrix_scale(content, factor);
+  return *this;
+}
+
+Matrix operator*(const double factor,const Matrix& mat){
+  Matrix tmp = mat;
+  tmp*=factor;
+  return tmp;
+}
+
+Matrix operator*(const Matrix &mat,const double factor){
+  return factor*mat;
+}
Index: src/Matrix.hpp
===================================================================
--- src/Matrix.hpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
+++ src/Matrix.hpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -0,0 +1,52 @@
+/*
+ * Matrix.hpp
+ *
+ *  Created on: Jun 25, 2010
+ *      Author: crueger
+ */
+
+#ifndef MATRIX_HPP_
+#define MATRIX_HPP_
+
+#include <gsl/gsl_matrix.h>
+#include "defs.hpp"
+
+class Vector;
+
+class Matrix
+{
+  friend Vector operator*(const Matrix&,const Vector&);
+public:
+  Matrix();
+  Matrix(const double*);
+  Matrix(const Matrix&);
+  virtual ~Matrix();
+
+  double &at(size_t i, size_t j);
+
+  double determinant();
+
+  Matrix invert();
+
+  // operators
+  Matrix &operator=(const Matrix&);
+
+  Matrix &operator+=(const Matrix&);
+  Matrix &operator-=(const Matrix&);
+  Matrix &operator*=(const Matrix&);
+
+  Matrix &operator*=(const double);
+
+  Matrix operator+(const Matrix&) const;
+  Matrix operator-(const Matrix&) const;
+  Matrix operator*(const Matrix&) const;
+
+private:
+  Matrix(gsl_matrix*);
+  gsl_matrix *content;
+};
+
+Matrix operator*(const double,const Matrix&);
+Matrix operator*(const Matrix&,const double);
+
+#endif /* MATRIX_HPP_ */
Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/analysis_correlation.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -19,4 +19,5 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
@@ -135,6 +136,8 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
+      Matrix FullMatrix = Matrix(FullMatrix_double);
+      Matrix FullInverseMatrix = FullMatrix.invert();
+      delete[](FullMatrix_double);
       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
@@ -176,6 +179,4 @@
             }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
   }
@@ -246,6 +247,8 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
+      Matrix FullMatrix = Matrix(FullMatrix_double);
+      Matrix FullInverseMatrix = FullMatrix.invert();
+      delete[] FullMatrix_double;
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
@@ -267,6 +270,4 @@
           }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
 
@@ -352,6 +353,8 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
+      Matrix FullMatrix = Matrix(FullMatrix_double);
+      Matrix FullInverseMatrix = FullMatrix.invert();
+      delete[](FullMatrix_double);
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
@@ -381,6 +384,4 @@
           }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
 
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/boundary.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -22,4 +22,5 @@
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
 
 #include<gsl/gsl_poly.h>
@@ -764,7 +765,9 @@
   int N[NDIM];
   int n[NDIM];
-  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-  double Rotations[NDIM*NDIM];
-  double *MInverse = InverseMatrix(M);
+  double *M_double =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
+  Matrix M = Matrix(M_double);
+  delete[](M_double);
+  Matrix Rotations;
+  Matrix MInverse = M.invert();
   Vector AtomTranslations;
   Vector FillerTranslations;
@@ -800,5 +803,5 @@
   // calculate filler grid in [0,1]^3
   FillerDistance = Vector(distance[0], distance[1], distance[2]);
-  FillerDistance.InverseMatrixMultiplication(M);
+  FillerDistance.MatrixMultiplication(MInverse);
   for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance[i]);
@@ -835,13 +838,13 @@
             }
 
-            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
-            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
-            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
-            Rotations[7] =               sin(phi[1])                                                ;
-            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            Rotations.at(0,0) =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
+            Rotations.at(0,1) =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
+            Rotations.at(0,2) =               cos(phi[1])*sin(phi[2])                                     ;
+            Rotations.at(1,0) = - sin(phi[0])*cos(phi[1])                                                ;
+            Rotations.at(1,1) =   cos(phi[0])*cos(phi[1])                                                ;
+            Rotations.at(1,2) =               sin(phi[1])                                                ;
+            Rotations.at(2,0) = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
+            Rotations.at(2,1) = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
+            Rotations.at(2,2) =               cos(phi[1])*cos(phi[2])                                     ;
           }
 
@@ -892,6 +895,4 @@
             }
       }
-  delete[](M);
-  delete[](MInverse);
 
   return Filling;
Index: src/ellipsoid.cpp
===================================================================
--- src/ellipsoid.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/ellipsoid.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -21,4 +21,5 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 
@@ -34,5 +35,5 @@
   Vector helper, RefPoint;
   double distance = -1.;
-  double Matrix[NDIM*NDIM];
+  Matrix Matrix;
   double InverseLength[3];
   double psi,theta,phi; // euler angles in ZX'Z'' convention
@@ -51,13 +52,13 @@
   theta = EllipsoidAngle[1];
   phi = EllipsoidAngle[2];
-  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
-  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
-  Matrix[2] = sin(psi)*sin(theta);
-  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
-  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
-  Matrix[5] = -cos(psi)*sin(theta);
-  Matrix[6] = sin(theta)*sin(phi);
-  Matrix[7] = sin(theta)*cos(phi);
-  Matrix[8] = cos(theta);
+  Matrix.at(0,0) = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
+  Matrix.at(1,0) = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
+  Matrix.at(2,0) = sin(psi)*sin(theta);
+  Matrix.at(0,1) = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
+  Matrix.at(1,1) = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
+  Matrix.at(2,1) = -cos(psi)*sin(theta);
+  Matrix.at(0,2) = sin(theta)*sin(phi);
+  Matrix.at(1,2) = sin(theta)*cos(phi);
+  Matrix.at(2,2) = cos(theta);
   helper.MatrixMultiplication(Matrix);
   helper.ScaleAll(InverseLength);
@@ -73,13 +74,13 @@
   phi = -EllipsoidAngle[2];
   helper.ScaleAll(EllipsoidLength);
-  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
-  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
-  Matrix[2] = sin(psi)*sin(theta);
-  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
-  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
-  Matrix[5] = -cos(psi)*sin(theta);
-  Matrix[6] = sin(theta)*sin(phi);
-  Matrix[7] = sin(theta)*cos(phi);
-  Matrix[8] = cos(theta);
+  Matrix.at(0,0) = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
+  Matrix.at(1,0) = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
+  Matrix.at(2,0) = sin(psi)*sin(theta);
+  Matrix.at(0,1) = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
+  Matrix.at(1,1) = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
+  Matrix.at(2,1) = -cos(psi)*sin(theta);
+  Matrix.at(0,2) = sin(theta)*sin(phi);
+  Matrix.at(1,2) = sin(theta)*cos(phi);
+  Matrix.at(2,2) = cos(theta);
   helper.MatrixMultiplication(Matrix);
   //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/molecule.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -27,4 +27,5 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "World.hpp"
 #include "Plane.hpp"
@@ -284,5 +285,5 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
-  double *matrix = NULL;
+  Matrix matrix;
   bond *Binder = NULL;
   double * const cell_size = World::getInstance().getDomain();
@@ -307,8 +308,9 @@
       } // (signs are correct, was tested!)
     }
-    matrix = ReturnFullMatrixforSymmetric(cell_size);
+    double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
+    matrix = Matrix(matrix_double);
+    delete[](matrix_double);
     Orthovector1.MatrixMultiplication(matrix);
     InBondvector -= Orthovector1; // subtract just the additional translation
-    delete[](matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -541,5 +543,4 @@
       break;
   }
-  delete[](matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
Index: src/molecule_fragmentation.cpp
===================================================================
--- src/molecule_fragmentation.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/molecule_fragmentation.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -22,4 +22,5 @@
 #include "periodentafel.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -1709,5 +1710,7 @@
   atom *OtherWalker = NULL;
   double * const cell_size = World::getInstance().getDomain();
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
+  double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
+  Matrix matrix = Matrix(matrix_double);
+  delete[](matrix_double);
   enum Shading *ColorList = NULL;
   double tmp;
@@ -1780,5 +1783,4 @@
   delete(AtomStack);
   delete[](ColorList);
-  delete[](matrix);
   DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
 };
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/molecule_geometry.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -19,4 +19,5 @@
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
 #include <boost/foreach.hpp>
 
@@ -154,7 +155,8 @@
 
   const double *cell_size = World::getInstance().getDomain();
-  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *M_double = ReturnFullMatrixforSymmetric(cell_size);
+  Matrix M = Matrix(M_double);
+  delete[](M_double);
   a->MatrixMultiplication(M);
-  delete[](M);
 
   return a;
@@ -274,6 +276,8 @@
 {
   double * const cell_size = World::getInstance().getDomain();
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
-  double *inversematrix = InverseMatrix(matrix);
+  double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
+  Matrix matrix = Matrix(matrix_double);
+  delete[](matrix_double);
+  Matrix inversematrix = matrix.invert();
   double tmp;
   bool flag;
@@ -324,6 +328,4 @@
     }
   } while (!flag);
-  delete[](matrix);
-  delete[](inversematrix);
 
   Center.Scale(1./static_cast<double>(getAtomCount()));
@@ -387,5 +389,6 @@
     DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
     // the eigenvectors specify the transformation matrix
-    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
+    Matrix M = Matrix(evec->data);
+    ActOnAllVectors( &Vector::MatrixMultiplication, static_cast<const Matrix>(M));
     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
 
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/vector.cpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -8,10 +8,14 @@
 
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
 #include "Helpers/Assert.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "Exceptions/MathException.hpp"
 
 #include <iostream>
+#include <gsl/gsl_blas.h>
+
 
 using namespace std;
@@ -34,7 +38,5 @@
 {
   content = gsl_vector_alloc(NDIM);
-  gsl_vector_set(content,0,src[0]);
-  gsl_vector_set(content,1,src[1]);
-  gsl_vector_set(content,2,src[2]);
+  gsl_vector_memcpy(content, src.content);
 }
 
@@ -49,4 +51,8 @@
 };
 
+Vector::Vector(gsl_vector *_content) :
+  content(_content)
+{}
+
 /**
  * Assignment operator
@@ -55,7 +61,5 @@
   // check for self assignment
   if(&src!=this){
-    gsl_vector_set(content,0,src[0]);
-    gsl_vector_set(content,1,src[1]);
-    gsl_vector_set(content,2,src[2]);
+    gsl_vector_memcpy(content, src.content);
   }
   return *this;
@@ -101,16 +105,17 @@
 double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
 {
-  double res = distance(y), tmp, matrix[NDIM*NDIM];
+  double res = distance(y), tmp;
+  Matrix matrix;
     Vector Shiftedy, TranslationVector;
     int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
+    matrix.at(0,0) = cell_size[0];
+    matrix.at(1,0) = cell_size[1];
+    matrix.at(2,0) = cell_size[3];
+    matrix.at(0,1) = cell_size[1];
+    matrix.at(1,1) = cell_size[2];
+    matrix.at(2,1) = cell_size[4];
+    matrix.at(0,2) = cell_size[3];
+    matrix.at(1,2) = cell_size[4];
+    matrix.at(2,2) = cell_size[5];
     // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     for (N[0]=-1;N[0]<=1;N[0]++)
@@ -138,16 +143,17 @@
 double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
 {
-  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
+  double res = DistanceSquared(y), tmp;
+  Matrix matrix;
     Vector Shiftedy, TranslationVector;
     int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
+    matrix.at(0,0) = cell_size[0];
+    matrix.at(1,0) = cell_size[1];
+    matrix.at(2,0) = cell_size[3];
+    matrix.at(0,1) = cell_size[1];
+    matrix.at(1,1) = cell_size[2];
+    matrix.at(2,1) = cell_size[4];
+    matrix.at(0,2) = cell_size[3];
+    matrix.at(1,2) = cell_size[4];
+    matrix.at(2,2) = cell_size[5];
     // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     for (N[0]=-1;N[0]<=1;N[0]++)
@@ -172,6 +178,7 @@
  * Tries to translate a vector into each adjacent neighbouring cell.
  */
-void Vector::KeepPeriodic(const double * const matrix)
-{
+void Vector::KeepPeriodic(const double * const _matrix)
+{
+  Matrix matrix = Matrix(_matrix);
   //  int N[NDIM];
   //  bool flag = false;
@@ -181,5 +188,5 @@
   //  Output(out);
   //  Log() << Verbose(0) << endl;
-    InverseMatrixMultiplication(matrix);
+    MatrixMultiplication(matrix.invert());
     for(int i=NDIM;i--;) { // correct periodically
       if (at(i) < 0) {  // get every coefficient into the interval [0,1)
@@ -203,6 +210,5 @@
 {
   double res = 0.;
-  for (int i=NDIM;i--;)
-    res += at(i)*y[i];
+  gsl_blas_ddot(content, y.content, &res);
   return (res);
 };
@@ -461,4 +467,16 @@
 };
 
+Vector &Vector::operator*=(const Matrix &mat){
+  (*this) = mat*(*this);
+  return *this;
+}
+
+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, 0.0, res);
+  return Vector(res);
+}
+
+
 /** Factors given vector \a a times \a m.
  * \param a vector
@@ -508,6 +526,5 @@
 void Vector::Scale(const double factor)
 {
-  for (int i=NDIM;i--;)
-    at(i) *= factor;
+  gsl_vector_scale(content,factor);
 };
 
@@ -516,6 +533,8 @@
  * \param *Minv inverse matrix
  */
-void Vector::WrapPeriodically(const double * const M, const double * const Minv)
-{
+void Vector::WrapPeriodically(const double * const _M, const double * const _Minv)
+{
+  Matrix M = Matrix(_M);
+  Matrix Minv = Matrix(_Minv);
   MatrixMultiplication(Minv);
   // truncate to [0,1] for each axis
@@ -550,12 +569,7 @@
  * \param *matrix NDIM_NDIM array
  */
-void Vector::MatrixMultiplication(const double * const M)
-{
-  Vector tmp;
-  // do the matrix multiplication
-  for(int i=NDIM;i--;)
-    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
-
-  (*this) = tmp;
+void Vector::MatrixMultiplication(const Matrix &M)
+{
+  (*this) *= M;
 };
 
@@ -565,4 +579,5 @@
 bool Vector::InverseMatrixMultiplication(const double * const A)
 {
+  /*
   double B[NDIM*NDIM];
   double detA = RDET3(A);
@@ -586,4 +601,13 @@
     return true;
   } else {
+    return false;
+  }
+  */
+  Matrix mat = Matrix(A);
+  try{
+    (*this) *= mat.invert();
+    return true;
+  }
+  catch(MathException &excpt){
     return false;
   }
@@ -679,6 +703,5 @@
 void Vector::AddVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) += y[i];
+  gsl_vector_add(content, y.content);
 }
 
@@ -688,6 +711,5 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) -= y[i];
+  gsl_vector_sub(content, y.content);
 }
 
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision f0f1cc841801082b74fbe0a203cc8ebbbc78b481)
+++ src/vector.hpp	(revision c7b39aba0c5bbf0cd82e49fb3793c618e6374c7e)
@@ -24,4 +24,5 @@
 
 class Vector;
+class Matrix;
 
 typedef std::vector<Vector> pointset;
@@ -31,6 +32,6 @@
  */
 class Vector : public Space{
+  friend Vector operator*(const Matrix&,const Vector&);
 public:
-
   Vector();
   Vector(const double x1, const double x2, const double x3);
@@ -59,5 +60,5 @@
   void ScaleAll(const double *factor);
   void Scale(const double factor);
-  void MatrixMultiplication(const double * const M);
+  void MatrixMultiplication(const Matrix &M);
   bool InverseMatrixMultiplication(const double * const M);
   void KeepPeriodic(const double * const matrix);
@@ -96,4 +97,5 @@
   Vector const operator+(const Vector& b) const;
   Vector const operator-(const Vector& b) const;
+  Vector& operator*=(const Matrix&);
 
   // Methods inherited from Space
@@ -104,4 +106,5 @@
 
 private:
+  Vector(gsl_vector *);
   gsl_vector *content;
 
@@ -118,4 +121,6 @@
 Vector const operator*(const Vector& a, const double m);
 Vector const operator*(const double m, const Vector& a);
+Vector operator*(const Matrix&,const Vector&);
+
 
 #endif /*VECTOR_HPP_*/
