Index: src/LinearAlgebra/gslmatrix.cpp
===================================================================
--- src/LinearAlgebra/gslmatrix.cpp	(revision bcf653d9a0b46f52fc007e0fbb11d60cc4aa49a2)
+++ src/LinearAlgebra/gslmatrix.cpp	(revision cca9efb7bf18583f35d0380acb97e9034c42bbe7)
@@ -27,5 +27,10 @@
 
 #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.
@@ -46,4 +51,19 @@
   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;
 };
 
@@ -90,5 +110,5 @@
  * \return (m,n)-th element of matrix
  */
-double GSLMatrix::Get(size_t m, size_t n)
+double GSLMatrix::Get(size_t m, size_t n) const
 {
   return gsl_matrix_get (matrix, m, n);
@@ -121,8 +141,24 @@
  * \return const pointer to \a m-th element
  */
-const double *GSLMatrix::const_Pointer(size_t m, size_t n)
+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 =============================== */
@@ -204,5 +240,5 @@
  * \return true - is null, false - else
  */
-bool GSLMatrix::IsNull()
+bool GSLMatrix::IsNull() const
 {
   return gsl_matrix_isnull (matrix);
@@ -212,5 +248,5 @@
  * \return true - is positive, false - else
  */
-bool GSLMatrix::IsPositive()
+bool GSLMatrix::IsPositive() const
 {
   return gsl_matrix_ispos (matrix);
@@ -220,5 +256,5 @@
  * \return true - is negative, false - else
  */
-bool GSLMatrix::IsNegative()
+bool GSLMatrix::IsNegative() const
 {
   return gsl_matrix_isneg (matrix);
@@ -228,5 +264,5 @@
  * \return true - is non-negative, false - else
  */
-bool GSLMatrix::IsNonNegative()
+bool GSLMatrix::IsNonNegative() const
 {
   return gsl_matrix_isnonneg (matrix);
@@ -237,5 +273,5 @@
  * \return true - matrix is positive-definite, false - else
  */
-bool GSLMatrix::IsPositiveDefinite()
+bool GSLMatrix::IsPositiveDefinite() const
 {
   if (rows != columns)  // only possible for square matrices.
@@ -249,5 +285,5 @@
  * if matrix is square, uses LU decomposition.
  */
-double GSLMatrix::Determinant()
+double GSLMatrix::Determinant() const
 {
   int signum = 0;
@@ -259,2 +295,65 @@
 };
 
+
+/* ============================ 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;
+}
+
