source: src/LinearAlgebra/Matrix.cpp@ 3bc926

/*
 * 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.
 */

/*
 * Matrix.cpp
 *
 *  Created on: Jun 25, 2010
 *      Author: crueger
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Helpers/MemDebug.hpp"

#include "LinearAlgebra/Matrix.hpp"
#include "Helpers/Assert.hpp"
#include "Exceptions/NotInvertibleException.hpp"
#include "Helpers/fast_functions.hpp"
#include "Helpers/Assert.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "LinearAlgebra/MatrixContent.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>

using namespace std;

Matrix::Matrix()
{
  content = new MatrixContent(NDIM, NDIM);
  createViews();
}

Matrix::Matrix(const double* src)
{
  content = new MatrixContent(NDIM, NDIM, src);
  createViews();
}

Matrix::Matrix(const Matrix &src)
{
  content = new MatrixContent(src.content);
  createViews();
}

Matrix::Matrix(const MatrixContent &src)
{
  content = new MatrixContent(src);
  createViews();
}

Matrix::Matrix(MatrixContent* _content)
{
  content = new MatrixContent(_content);
  createViews();
}

Matrix::~Matrix()
{
  // 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 Matrix::createViews(){
  // create row views
  for(int i=NDIM;i--;){
    VectorContent *rowContent = new VectorViewContent(gsl_matrix_row(content->content,i));
    rows_ptr[i] = new Vector(rowContent);
  }
  // create column views
  for(int i=NDIM;i--;){
    VectorContent *columnContent = new VectorViewContent(gsl_matrix_column(content->content,i));
    columns_ptr[i] = new Vector(columnContent);
  }
  // create diagonal view
  VectorContent *diagonalContent = new VectorViewContent(gsl_matrix_diagonal(content->content));
  diagonal_ptr = new Vector(diagonalContent);
}

void Matrix::setIdentity(){
  content->setIdentity();
}

void Matrix::setZero(){
  content->setZero();
}

void Matrix::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));
}

Matrix &Matrix::operator=(const Matrix &src)
{
  // MatrixContent checks for self-assignment
  *content = *(src.content);
  return *this;
}

const Matrix &Matrix::operator+=(const Matrix &rhs)
{
  *content += *(rhs.content);
  return *this;
}

const Matrix &Matrix::operator-=(const Matrix &rhs)
{
  *content -= *(rhs.content);
  return *this;
}

const Matrix &Matrix::operator*=(const Matrix &rhs)
{
  (*content) *= (*rhs.content);
  return *this;
}

const Matrix Matrix::operator+(const Matrix &rhs) const
{
  Matrix tmp = *this;
  tmp+=rhs;
  return tmp;
}

const Matrix Matrix::operator-(const Matrix &rhs) const
{
  Matrix tmp = *this;
  tmp-=rhs;
  return tmp;
}

const Matrix Matrix::operator*(const Matrix &rhs) const
{
  Matrix tmp(content);
  tmp *= rhs;
  return tmp;
}

double &Matrix::at(size_t i, size_t j)
{
  return content->at(i,j);
}

const double Matrix::at(size_t i, size_t j) const
{
  return content->at(i,j);
}

Vector &Matrix::row(size_t i)
{
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  return *rows_ptr[i];
}

const Vector &Matrix::row(size_t i) const
{
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  return *rows_ptr[i];
}

Vector &Matrix::column(size_t i)
{
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  return *columns_ptr[i];
}

const Vector &Matrix::column(size_t i) const
{
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  return *columns_ptr[i];
}

Vector &Matrix::diagonal()
{
  return *diagonal_ptr;
}

const Vector &Matrix::diagonal() const
{
  return *diagonal_ptr;
}

void Matrix::set(size_t i, size_t j, const double value)
{
  content->set(i,j, value);
}

double Matrix::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);
}


Matrix Matrix::invert() const{
  double det = determinant();
  if(fabs(det)<MYEPSILON){
    throw NotInvertibleException(__FILE__,__LINE__);
  }

  double detReci = 1./det;
  Matrix 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;
}

Matrix Matrix::transpose() const
{
  std::cout << "const Matrix::transpose()." << std::endl;
  Matrix res = Matrix(const_cast<const MatrixContent *>(content)->transpose());
  return res;
}

Matrix &Matrix::transpose()
{
  std::cout << "Matrix::transpose()." << std::endl;
  content->transpose();
  return *this;
}

Vector Matrix::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 Matrix &Matrix::operator*=(const double factor)
    {
  *content *= factor;
  return *this;
}

const Matrix operator*(const double factor,const Matrix& mat)
{
  Matrix tmp = mat;
  return tmp *= factor;
}

const Matrix operator*(const Matrix &mat,const double factor)
{
  return factor*mat;
}

bool Matrix::operator==(const Matrix &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
 */
Matrix ReturnFullMatrixforSymmetric(const double * const symm)
{
  Matrix 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 Matrix &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 Matrix &mat,const Vector &vec)
{
  return (*mat.content) * vec;
}

Vector &operator*=(Vector& lhs,const Matrix &mat)
{
  lhs = mat*lhs;
  return lhs;
}