source: src/LinearAlgebra/Matrix.cpp@ bf3817

/*
 * 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 "VectorContent.hpp"
#include "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();
  content->content = gsl_matrix_calloc(NDIM, NDIM);
  createViews();
}

Matrix::Matrix(const double* src){
  content = new MatrixContent();
  content->content = gsl_matrix_alloc(NDIM, NDIM);
  set(0,0, src[0]);
  set(1,0, src[1]);
  set(2,0, src[2]);

  set(0,1, src[3]);
  set(1,1, src[4]);
  set(2,1, src[5]);

  set(0,2, src[6]);
  set(1,2, src[7]);
  set(2,2, src[8]);
  createViews();
}

Matrix::Matrix(const Matrix &src){
  content = new MatrixContent();
  content->content = gsl_matrix_alloc(NDIM, NDIM);
  gsl_matrix_memcpy(content->content,src.content->content);
  createViews();
}

Matrix::Matrix(MatrixContent* _content) :
  content(_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;
  gsl_matrix_free(content->content);
  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::one(){
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      set(i,j,i==j);
    }
  }
}

void Matrix::zero(){
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      set(i,j,0.);
    }
  }
}

Matrix &Matrix::operator=(const Matrix &src){
  if(&src!=this){
    gsl_matrix_memcpy(content->content,src.content->content);
  }
  return *this;
}

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

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

const Matrix &Matrix::operator*=(const Matrix &rhs){
  (*this) = (*this)*rhs;
  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{
  gsl_matrix *res = gsl_matrix_alloc(NDIM, NDIM);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content->content, rhs.content->content, 0.0, res);
  MatrixContent *content= new MatrixContent();
  content->content = res;
  return Matrix(content);
}

double &Matrix::at(size_t i, size_t j){
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
  return *gsl_matrix_ptr (content->content, i, j);
}

const double Matrix::at(size_t i, size_t j) const{
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
  return gsl_matrix_get(content->content, 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){
  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
  gsl_matrix_set(content->content,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::transpose() const
{
  Matrix copy(*this);
  copy.transpose();
  return copy;
}


void Matrix::transpose()
{
  double tmp;
  for (int i=0;i<NDIM;i++)
    for (int j=i+1;j<NDIM;j++) {
      tmp = at(j,i);
      at(i,j) = tmp;
      at(j,i) = tmp;
    }
}


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

Vector Matrix::transformToEigenbasis()
{
  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
  gsl_vector *eval = gsl_vector_alloc(NDIM);
  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
  gsl_eigen_symmv(content->content, eval, evec, T);
  gsl_eigen_symmv_free(T);
  gsl_matrix_memcpy(content->content, evec);
  Vector evalues(gsl_vector_get(eval,0), gsl_vector_get(eval,1), gsl_vector_get(eval,2));
  return evalues;
}

const Matrix &Matrix::operator*=(const double factor){
  gsl_matrix_scale(content->content, factor);
  return *this;
}

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

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

bool Matrix::operator==(const Matrix &rhs) const{
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
        return false;
      }
    }
  }
  return true;
}

/** 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){
  gsl_vector *res = gsl_vector_calloc(NDIM);
  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content->content, vec.content->content, 0.0, res);
  VectorContent *content = new VectorContent();
  content->content = res;
  return Vector(content);
}

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