source: src/LinearAlgebra/Eigenspace.cpp@ 9c5296

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