source: src/LinearAlgebra/Eigenspace.cpp@ 7d059d

/*
 * 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!");

  createEigenvectorViews();
}


/** 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 );
    Log() << Verbose(1) << iter << " th Eigenvector is " << *ev << std::endl;
  }
}


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


const VectorContent & Eigenspace::getEigenvector(size_t i) const
{
  return *(Eigenvectors[i]);
}

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


/** Sets the eigenvectors and -values.
 *
 * We only implement this for eigenpairs, not eigenvector and eigenvalue
 * distinctly to avoid mismatches, i.e. old eigenvalues used because only new
 * eigenvectors have been given.
 *
 * @param CurrentEigenvectors eigenvectors to set to
 * @param CurrentEigenvalues eigenvalues to set to
 */
void Eigenspace::setEigenpairs(const eigenvectorset &CurrentEigenvectors, const eigenvalueset &CurrentEigenvalues)
{
  ASSERT(CurrentEigenvectors.size() == Eigenvectors.size(),
      "Eigenspace::setEigenpairs() - "
      +toString(CurrentEigenvectors.size())+" x CurrentEigenvectors is not the same as "
      +toString(Eigenvectors.size())+"x Eigenvectors!");
  // set eigenvectors
  std::vector<boost::shared_ptr<VectorContent> >::const_iterator newev = CurrentEigenvectors.begin();
  for(std::vector<boost::shared_ptr<VectorContent> >::iterator ev = Eigenvectors.begin();
    ev != Eigenvectors.end();
    ++ev, ++newev) {
    *ev = *newev;
  }
  // eigenvalues can be just copied (integral type)
  Eigenvalues = CurrentEigenvalues;
}


/** Output operator for Eigenspace.
 *
 * @param ost output stream
 * @param _s Eigenspace to print: Indices and contained matrix
 * @return given output stream
 */
std::ostream & operator<<(std::ostream &ost, const Eigenspace &_s)
{
  BOOST_FOREACH(size_t iter, _s.Indices) {
    ost << iter << " ";
  }
  ost << " with matrix " << _s.EigenspaceMatrix;
  return ost;
}