source: src/LinearAlgebra/Subspace.cpp@ a06042

/*
 * 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(_FullSpace.getIndices().size(), _s.size()),
  ProjectFromSubspace(_s.size(), _FullSpace.getIndices().size()),
  FullSpace(_FullSpace)
{
  // create projection matrices
  createProjectionMatrices();

  // create eigenspace matrices
  projectFullSpaceMatrixToSubspace();
}


/** 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;
}


/** Creates the projection matrix from Subspace::FullSpace::EigenvectorMatrix.
 * We simply copy the respective eigenvectors from Subspace::FullSpace into
 * Subspace::Eigenvectors.
 */
void Subspace::createProjectionMatrices()
{
  // first ProjectToSubspace

  // generate column vectors
  std::vector< boost::shared_ptr<VectorContent> > ColumnVectors;
  for (size_t i = 0; i < Indices.size(); ++i) {
    boost::shared_ptr<VectorContent> p(ProjectToSubspace.getColumnVector(i));
    ColumnVectors.push_back( p );
  }

  // then copy from real eigenvectors
  size_t localindex = 0;
  BOOST_FOREACH(size_t iter, Indices) {
    *ColumnVectors[localindex] = FullSpace.getEigenvector(iter);
    localindex++;
  }
  Log() << Verbose(2) << "ProjectionToSubspace matrix is " << ProjectToSubspace << std::endl;

  // then ProjectFromSubspace is just transposed
  ProjectFromSubspace = ((const MatrixContent)ProjectToSubspace).transpose();
  Log() << Verbose(2) << "ProjectFromSubspace matrix is " << ProjectFromSubspace << std::endl;
}


/** Creates the subspace matrix by projecting down the FullSpace::EigenspaceMatrix.
 *  We just perform \f$M = P^t \cdot A \cdot P\f$, when P are the projection matrix,
 *  A the full matrix and M the subspace matrix.
 */
void Subspace::projectFullSpaceMatrixToSubspace()
{
  // construct small matrix
  EigenspaceMatrix = ProjectFromSubspace * FullSpace.getEigenspaceMatrix() * ProjectToSubspace;
  Log() << Verbose(2) << "EigenspaceMatrix matrix is " << EigenspaceMatrix << std::endl;
}