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