Index: src/linearsystemofequations.cpp
===================================================================
--- src/linearsystemofequations.cpp	(revision f6061053393b67472b604538fdfdae9f6563328c)
+++ src/linearsystemofequations.cpp	(revision f6061053393b67472b604538fdfdae9f6563328c)
@@ -0,0 +1,129 @@
+/*
+ * linearsystemofequations.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#include "defs.hpp"
+#include "gslmatrix.hpp"
+#include "gslvector.hpp"
+#include "linearsystemofequations.hpp"
+#include "logger.hpp"
+#include "vector.hpp"
+
+#include <cassert>
+#include <gsl/gsl_permutation.h>
+
+/** Constructor for class LinearSystemOfEquations.
+ * Allocates Vector and Matrix classes.
+ * \param m column dimension
+ * \param n row dimension
+ */
+LinearSystemOfEquations::LinearSystemOfEquations(int m, int n) : rows(m), columns(n), IsSymmetric(false)
+{
+  A = new GSLMatrix(m, n);
+  b = new GSLVector(m);
+  x = new GSLVector(n);
+};
+
+/** Destructor for class LinearSystemOfEquations.
+ * Destructs Vector and Matrix classes.
+ */
+LinearSystemOfEquations::~LinearSystemOfEquations()
+{
+  delete(A);
+  delete(b);
+  delete(x);
+};
+
+/** Sets whether matrix is to be regarded as symmetric.
+ * Note that we do not check whether it really is, just take upper diagonal.
+ * \param symmetric true or false
+ */
+bool LinearSystemOfEquations::SetSymmetric(bool symmetric)
+{
+  assert (rows == columns && "Rows and columns don't have equal size! Setting symmetric not possible.");
+  return (IsSymmetric = symmetric);
+};
+
+/** Initializes vector b to the components of the given vector.
+ * \param *x Vector with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::Setb(Vector *x)
+{
+  assert ( columns == NDIM && "Vector class is always three-dimensional, unlike this LEqS!");
+  b->SetFromDoubleArray(x->x);
+};
+
+/** Initializes vector b to the components of the given vector.
+ * \param *x array with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::Setb(double *x)
+{
+  b->SetFromDoubleArray(x);
+};
+
+/** Initializes matrix a to the components of the given array.
+ * note that sort order should be
+ * \param *x array with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::SetA(double *x)
+{
+  A->SetFromDoubleArray(x);
+};
+
+/** Returns the solution vector x \f$A \cdot x = b\f$ as an array of doubles.
+ * \param *array pointer allocated array for solution on return (no bounds check, dimension must match)
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::GetSolutionAsArray(double *&array)
+{
+  bool status = Solve();
+
+  // copy solution
+  for (size_t i=0;i<x->dimension;i++) {
+    array[i] = x->Get(i);
+  }
+  return status;
+};
+
+/** Returns the solution vector x \f$A \cdot x = b\f$ as an array of doubles.
+ * \param &x solution vector on return (must be 3-dimensional)
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::GetSolutionAsVector(Vector &v)
+{
+  assert(rows == NDIM && "Solution can only be returned as vector if number of columns is NDIM.");
+  bool status = Solve();
+
+  // copy solution
+  for (size_t i=0;i<x->dimension;i++)
+    v.x[i] = x->Get(i);
+  return status;
+};
+
+/** Solves the given system of \f$A \cdot x = b\f$.
+ * Use either LU or Householder decomposition.
+ * Solution is stored in LinearSystemOfEquations::x
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::Solve()
+{
+  // calculate solution
+  int s;
+  if (IsSymmetric) { // use LU
+    gsl_permutation * p = gsl_permutation_alloc (x->dimension);
+    gsl_linalg_LU_decomp (A->matrix, p, &s);
+    gsl_linalg_LU_solve (A->matrix, p, b->vector, x->vector);
+    gsl_permutation_free (p);
+  } else {  // use Householder
+    //GSLMatrix *backup = new GSLMatrix(rows,columns);
+    //*backup = *A;
+    gsl_linalg_HH_solve (A->matrix, b->vector, x->vector);
+    //*A = *backup;
+    //delete(backup);
+  }
+  return true;
+};
+
