source: src/LinearAlgebra/LinearSystemOfEquations.cpp@ 1d5a871

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * LinearSystemOfEquations.cpp
 *
 *  Created on: Jan 8, 2010
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "Helpers/defs.hpp"
#include "LinearAlgebra/LinearSystemOfEquations.hpp"
#include "LinearAlgebra/MatrixContent.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/logger.hpp"
#include "LinearAlgebra/Vector.hpp"

#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 MatrixContent(m, n);
  b = new VectorContent(m);
  x = new VectorContent(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->setFromVector(*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->at(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[i] = x->at(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->content, p, &s);
    gsl_linalg_LU_solve (A->content, p, b->content, x->content);
    gsl_permutation_free (p);
  } else {  // use Householder
    //MatrixContent *backup = new MatrixContent(rows,columns);
    //*backup = *A;
    gsl_linalg_HH_solve (A->content, b->content, x->content);
    //*A = *backup;
    //delete(backup);
  }
  return true;
};