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

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


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

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

class MatrixContent
{
  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;

  const size_t getRows() const;
  const size_t getColumns() 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;
  const size_t rows;      // store for internal purposes
  const size_t columns;  // store for internal purposes

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