source: src/LinearAlgebra/gslmatrix.cpp@ 670cd1

/*
 * 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 <gsl/gsl_linalg.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);
};

/** 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)
{
  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)
{
  return gsl_matrix_const_ptr (matrix, m, n);
};

/* ========================== 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()
{
  return gsl_matrix_isnull (matrix);
};

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

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

/** Checks whether matrix' elements are strictly non-negative.
 * \return true - is non-negative, false - else
 */
bool GSLMatrix::IsNonNegative()
{
  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()
{
  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()
{
  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); 
};