Index: src/Exceptions/NotInvertibleException.cpp
===================================================================
--- src/Exceptions/NotInvertibleException.cpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
+++ src/Exceptions/NotInvertibleException.cpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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 325390d0efa70ca974c1bc1c47244bd26c28bc39)
+++ src/Exceptions/NotInvertibleException.hpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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 923b6cca0028d3b4eea5fdb75a369201bd898f35)
+++ src/Makefile.am	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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
@@ -124,17 +145,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} \
@@ -146,5 +156,4 @@
   ${HELPERSOURCE} \
   ${LEGACYSOURCE} \
-  ${EXCEPTIONSOURCE} \
   bond.cpp \
   bondgraph.cpp \
@@ -192,5 +201,4 @@
   ${PATTERNHEADER} \
   ${DESCRIPTORHEADER} \
-  ${EXCEPTIONHEADER} \
   ${LEGACYHEADER} \
   bond.hpp \
Index: src/Matrix.cpp
===================================================================
--- src/Matrix.cpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
+++ src/Matrix.cpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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 325390d0efa70ca974c1bc1c47244bd26c28bc39)
+++ src/Matrix.hpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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/vector.cpp
===================================================================
--- src/vector.cpp	(revision 923b6cca0028d3b4eea5fdb75a369201bd898f35)
+++ src/vector.cpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -8,8 +8,10 @@
 
 #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>
@@ -48,4 +50,8 @@
   gsl_vector_set(content,2,x3);
 };
+
+Vector::Vector(gsl_vector *_content) :
+  content(_content)
+{}
 
 /**
@@ -458,4 +464,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
@@ -548,10 +566,5 @@
 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;
+  (*this) *= Matrix(M);
 };
 
@@ -561,4 +574,5 @@
 bool Vector::InverseMatrixMultiplication(const double * const A)
 {
+  /*
   double B[NDIM*NDIM];
   double detA = RDET3(A);
@@ -582,4 +596,13 @@
     return true;
   } else {
+    return false;
+  }
+  */
+  Matrix mat = Matrix(A);
+  try{
+    (*this) *= mat.invert();
+    return true;
+  }
+  catch(MathException &excpt){
     return false;
   }
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 923b6cca0028d3b4eea5fdb75a369201bd898f35)
+++ src/vector.hpp	(revision 325390d0efa70ca974c1bc1c47244bd26c28bc39)
@@ -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);
@@ -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_*/
