source: src/LinearAlgebra/MatrixContent.cpp@ ef0c6b

/*
 * MatrixContent.cpp
 *
 *  Created on: Nov 14, 2010
 *      Author: heber
 */


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

#include "Helpers/MemDebug.hpp"

#include "LinearAlgebra/RealSpaceMatrix.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_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <cmath>
#include <cassert>
#include <iostream>
#include <set>

using namespace std;


/** Constructor for class MatrixContent.
 * \param rows number of rows
 * \param columns number of columns
 */
MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
    rows(_rows),
    columns(_columns)
{
  content = gsl_matrix_calloc(rows, columns);
}

/** Constructor of class VectorContent.
 * We need this MatrixBaseCase for the VectorContentView class.
 * There no content should be allocated, as it is just a view with an internal
 * gsl_vector_view. Hence, MatrixBaseCase is just dummy class to give the
 * constructor a unique signature.
 * \param MatrixBaseCase
 */
MatrixContent::MatrixContent(size_t _rows, size_t _columns, MatrixBaseCase) :
    rows(_rows),
    columns(_columns)
{}

/** Constructor for class MatrixContent.
 * \param rows number of rows
 * \param columns number of columns
 * \param *src array with components to initialize matrix with
 */
MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
    rows(_rows),
    columns(_columns)
{
  content = gsl_matrix_calloc(rows, columns);
  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]);
}

/** Constructor for class MatrixContent.
 * We embed the given gls_matrix pointer within this class and set it to NULL
 * afterwards.
 * \param *src source gsl_matrix vector to embed within this class
 */
MatrixContent::MatrixContent(gsl_matrix *&src) :
    rows(src->size1),
    columns(src->size2)
{
  content = gsl_matrix_alloc(src->size1, src->size2);
  gsl_matrix_memcpy(content,src);
//  content = src;
//  src = NULL;
}

/** Copy constructor for class MatrixContent.
 * \param &src reference to source MatrixContent
 */
MatrixContent::MatrixContent(const MatrixContent &src) :
    rows(src.rows),
    columns(src.columns)
{
  content = gsl_matrix_alloc(src.rows, src.columns);
  gsl_matrix_memcpy(content,src.content);
}

/** Copy constructor for class MatrixContent.
 * \param *src pointer to source MatrixContent
 */
MatrixContent::MatrixContent(const MatrixContent *src) :
    rows(src->rows),
    columns(src->columns)
{
  ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
  content = gsl_matrix_alloc(src->rows, src->columns);
  gsl_matrix_memcpy(content,src->content);
}

/** Destructor for class MatrixContent.
 */
MatrixContent::~MatrixContent()
{
  gsl_matrix_free(content);
}

/** Set matrix to identity.
 */
void MatrixContent::setIdentity()
{
  for(int i=rows;i--;){
    for(int j=columns;j--;){
      set(i,j,i==j);
    }
  }
}

/** Set all matrix components to zero.
 */
void MatrixContent::setZero()
{
  for(int i=rows;i--;){
    for(int j=columns;j--;){
      set(i,j,0.);
    }
  }
}

/** Set all matrix components to a given value.
 * \param _value value to set each component to
 */
void MatrixContent::setValue(double _value)
{
  for(int i=rows;i--;){
    for(int j=columns;j--;){
      set(i,j,_value);
    }
  }
}

/** Copy operator for MatrixContent with self-assignment check.
 * \param &src matrix to compare to
 * \return reference to this
 */
MatrixContent &MatrixContent::operator=(const MatrixContent &src)
{
  if(&src!=this){
    gsl_matrix_memcpy(content,src.content);
  }
  return *this;
}

/** Addition operator.
 * \param &rhs matrix to add
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
{
  gsl_matrix_add(content, rhs.content);
  return *this;
}

/** Subtraction operator.
 * \param &rhs matrix to subtract
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
    {
  gsl_matrix_sub(content, rhs.content);
  return *this;
}

/** Multiplication operator.
 * Note that here matrix have to have same dimensions.
 * \param &rhs matrix to multiply with
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
{
  ASSERT(rows == rhs.rows,
      "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
  ASSERT(columns == rhs.columns,
      "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
  (*this) = (*this)*rhs;
  return *this;
}

/** Multiplication with copy operator.
 * \param &rhs matrix to multiply with
 * \return reference to newly allocated MatrixContent
 */
