source: LinearAlgebra/src/LinearAlgebra/MatrixContent.hpp@ c5038e

/*
 * MatrixContent.hpp
 *
 *  Created on: Jul 2, 2010
 *      Author: crueger
 */

#ifndef MATRIXCONTENT_HPP_
#define MATRIXCONTENT_HPP_

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

#include "CodePatterns/Assert.hpp"

/** MatrixContent is a wrapper for gsl_matrix.
 *
 * The GNU Scientific Library contaisn some very well written routines for
 * linear algebra problems. However, it's syntax is C and hence it does not
 * lend itself to readable written code, i.e. C = A * B, where A, B, and C
 * are all matrices. Writing code this way is very convenient, concise and
 * also in the same manner as one would type in MatLab.
 * However, with C++ and its feature to overload functions and its operator
 * functions such syntax is easy to create.
 *
 * Hence, this class is a C++ wrapper to gsl_matrix. There already some other
 * libraries, however they fall short for the following reasons:
 * -# gslwrap: not very well written and uses floats instead of doubles
 * -# gslcpp: last change is from 2007 and only very few commits
 * -# o2scl: basically, the same functionality as gsl only written in C++,
 *    however it uses GPLv3 license which is inconvenient for MoleCuilder.
 *
 * <h1>Howto use MatrixContent</h1>
 *
 * One should not use MatrixContent directly but either have it as a member
 * variable in a specialized class or inherit from it.
 *
 */

#include <gsl/gsl_matrix.h>

#include <iosfwd>

#include "MatrixVector_ops.hpp"

class VectorContent;
namespace boost {
  namespace serialization {
    class access;
  }
}

/** Dummy structure to create a unique signature.
 *
 */
struct MatrixBaseCase{};

/** This class safely stores matrix dimensions away.
 *
 * This way we simulate const-ness of the dimensions while still allowing
 * serialization to modify them.
 *
 */
struct MatrixDimension {
public:
  MatrixDimension(size_t _rows, size_t _columns) : rows(_rows), columns(_columns) {}

  size_t getRows() const {return rows;}
  size_t getColumns() const {return columns;}

private:
  void setRows(size_t _rows) {rows = _rows;}
  void setColumns(size_t _columns) {columns = _columns;}

  friend class boost::serialization::access;
  // serialization (due to gsl_vector separate versions needed)
  template<class Archive>
  void serialize(Archive & ar, const unsigned int version)
  {
    ar & rows;
    ar & columns;
  }

private:
  size_t rows;      // store for internal purposes
  size_t columns;  // store for internal purposes
};

class MatrixContent : public MatrixDimension
{
  friend std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat);
  friend class RealSpaceMatrix;
  friend class LinearSystemOfEquations;

  // matrix vector products
  friend VectorContent const operator*(const VectorContent& vec, const MatrixContent& mat);
  friend VectorContent const operator*(const MatrixContent& mat, const VectorContent& vec);
  friend MatrixContent const operator+(const MatrixContent& a, const MatrixContent& b);
  friend MatrixContent const operator-(const MatrixContent& a, const MatrixContent& b);

  // unit tests
  friend class MatrixContentTest;
  friend class MatrixContentSymmetricTest;
public:
  MatrixContent(size_t rows, size_t columns);
  MatrixContent(size_t _rows, size_t _columns, MatrixBaseCase);
  MatrixContent(size_t rows, size_t columns, const double *src);
  MatrixContent(size_t rows, size_t columns, std::istream &inputstream);
  MatrixContent(gsl_matrix *&src);
  MatrixContent(const MatrixContent &src);
  MatrixContent(const MatrixContent *src);
  virtual ~MatrixContent();

  // Accessing
  double &at(size_t i, size_t j);
  const double at(size_t i, size_t j) const;
  void set(size_t i, size_t j, const double value);

  // Initializing
  void setFromDoubleArray(double * x);
  void setIdentity();
  void setValue(double _value);
  void setZero();

  // Exchanging elements
  bool SwapRows(size_t i, size_t j);
  bool SwapColumns(size_t i, size_t j);
  bool SwapRowColumn(size_t i, size_t j);

  // Transformations
  MatrixContent transposed() const;
  MatrixContent &transpose();
  gsl_vector* transformToEigenbasis();
  void sortEigenbasis(gsl_vector *eigenvalues);
  void performSingularValueDecomposition(MatrixContent &V, VectorContent &S);
  VectorContent solveOverdeterminedLinearEquation(const VectorContent &b);
  VectorContent solveOverdeterminedLinearEquation(MatrixContent &V, VectorContent &S, const VectorContent &b);

  // Properties
  bool IsNull() const;
  bool IsPositive() const;
  bool IsNegative() const;
  bool IsNonNegative() const;
  bool IsPositiveDefinite() const;
  double Determinant() const;

  // Operators
  MatrixContent &operator=(const MatrixContent &src);
  const MatrixContent &operator+=(const MatrixContent &rhs);
  const MatrixContent &operator-=(const MatrixContent &rhs);
  const MatrixContent &operator*=(const MatrixContent &rhs);
  const MatrixContent operator*(const MatrixContent &rhs) const;
  const MatrixContent &operator&=(const MatrixContent &rhs);
  const MatrixContent operator&(const MatrixContent &rhs) const;
  const MatrixContent &operator*=(const double factor);
  bool operator==(const MatrixContent &rhs) const;

  VectorContent *getColumnVector(size_t column) const;
  VectorContent *getRowVector(size_t row) const;
  VectorContent *getDiagonalVector() const;

  static std::pair<size_t, size_t> preparseMatrixDimensions(std::istream &inputstream);
  void write(std::ostream &ost) const;

protected:
  double *Pointer(size_t m, size_t n);
  const double *const_Pointer(size_t m, size_t n) const;

  gsl_matrix * content;

private:
  bool free_content_on_exit;
};

std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat);
const MatrixContent operator*(const double factor,const MatrixContent& mat);
const MatrixContent operator*(const MatrixContent &mat,const double factor);
MatrixContent const operator+(const MatrixContent& a, const MatrixContent& b);
MatrixContent const operator-(const MatrixContent& a, const MatrixContent& b);

/** Matrix view.
 * Extension of MatrixContent to contain not a gsl_matrix but only a view on a
 * gsl_matrix
 *
 * We need the above MatrixBaseCase here:
 * content, i.e. the gsl_matrix, must not be allocated, as it is just a view
 * with an internal gsl_matrix_view. Hence, MatrixBaseCase is just dummy class
 * to give the constructor a unique signature.
 */
struct MatrixViewContent : public MatrixContent
{
  MatrixViewContent(gsl_matrix_view _view) :
    MatrixContent(_view.matrix.size1, _view.matrix.size2, MatrixBaseCase()),
    view(_view)
  {
    content = &view.matrix;
  }
  ~MatrixViewContent(){
    content=0;
  }
  gsl_matrix_view view;
};

#endif /* MATRIXCONTENT_HPP_ */