Index: src/LinearAlgebra/Eigenspace.cpp
===================================================================
--- src/LinearAlgebra/Eigenspace.cpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
+++ src/LinearAlgebra/Eigenspace.cpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
@@ -0,0 +1,147 @@
+/*
+ * 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.
+ */
+
+/*
+ * Eigenspace.cpp
+ *
+ *  Created on: Nov 22, 2010
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Helpers/MemDebug.hpp"
+
+#include <iostream>
+#include <boost/foreach.hpp>
+
+#include "Helpers/Assert.hpp"
+#include "Helpers/Log.hpp"
+#include "Helpers/Verbose.hpp"
+
+#include "Eigenspace.hpp"
+
+
+/** Constructor for Eigenspace.
+ * @param _Indices index set for columns and rows of EigenspaceMatrix.
+ */
+Eigenspace::Eigenspace(const indexset &_Indices) :
+  Indices(_Indices),
+  EigenspaceMatrix(_Indices.size(),_Indices.size()),
+  EigenvectorMatrix(_Indices.size(),_Indices.size())
+{
+  createEigenvectorViews();
+}
+
+
+/** Constructor for Eigenspace.
+ * @param _Indices index set for columns and rows of EigenspaceMatrix.
+ * @param _matrix matrix to use
+ */
+Eigenspace::Eigenspace(const indexset &_Indices, const MatrixContent &_matrix) :
+  Indices(_Indices),
+  EigenspaceMatrix(_matrix),
+  EigenvectorMatrix(_matrix.getRows(), _matrix.getColumns())
+{
+  ASSERT(_Indices.size() == _matrix.getRows(),
+      "Eigenspace::Eigenspace() - number of indices "+toString(_Indices.size())
+      +" and row dimension of matrix "+toString(_matrix.getRows())+" don't match!");
+  ASSERT(_Indices.size() == _matrix.getColumns(),
+      "Eigenspace::Eigenspace() - number of indices "+toString(_Indices.size())
+      +" and column dimension of matrix "+toString(_matrix.getColumns())+" don't match!");
+}
+
+
+/** Destructor for Eigenspace.
+ *
+ */
+Eigenspace::~Eigenspace()
+{}
+
+
+/** Create a number of eigenvectors in Eigenspace::Eigenvectors.
+ * These are views on the resepctive column vectors in
+ * Eigenspace::EigenvectorMatrix.
+ */
+void Eigenspace::createEigenvectorViews()
+{
+  const size_t ColumnCount = EigenspaceMatrix.getColumns();
+  for (size_t iter = 0; iter < ColumnCount; ++iter) {
+    boost::shared_ptr<VectorContent> ev(EigenspaceMatrix.getColumnVector(iter));
+    Eigenvectors.push_back( ev );
+  }
+}
+
+
+/** Obtains the eigenvectors to Eigenspace::EigenspaceMatrix.
+ * We diagonalize the matrix and set the Eigenspace::Eigenvectors and
+ * Eigenspace::Eigenvalues.
+ */
+void Eigenspace::calculateEigenSubspace()
+{
+  EigenvectorMatrix = EigenspaceMatrix;
+  VectorContent *Eigenvalues = new VectorContent(EigenvectorMatrix.transformToEigenbasis());
+  Log() << Verbose(2) << "Eigenvector matrix is " << EigenvectorMatrix << std::endl;
+  Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
+}
+
+
+/** Getter for constant reference to Eigenspace::EigenspaceMatrix.
+ *
+ * @return Eigenspace::EigenspaceMatrix
+ */
+const MatrixContent & Eigenspace::getEigenspaceMatrix() const
+{
+  return EigenspaceMatrix;
+}
+
+
+/** Setter for Eigenspace::EigenspaceMatrix.
+ *
+ * @param _content matrix to use
+ */
+void Eigenspace::setEigenspaceMatrix(const MatrixContent &_content)
+{
+  EigenspaceMatrix = _content;
+  EigenvectorMatrix.setZero();
+}
+
+/** Checks whether each index of _indexset is contained in this subspace:Indices.
+ *
+ * @param _subspace subspace set to check
+ * @return true - _index is contained, false - at least index is not contained
+ */
+bool Eigenspace::contains(const Eigenspace &_subspace)
+{
+  BOOST_FOREACH(size_t iter, _subspace.getIndices()) {
+    if (Indices.count(iter) == 0) {  // index not present, then not contained
+      return false;
+    }
+  }
+  return true;
+}
+
+
+/** Getter for Eigenspace::Indices.
+ *
+ * @return indexset of this eigenspace
+ */
+const Eigenspace::indexset & Eigenspace::getIndices() const
+{
+  return Indices;
+}
+
+std::ostream & operator<<(std::ostream &ost, const Eigenspace &_s)
+{
+  BOOST_FOREACH(size_t iter, _s.getIndices()) {
+    ost << iter << " ";
+  }
+  return ost;
+}
Index: src/LinearAlgebra/Eigenspace.hpp
===================================================================
--- src/LinearAlgebra/Eigenspace.hpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
+++ src/LinearAlgebra/Eigenspace.hpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
@@ -0,0 +1,63 @@
+/*
+ * Eigenspace.hpp
+ *
+ *  Created on: Nov 22, 2010
+ *      Author: heber
+ */
+
+#ifndef EIGENSPACE_HPP_
+#define EIGENSPACE_HPP_
+
+#include <boost/shared_ptr.hpp>
+
+#include <iosfwd>
+#include <set>
+#include <vector>
+#include "MatrixContent.hpp"
+#include "VectorContent.hpp"
+
+/** Container for a matrix and its eigenvectors and -values.
+ *
+ *  Note that a set of eigenvectors always forms a subspace within a vector
+ *  space, the so-called eigenspace. Hence, the name of this class is
+ *  eigenspace.
+ *
+ *  Eigenvector are stored internally as a matrix. A returned set of
+ *  eigenvectors is just a set of VectorViewContent on row vectors in
+ *  this matrix.
+ */
+class Eigenspace
+{
+public:
+  typedef std::vector< boost::shared_ptr<VectorContent> > eigenvectorset;
+  typedef std::vector< double > eigenvalueset;
+  typedef std::set<size_t> indexset;
+
+  Eigenspace(const indexset &_Indices);
+  Eigenspace(const indexset &_Indices, const MatrixContent &_matrix);
+  virtual ~Eigenspace();
+
+  void calculateEigenSubspace();
+
+  // accessors
+  const MatrixContent & getEigenspaceMatrix() const;
+  const indexset & getIndices() const;
+  void setEigenspaceMatrix(const MatrixContent &_content);
+
+  bool contains(const Eigenspace &_indexset);
+
+protected:
+  const indexset Indices;
+  eigenvectorset Eigenvectors;
+  eigenvalueset Eigenvalues;
+  MatrixContent EigenspaceMatrix;
+  MatrixContent EigenvectorMatrix;
+
+private:
+  void createEigenvectorViews();
+
+};
+
+std::ostream & operator<<(std::ostream &ost, const Eigenspace &_s);
+
+#endif /* EIGENSPACE_HPP_ */
Index: src/LinearAlgebra/Makefile.am
===================================================================
--- src/LinearAlgebra/Makefile.am	(revision 74237190563beadc942e60dbd56f5fd7b388031b)
+++ src/LinearAlgebra/Makefile.am	(revision 9c529645733d09290760e9ebeff51081f40223b8)
@@ -12,4 +12,5 @@
   ${HELPERSOURCE} \
   BoxVector.cpp \