const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
{
  gsl_matrix *res = gsl_matrix_alloc(rows, rhs.columns);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
  // gsl_matrix is taken over by constructor, hence no free
  MatrixContent tmp(res);
  gsl_matrix_free(res);
  return tmp;
}

/* ========================== Accessing =============================== */

/** Accessor for manipulating component (i,j).
 * \param i row number
 * \param j column number
 * \return reference to component (i,j)
 */
double &MatrixContent::at(size_t i, size_t j)
{
  ASSERT((i>=0) && (i<rows),
      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
  ASSERT((j>=0) && (j<columns),
      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
  return *gsl_matrix_ptr (content, i, j);
}

/** Constant accessor for (value of) component (i,j).
 * \param i row number
 * \param j column number
 * \return const component (i,j)
 */
const double MatrixContent::at(size_t i, size_t j) const
{
  ASSERT((i>=0) && (i<rows),
      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
  ASSERT((j>=0) && (j<columns),
      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
  return gsl_matrix_get(content, i, j);
}

/** These functions return a pointer to the \a m-th element of a matrix.
 *  If \a m or \a n lies outside the allowed range of 0 to MatrixContent::dimension-1 then the error handler is invoked and a null pointer is returned.
 * \param m index
 * \return pointer to \a m-th element
 */
double *MatrixContent::Pointer(size_t m, size_t n)
{
  return gsl_matrix_ptr (content, m, n);
};

/** These functions return a constant pointer to the \a m-th element of a matrix.
 *  If \a m or \a n lies outside the allowed range of 0 to MatrixContent::dimension-1 then the error handler is invoked and a null pointer is returned.
 * \param m index
 * \return const pointer to \a m-th element
 */
const double *MatrixContent::const_Pointer(size_t m, size_t n) const
{
  return gsl_matrix_const_ptr (content, m, n);
};

/* ========================== Initializing =============================== */

/** Setter for component (i,j).
 * \param i row numbr
 * \param j column numnber
 * \param value value to set componnt (i,j) to
 */
void MatrixContent::set(size_t i, size_t j, const double value)
{
  ASSERT((i>=0) && (i<rows),
      "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
  ASSERT((j>=0) && (j<columns),
      "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
  gsl_matrix_set(content,i,j,value);
}

/** This function sets the matrix from a double array.
 * Creates a matrix view of the array and performs a memcopy.
 * \param *x array of values (no dimension check is performed)
 */
void MatrixContent::setFromDoubleArray(double * x)
{
  gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
  gsl_matrix_memcpy (content, &m.matrix);
};

/* ====================== Exchanging elements ============================ */
/** This function exchanges the \a i-th and \a j-th row of the matrix in-place.
 * \param i i-th row to swap with ...
 * \param j ... j-th row to swap against
 */
bool MatrixContent::SwapRows(size_t i, size_t j)
{
  return (gsl_matrix_swap_rows (content, i, j) == GSL_SUCCESS);
};

/** This function exchanges the \a i-th and \a j-th column of the matrix in-place.
 * \param i i-th column to swap with ...
 * \param j ... j-th column to swap against
 */
bool MatrixContent::SwapColumns(size_t i, size_t j)
{
  return (gsl_matrix_swap_columns (content, i, j) == GSL_SUCCESS);
};

/** This function exchanges the \a i-th row and \a j-th column of the matrix in-place.
 * The matrix must be square for this operation to be possible.
 * \param i i-th row to swap with ...
 * \param j ... j-th column to swap against
 */
bool MatrixContent::SwapRowColumn(size_t i, size_t j)
{
  assert (rows == columns && "The matrix must be square for swapping row against column to be possible.");
  return (gsl_matrix_swap_rowcol (content, i, j) == GSL_SUCCESS);
};

/** Return transposed matrix.
 * \return new matrix that is transposed of this.
 */
MatrixContent MatrixContent::transpose() const
{
  gsl_matrix *res = gsl_matrix_alloc(columns, rows);  // column and row dimensions exchanged!
  gsl_matrix_transpose_memcpy(res, content);
  MatrixContent newContent(res);
  gsl_matrix_free(res);
  return newContent;
}

/** Turn this matrix into its transposed.
 * Note that this is only possible if rows == columns.
 */
MatrixContent &MatrixContent::transpose()
{
  ASSERT( rows == columns,
      "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
  double tmp;
  for (size_t i=0;i<rows;i++)
    for (size_t j=i+1;j<rows;j++) {
      tmp = at(j,i);
      at(j,i) = at(i,j);
      at(i,j) = tmp;
    }
  return *this;
}

/** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
 * Note that we only return real-space part in case of non-symmetric matrix.
 * \warn return vector has to be freed'd
 * TODO: encapsulate return value in boost::shared_ptr or in VectorContent.
 * \return gsl_vector pointer to vector of eigenvalues
 */
gsl_vector* MatrixContent::transformToEigenbasis()
{
  if (rows == columns) { // symmetric
    gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
    gsl_vector *eval = gsl_vector_alloc(rows);
    gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
    gsl_eigen_symmv(content, eval, evec, T);
    gsl_eigen_symmv_free(T);
    gsl_matrix_memcpy(content, evec);
    gsl_matrix_free(evec);
    return eval;
  } else { // non-symmetric
    //ASSERT(false, "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");

    // blow up gsl_matrix in content to square matrix, fill other components with zero
    const size_t greaterDimension = rows > columns ? rows : columns;
    gsl_matrix *content_square = gsl_matrix_alloc(greaterDimension, greaterDimension);
    for (size_t i=0; i<greaterDimension; i++) {
      for (size_t j=0; j<greaterDimension; j++) {
        const double value = ((i < rows) && (j < columns)) ? gsl_matrix_get(content,i,j) : 0.;
        gsl_matrix_set(content_square, i,j, value);
      }
    }

    // show squared matrix
    MatrixContent *ContentSquare = new MatrixViewContent(gsl_matrix_submatrix(content,0,0,content->size1, content->size2));
    std::cout << "The squared matrix is " << *ContentSquare << std::endl;
    delete ContentSquare;

    // solve eigenvalue problem
    gsl_eigen_nonsymmv_workspace *T = gsl_eigen_nonsymmv_alloc(rows);
    gsl_vector_complex *eval = gsl_vector_complex_alloc(greaterDimension);
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(greaterDimension, greaterDimension);
    gsl_eigen_nonsymmv(content_square, eval, evec, T);
    gsl_eigen_nonsymmv_free(T);

    // copy eigenvectors real-parts into content_square and ...
    // ... show complex-valued eigenvector matrix
    std::cout << "Resulting eigenvector matrix is [";
    for (size_t i=0; i<greaterDimension; i++) {
      for (size_t j=0; j<greaterDimension; j++) {
        std::cout << "(" << GSL_REAL(gsl_matrix_complex_get(evec,i,j))
            << "," << GSL_IMAG(gsl_matrix_complex_get(evec,i,j)) << ")";
        if (j < greaterDimension-1)
          std::cout << " ";
        gsl_matrix_set(content_square, i,j, GSL_REAL(gsl_matrix_complex_get(evec,i,j)));
      }
      if (i < greaterDimension-1)
        std::cout << "; ";
    }
    std::cout << "]" << std::endl;

    // scan for zero rows and columns to drop
    std::set<size_t> RowDropList;
    std::set<size_t> ColumnDropList;
    for (size_t i=0; i<greaterDimension; i++) { // only copy real space part
      double Rsum = 0.;
      double Csum = 0.;
      for (size_t j=0; j<greaterDimension; j++) {
        Rsum += fabs(GSL_REAL(gsl_matrix_complex_get(evec,i,j)));
        Csum += fabs(GSL_REAL(gsl_matrix_complex_get(evec,j,i)));
      }
      Rsum /= (double)greaterDimension;
      Csum /= (double)greaterDimension;
      if (Rsum < MYEPSILON)
        RowDropList.insert(i);
      if (Csum < MYEPSILON)
        ColumnDropList.insert(i);
    }

    // copy real-parts of complex eigenvalues and eigenvectors
    gsl_vector *eval_real = gsl_vector_alloc(greaterDimension);
    size_t I=0;
    size_t J=0;
    for (size_t i=0; i<greaterDimension; i++) { // only copy real space part
      if (RowDropList.find(i) == RowDropList.end()) {
        for (size_t j=0; j<greaterDimension; j++) {
          if (ColumnDropList.find(j) == ColumnDropList.end()) {
            if (fabs(GSL_IMAG(gsl_matrix_complex_get(evec,i,j))) > MYEPSILON)
              std::cerr << "MatrixContent::transformToEigenbasis() - WARNING: eigenvectors are complex-valued!" << std::endl;
            gsl_matrix_set(content, I,J, GSL_REAL(gsl_matrix_complex_get(evec,i,j)));
            J++;
          }
        }
        if (fabs(GSL_IMAG(gsl_vector_complex_get(eval,I))) > MYEPSILON)
          std::cerr << "MatrixContent::transformToEigenbasis() - WARNING: eigenvectors are complex-valued!" << std::endl;
        gsl_vector_set(eval_real, I, GSL_REAL(gsl_vector_complex_get(eval, i)));
        I++;
      }
    }
    gsl_matrix_complex_free(evec);
    gsl_vector_complex_free(eval);
    return eval_real;
  }
}

/* ============================ Properties ============================== */
/** Checks whether matrix' elements are strictly null.
 * \return true - is null, false - else
 */
bool MatrixContent::IsNull() const
{
  return gsl_matrix_isnull (content);
};

/** Checks whether matrix' elements are strictly positive.
 * \return true - is positive, false - else
 */
bool MatrixContent::IsPositive() const
{
  return gsl_matrix_ispos (content);
};

/** Checks whether matrix' elements are strictly negative.
 * \return true - is negative, false - else
 */
bool MatrixContent::IsNegative() const
{
  return gsl_matrix_isneg (content);
};

/** Checks whether matrix' elements are strictly non-negative.
 * \return true - is non-negative, false - else
 */
bool MatrixContent::IsNonNegative() const
{
  return gsl_matrix_isnonneg (content);
};

/** This function performs a Cholesky decomposition to determine whether matrix is positive definite.
 * We check whether GSL returns GSL_EDOM as error, indicating that decomposition failed due to matrix not being positive-definite.
 * \return true - matrix is positive-definite, false - else
 */
bool MatrixContent::IsPositiveDefinite() const
{
  if (rows != columns)  // only possible for square matrices.
    return false;
  else
    return (gsl_linalg_cholesky_decomp (content) != GSL_EDOM);
};


/** Calculates the determinant of the matrix.
 * if matrix is square, uses LU decomposition.
 */
double MatrixContent::Determinant() const
{
  int signum = 0;
  assert (rows == columns && "Determinant can only be calculated for square matrices.");
  gsl_permutation *p = gsl_permutation_alloc(rows);
  gsl_linalg_LU_decomp(content, p, &signum);
  gsl_permutation_free(p);
  return gsl_linalg_LU_det(content, signum);
};

/* ============================= Operators =============================== */

/** Scalar multiplication operator.
 * \param factor factor to scale with
 */
const MatrixContent &MatrixContent::operator*=(const double factor)
{
  gsl_matrix_scale(content, factor);
  return *this;
}

/** Scalar multiplication and copy operator.
 * \param factor factor to scale with
 * \param &mat MatrixContent to scale
 * \return copied and scaled MatrixContent
 */
const MatrixContent operator*(const double factor,const MatrixContent& mat)
{
  MatrixContent tmp = mat;
  tmp*=factor;
  return tmp;
}

/** Scalar multiplication and copy operator (with operands exchanged).
 * \param &mat MatrixContent to scale
 * \param factor factor to scale with
 * \return copied and scaled MatrixContent
 */
const MatrixContent operator*(const MatrixContent &mat,const double factor)
{
  return factor*mat;
}

/** Equality operator.
 * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
 * \param &rhs MatrixContent to checks against
 */
bool MatrixContent::operator==(const MatrixContent &rhs) const
    {
  if ((rows == rhs.rows) && (columns == rhs.columns)) {
    for(int i=rows;i--;){
      for(int j=columns;j--;){
        if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
          return false;
        }
      }
    }
    return true;
  }
  return false;
}

Vector operator*(const MatrixContent &mat,const Vector &vec)
{
  Vector result;
  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, result.content->content);
  return result;
}

std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat)
{
  ost << "[";
  for (size_t i=0;i<mat.rows;i++) {
    for (size_t j=0;j<mat.columns;j++) {
      ost << mat.at(i,j);
      if (j != mat.columns-1)
        ost << " ";
    }
    if (i != mat.rows-1)
      ost << "; ";
  }
  ost << "]";
  return ost;
}