source: src/LinearAlgebra/MatrixContent.cpp@ bbf1bd

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

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 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, 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.
 * \warn return vector has to be freed'd
 * \return gsl_vector pointer to vector of eigenvalues
 */
gsl_vector* MatrixContent::transformToEigenbasis()
{
  ASSERT (rows == columns,
      "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
  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;
}

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