+  Eigenspace.cpp \
   leastsquaremin.cpp \
   Line.cpp \
@@ -22,4 +23,5 @@
   RealSpaceMatrix.cpp \
   Space.cpp \
+  Subspace.cpp \
   vector_ops.cpp \
   Vector.cpp \
@@ -29,4 +31,5 @@
 LINALGHEADER = \
   BoxVector.hpp \
+  Eigenspace.hpp \
   leastsquaremin.hpp \
   Line.hpp \
@@ -39,4 +42,5 @@
   Plane.hpp \
   Space.hpp \
+  Subspace.hpp \
   vector_ops.hpp \
   Vector.hpp \
Index: src/LinearAlgebra/Subspace.cpp
===================================================================
--- src/LinearAlgebra/Subspace.cpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
+++ src/LinearAlgebra/Subspace.cpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
@@ -0,0 +1,213 @@
+/*
+ * 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.
+ */
+
+/*
+ * Subspace.cpp
+ *
+ *  Created on: Nov 22, 2010
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Helpers/MemDebug.hpp"
+
+#include <boost/shared_ptr.hpp>
+#include <boost/foreach.hpp>
+
+#include "Helpers/Assert.hpp"
+#include "Helpers/Log.hpp"
+#include "Helpers/Verbose.hpp"
+
+#include "Eigenspace.hpp"
+#include "Subspace.hpp"
+
+/** Constructor for class Subspace.
+ *
+ * @param _s indices that make up this subspace
+ * @param _FullSpace reference to the full eigenspace
+ */
+Subspace::Subspace(indexset & _s, Eigenspace &_FullSpace) :
+  Eigenspace(_s),
+  ProjectToSubspace(_s.size(), _FullSpace.getIndices().size()),
+  ProjectFromSubspace(_FullSpace.getIndices().size(), _s.size()),
+  FullSpace(_FullSpace)
+{}
+
+
+/** Destructor for class Subspace.
+ *
+ */
+Subspace::~Subspace()
+{}
+
+
+/** Diagonalizes the subspace matrix.
+ *
+ * @param fullmatrix full matrix to project into subspace and solve
+ */
+void Subspace::calculateEigenSubspace()
+{
+  getSubspacematrixFromBigmatrix(FullSpace.getEigenspaceMatrix());
+}
+
+
+/** Projects a given full matrix down into the subspace of this class.
+ *
+ * @param bigmatrix full matrix to project into subspace
+ */
+void Subspace::getSubspacematrixFromBigmatrix(const MatrixContent & bigmatrix)
+{
+  // check whether subsystem is big enough for both index sets
+  ASSERT(Indices.size() <= bigmatrix.getRows(),
+      "embedEigenspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
+      +" than needed by index set "
+      +toString(Indices.size())+"!");
+  ASSERT(Indices.size() <= bigmatrix.getColumns(),
+      "embedEigenspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
+      +" than needed by index set "
+      +toString(Indices.size())+"!");
+
+  // construct small matrix
+  MatrixContent *subsystem =  new MatrixContent(Indices.size(), Indices.size());
+  size_t localrow = 0; // local row indices for the subsystem
+  size_t localcolumn = 0;
+  BOOST_FOREACH( size_t rowindex, Indices) {
+    localcolumn = 0;
+    BOOST_FOREACH( size_t columnindex, Indices) {
+      ASSERT((rowindex < bigmatrix.getRows()) && (columnindex < bigmatrix.getColumns()),
+          "current index pair ("
+          +toString(rowindex)+","+toString(columnindex)
+          +") is out of bounds of bigmatrix ("
+          +toString(bigmatrix.getRows())+","+toString(bigmatrix.getColumns())
+          +")");
+      subsystem->at(localrow,localcolumn) = (*Eigenvectors[rowindex]) * (bigmatrix * (*Eigenvectors[columnindex]));
+      localcolumn++;
+    }
+    localrow++;
+  }
+}
+
+
+void Subspace::correctEigenvectorsFromSubIndices()
+{
+}
+
+
+/** Creates the inverse of LocalToGlobal.
+ * Mapping is one-to-one and onto, hence it's simply turning around the map.
+ */
+void Subspace::invertLocalToGlobalMapping()
+{
+  GlobalToLocal.clear();
+  for (mapping::const_iterator iter = LocalToGlobal.begin();
+      iter != LocalToGlobal.end();
+      ++iter) {
+    GlobalToLocal.insert( make_pair(iter->second, iter->first) );
+  }
+}
+
+
+/** Project calculated Eigenvectors into full space.
+ *
+ * @return set of projected eigenvectors.
+ */
+Eigenspace::eigenvectorset Subspace::getEigenvectorsInFullSpace()
+{
+  // check whether bigmatrix is at least as big as EigenspaceMatrix
+  ASSERT(Eigenvectors.size() > 0,
+      "embedEigenspaceMatrix() - no Eigenvectors given!");
+  ASSERT(EigenspaceMatrix.getRows() <= Eigenvectors[0]->getDimension(),
+      "embedEigenspaceMatrix() - EigenspaceMatrix has more rows "
+      +toString(EigenspaceMatrix.getRows())+" than eigenvectors "
+      +toString(Eigenvectors[0]->getDimension())+"!");
+  ASSERT(EigenspaceMatrix.getColumns() <= Eigenvectors.size(),
+      "embedEigenspaceMatrix() - EigenspaceMatrix has more columns "
+      +toString(EigenspaceMatrix.getColumns())+" than eigenvectors "
+      +toString(Eigenvectors.size())+"!");
+  // check whether EigenspaceMatrix is big enough for both index sets
+  ASSERT(EigenspaceMatrix.getColumns() == EigenspaceMatrix.getRows(),
+      "embedEigenspaceMatrix() - EigenspaceMatrix is not square "
+      +toString(EigenspaceMatrix.getRows())+" than needed by index set "
+      +toString(EigenspaceMatrix.getColumns())+"!");
+  ASSERT(Indices.size() == EigenspaceMatrix.getColumns(),
+      "embedEigenspaceMatrix() - EigenspaceMatrix has not the same number of columns "
+      +toString(EigenspaceMatrix.getColumns())+" compared to the index set "
+      +toString(Indices.size())+"!");
+
+  // construct intermediate matrix
+  MatrixContent *intermediatematrix = new MatrixContent(Eigenvectors[0]->getDimension(), Indices.size());
+  size_t localcolumn = 0;
+  BOOST_FOREACH(size_t columnindex, Indices) {
+    // create two vectors from each row and copy assign them
+    boost::shared_ptr<VectorContent> srceigenvector(Eigenvectors[columnindex]);
+    boost::shared_ptr<VectorContent> desteigenvector(intermediatematrix->getColumnVector(localcolumn));
+    *desteigenvector = *srceigenvector;
+    localcolumn++;
+  }
+  // matrix product with eigenbasis EigenspaceMatrix matrix
+  *intermediatematrix *= EigenspaceMatrix;
+
+  // and place at right columns into bigmatrix
+  MatrixContent *bigmatrix = new MatrixContent(Eigenvectors[0]->getDimension(), Eigenvectors.size());
+  bigmatrix->setZero();
+  localcolumn = 0;
+  BOOST_FOREACH(size_t columnindex, Indices) {
+    // create two vectors from each row and copy assign them
+    boost::shared_ptr<VectorContent> srceigenvector(intermediatematrix->getColumnVector(localcolumn));
+    boost::shared_ptr<VectorContent> desteigenvector(bigmatrix->getColumnVector(columnindex));
+    *desteigenvector = *srceigenvector;
+    localcolumn++;
+  }
+
+  return Eigenvectors;
+}
+
+
+/** Remove a subset of indices from the SubIndices.
+ *
+ * @param _s subset to remove
+ * @return true - removed, false - not found
+ */
+bool Subspace::removeSubset(boost::shared_ptr<Subspace> &_s)
+{
+  if (SubIndices.count(_s) != 0) {
+    SubIndices.erase(_s);
+    return true;
+  } else {
+    return false;
+  }
+}
+
+/** Add a subspace set to SubIndices.
+ *
+ * @param _s subset to add
+ * @return true - not present before, false - has already been present
+ */
+bool Subspace::addSubset(boost::shared_ptr<Subspace> & _s)
+{
+  if (SubIndices.count(_s) != 0)
+    return false;
+  else {
+    SubIndices.insert(_s);
+    return true;
+  }
+}
+
+
+/** Simply a getter for Eigenvectors, as they are stored as subspace eigenvectors.
+ *
+ * @return set of eigenvectors in subspace
+ */
+Eigenspace::eigenvectorset Subspace::getEigenvectorsInSubspace()
+{
+  return Eigenvectors;
+}
+
Index: src/LinearAlgebra/Subspace.hpp
===================================================================
--- src/LinearAlgebra/Subspace.hpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
+++ src/LinearAlgebra/Subspace.hpp	(revision 9c529645733d09290760e9ebeff51081f40223b8)
@@ -0,0 +1,63 @@
+/*
+ * Subspace.hpp
+ *
+ *  Created on: Nov 22, 2010
+ *      Author: heber
+ */
+
+#ifndef SUBSPACE_HPP_
+#define SUBSPACE_HPP_
+
+#include <map>
+#include <set>
+#include <vector>
+#include "Eigenspace.hpp"
+#include "MatrixContent.hpp"
+#include "VectorContent.hpp"
+
+/** A subset of eigenvectors from an Eigenspace.
+ *
+ * In this class we regard a sub set of eigenvectors of an Eigenspace
+ * which span a subspace of the eigenspace. This is used for diagonalization
+ * of the Eigenspace's matrix in linear-scaling, subspace decomposition
+ * schemes.
+ *
+ * Here, beyond the contents of Eigenspace, we need projection matrices from
+ * and to this subspace and also mappings from the global indices to the local
+ * indices, to identify local eigenvectors in this Subspace with their
+ * counterparts in the full Eigenspace.
+ *
+ */
+class Subspace : public Eigenspace
+{
+public:
+  typedef std::map<size_t, size_t> mapping;
+  typedef std::set< boost::shared_ptr<Subspace> > subset;
+
+  Subspace(indexset &_s, Eigenspace &_FullSpace);
+  ~Subspace();
+
+  // manipulate subsets
+  bool addSubset(boost::shared_ptr<Subspace> &_s);
+  bool removeSubset(boost::shared_ptr<Subspace> &_s);
+
+  void calculateEigenSubspace();
+  void correctEigenvectorsFromSubIndices();
+  eigenvectorset getEigenvectorsInFullSpace();
+  eigenvectorset getEigenvectorsInSubspace();
+
+private:
+
+  void invertLocalToGlobalMapping();
+  void getSubspacematrixFromBigmatrix(const MatrixContent & bigmatrix);
+
+  mapping LocalToGlobal;
+  mapping GlobalToLocal;
+  subset SubIndices;
+  MatrixContent ProjectToSubspace;
+  MatrixContent ProjectFromSubspace;
+  Eigenspace &FullSpace;
+};
+
+
+#endif /* SUBSPACE_HPP_ */
