source: src/LinearAlgebra/RealSpaceMatrix.cpp@ b88fe4

/*
 * 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 "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Assert.hpp"
#include "LinearAlgebra/defs.hpp"
#include "LinearAlgebra/Exceptions.hpp"
#include "LinearAlgebra/fast_functions.hpp"
#include "LinearAlgebra/MatrixContent.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.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>
#include <limits>

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 = content->getRowVector(i);
    rows_ptr[i] = new Vector(rowContent);
    ASSERT(rowContent == NULL, "RealSpaceMatrix::createViews() - rowContent was not taken over as supposed to happen!");
  }
  // create column views
  for(int i=NDIM;i--;){
    VectorContent *columnContent = content->getColumnVector(i);
    columns_ptr[i] = new Vector(columnContent);
    ASSERT(columnContent == NULL, "RealSpaceMatrix::createViews() - columnContent was not taken over as supposed to happen!");
  }
  // create diagonal view
  VectorContent *diagonalContent = content->getDiagonalVector();
  diagonal_ptr = new Vector(diagonalContent);
  ASSERT(diagonalContent == NULL, "RealSpaceMatrix::createViews() - diagonalContent was not taken over as supposed to happen!");
}

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));
//  double phi[NDIM];
//  phi[0] = x;
//  phi[0] = y;
//  phi[0] = z;
//  setRotation(phi);
}

void RealSpaceMatrix::setRotation(const double phi[NDIM])
{
  setRotation(phi[0], phi[1], phi[2]);
//  set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
//  set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
//  set(0,2,              cos(phi[1])*sin(phi[2])                                        );
//  set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
//  set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
//  set(1,2,              sin(phi[1])                                                    );
//  set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
//  set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
//  set(2,2,              cos(phi[1])*cos(phi[2])                                        );
}

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) <= LINALG_MYEPSILON()){
    throw NotInvertibleException() << LinearAlgebraDeterminant(det) << LinearAlgebraMatrixContent(this->content);
  }

  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
{
  RealSpaceMatrix res = RealSpaceMatrix(const_cast<const MatrixContent *>(content)->transpose());
  return res;
}

RealSpaceMatrix &RealSpaceMatrix::transpose()
{
  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;
}