source: src/LinearAlgebra/gslmatrix.cpp@ cca9ef

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

/*
 * gslmatrix.cpp
 *
 *  Created on: Jan 8, 2010
 *      Author: heber
 */

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

#include "Helpers/MemDebug.hpp"

using namespace std;

#include "LinearAlgebra/gslmatrix.hpp"
#include "Helpers/helpers.hpp"
#include "Helpers/fast_functions.hpp"

#include <cassert>
#include <iostream>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/** Constructor of class GSLMatrix.
 * Allocates GSL structures
 * \param m dimension of matrix
 */
GSLMatrix::GSLMatrix(size_t m, size_t n) : rows(m), columns(n)
{
  matrix = gsl_matrix_calloc(rows, columns);
};

/** Copy constructor of class GSLMatrix.
 * Allocates GSL structures and copies components from \a *src.
 * \param *src source matrix
 */
GSLMatrix::GSLMatrix(const GSLMatrix * const src) : rows(src->rows), columns(src->columns)
{
  matrix = gsl_matrix_alloc(rows, columns);
  gsl_matrix_memcpy (matrix, src->matrix);
};

/** Copy constructor of class GSLMatrix.
 * We take over pointer to gsl_matrix and set the parameter pointer to NULL
 * afterwards.
 * \param m row count
 * \param n column count
 * \param *&src source gsl_matrix
 */
GSLMatrix::GSLMatrix(size_t m, size_t n, gsl_matrix *&src) :
    rows(m),
    columns(n)
{
  matrix = src;
  src = NULL;
};

/** Destructor of class GSLMatrix.
 * Frees GSL structures
 */
GSLMatrix::~GSLMatrix()
{
  gsl_matrix_free(matrix);
  rows = 0;
  columns = 0;
};

/** Assignment operator.
 * \param &rhs right hand side
 * \return object itself
 */
GSLMatrix& GSLMatrix::operator=(const GSLMatrix& rhs)
{
  if (this == &rhs)   // not necessary here, but identity assignment check is generally a good idea
    return *this;
  assert(rows == rhs.rows && columns == rhs.columns && "Number of rows and columns do not match!");

  gsl_matrix_memcpy (matrix, rhs.matrix);

  return *this;
};

/* ============================ Accessing =============================== */
/** 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 GSLMatrix::SetFromDoubleArray(double * x)
{
  gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
  gsl_matrix_memcpy (matrix, &m.matrix);
};

/** This function returns the i-th element of a matrix.
 * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and 0 is returned.
 * \param m row index
 * \param n colum index
 * \return (m,n)-th element of matrix
 */
double GSLMatrix::Get(size_t m, size_t n) const
{
  return gsl_matrix_get (matrix, m, n);
};

/** This function sets the value of the \a m -th element of a matrix to \a x.
 *  If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked.
 * \param m row index
 * \param m column index
 * \param x value to set element (m,n) to
 */
void GSLMatrix::Set(size_t m, size_t n, double x)
{
  gsl_matrix_set (matrix, m, n, x);
};

/** 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 GSLMatrix::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 *GSLMatrix::Pointer(size_t m, size_t n)
{
  return gsl_matrix_ptr (matrix, 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 GSLMatrix::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 *GSLMatrix::const_Pointer(size_t m, size_t n) const
{
  return gsl_matrix_const_ptr (matrix, m, n);
};

/** Returns the number of rows of the matrix.
 * \return number of rows
 */
size_t GSLMatrix::getRowCount() const
{
  return rows;
}

/** Returns the number of columns of the matrix.
 * \return number of columns
 */
size_t GSLMatrix::getColumnCount() const
{
  return columns;
}

/* ========================== Initializing =============================== */
/** This function sets all the elements of the matrix to the value \a x.
 * \param *x
 */
void GSLMatrix::SetAll(double x)
{
  gsl_matrix_set_all (matrix, x);
};

/** This function sets all the elements of the matrix to zero.
 */
void GSLMatrix::SetZero()
{
  gsl_matrix_set_zero (matrix);
};

/** This function sets the elements of the matrix to the corresponding elements of the identity matrix, \f$m(i,j) = \delta(i,j)\f$, i.e. a unit diagonal with all off-diagonal elements zero.
 * This applies to both square and rectangular matrices.
 */
void GSLMatrix::SetIdentity()
{
  gsl_matrix_set_identity (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 GSLMatrix::SwapRows(size_t i, size_t j)
{
  return (gsl_matrix_swap_rows (matrix, 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 GSLMatrix::SwapColumns(size_t i, size_t j)
{
  return (gsl_matrix_swap_columns (matrix, 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 GSLMatrix::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 (matrix, i, j) == GSL_SUCCESS);
};

/** This function transposes the matrix.
 * Note that the function is extended to the non-square case, where the matrix is re-allocated and copied.
 */
bool GSLMatrix::Transpose()
{
  if (rows == columns)// if square, use GSL
    return (gsl_matrix_transpose (matrix) == GSL_SUCCESS);
  else { // otherwise we have to copy a bit
    gsl_matrix *dest = gsl_matrix_alloc(columns, rows);
    for (size_t i=0;i<rows; i++)
      for (size_t j=0;j<columns;j++) {
        gsl_matrix_set(dest, j,i, gsl_matrix_get(matrix, i,j) );
      }
    gsl_matrix_free(matrix);
    matrix = dest;
    flip(rows, columns);
    return true;
  }
};

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

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

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

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

/** 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 GSLMatrix::IsPositiveDefinite() const
{
  if (rows != columns)  // only possible for square matrices.
    return false;
  else
    return (gsl_linalg_cholesky_decomp (matrix) != GSL_EDOM);
};


/** Calculates the determinant of the matrix.
 * if matrix is square, uses LU decomposition.
 */
double GSLMatrix::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(matrix, p, &signum);
  gsl_permutation_free(p);
  return gsl_linalg_LU_det(matrix, signum); 
};


/* ============================ Properties ============================== */

const GSLMatrix &GSLMatrix::operator+=(const GSLMatrix &rhs)
{
  gsl_matrix_add(matrix, rhs.matrix);
  return *this;
}

const GSLMatrix &GSLMatrix::operator-=(const GSLMatrix &rhs)
{
  gsl_matrix_sub(matrix, rhs.matrix);
  return *this;
}

const GSLMatrix &GSLMatrix::operator*=(const GSLMatrix &rhs)
{
  (*this) = (*this)*rhs;
  return *this;
}

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

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

const GSLMatrix GSLMatrix::operator*(const GSLMatrix &rhs) const
{
  gsl_matrix *res = gsl_matrix_alloc(rhs.rows, rhs.columns);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix, rhs.matrix, 0.0, res);
  return GSLMatrix(rhs.rows, rhs.columns, res);
}

/** Print matrix in a matlab style manner.
 * \param &ost reference to output stream
 * \param &mat reference to matrix to print
 * \return reference to obtained output stream for concatenation
 */
ostream &operator<<(ostream &ost,const GSLMatrix &mat)
{
  ost << "[";
  for(size_t i = 0;i<mat.rows;++i){
    for(size_t j = 0; j<mat.columns;++j){
      ost << mat.Get(i,j);
      if(j!=mat.columns-1)
        ost << " ";
    }
    if(i!=mat.rows-1)
      ost << "; ";
  }
  ost << "]";
  return ost;
}