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;
+}
+
