/*
 * solver_test.cpp
 *
 *  Created on: 21.07.2010
 *      Author: Julian Iseringhausen
 */

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

#include <limits>

#ifdef HAVE_LAPACK
#include "solver/dgesv.hpp"
#include "solver/dsysv.hpp"
#endif
#include "solver/givens.hpp"
#include "solver/solver_regular.hpp"

#include "solver_test.hpp"

using namespace VMG;
using VMGTests::SolverTestSuite;

CPPUNIT_TEST_SUITE_REGISTRATION(SolverTestSuite);

void SolverTestSuite::setUp()
{
#ifdef HAVE_LAPACK
  dgesv = new DGESV<SolverRegular>();
  dsysv = new DSYSV<SolverRegular>();

  dgesv->Realloc(5);
  dsysv->Realloc(10);
#endif

  givens = new Givens<SolverRegular>();
  givens->Realloc(5);
}

void SolverTestSuite::tearDown()
{
  MG::Destroy();

#ifdef HAVE_LAPACK
  delete dgesv;
  delete dsysv;
#endif
  delete givens;
}

void SolverTestSuite::DGESVTest()
{
#ifdef HAVE_LAPACK
  dgesv->Mat(0, 0) = 1.0;
  dgesv->Mat(0, 1) = -4.0;
  dgesv->Mat(0, 2) = 0.0;
  dgesv->Mat(0, 3) = -5.0;
  dgesv->Mat(0, 4) = 0.0;

  dgesv->Mat(1, 0) = 0.0;
  dgesv->Mat(1, 1) = 2.0;
  dgesv->Mat(1, 2) = 2.0;
  dgesv->Mat(1, 3) = 5.0;
  dgesv->Mat(1, 4) = 4.0;

  dgesv->Mat(2, 0) = 4.0;
  dgesv->Mat(2, 1) = -16.0;
  dgesv->Mat(2, 2) = 0.0;
  dgesv->Mat(2, 3) = -25.0;
  dgesv->Mat(2, 4) = 0.0;

  dgesv->Mat(3, 0) = -1.0;
  dgesv->Mat(3, 1) = -2.0;
  dgesv->Mat(3, 2) = -4.0;
  dgesv->Mat(3, 3) = -5.0;
  dgesv->Mat(3, 4) = -8.0;

  dgesv->Mat(4, 0) = 0.0;
  dgesv->Mat(4, 1) = 0.0;
  dgesv->Mat(4, 2) = 2.0;
  dgesv->Mat(4, 3) = 0.0;
  dgesv->Mat(4, 4) = 5.0;

  dgesv->Rhs(0) = 12.5;
  dgesv->Rhs(1) = 8.0;
  dgesv->Rhs(2) = -10.0;
  dgesv->Rhs(3) = 3.5;
  dgesv->Rhs(4) = -16.0;

  dgesv->Compute();

  CPPUNIT_ASSERT_DOUBLES_EQUAL(  8.5, dgesv->Sol(0), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-16.0, dgesv->Sol(1), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-18.0, dgesv->Sol(2), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( 12.0, dgesv->Sol(3), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(  4.0, dgesv->Sol(4), 1.0e-12);
#endif
}


void SolverTestSuite::DSYSVTest()
{
#ifdef HAVE_LAPACK
  dsysv->Mat(0, 0) = 4.0;
  dsysv->Mat(0, 1) = -1.0;
  dsysv->Mat(0, 2) = -1.0;
  dsysv->Mat(0, 3) = -1.0;
  dsysv->Mat(0, 4) = 0.0;
  dsysv->Mat(0, 5) = 0.0;
  dsysv->Mat(0, 6) = -1.0;
  dsysv->Mat(0, 7) = 0.0;
  dsysv->Mat(0, 8) = 0.0;
  dsysv->Mat(0, 9) = 1.0;

  dsysv->Mat(1, 0) = -1.0;
  dsysv->Mat(1, 1) = 4.0;
  dsysv->Mat(1, 2) = -1.0;
  dsysv->Mat(1, 3) = 0.0;
  dsysv->Mat(1, 4) = -1.0;
  dsysv->Mat(1, 5) = 0.0;
  dsysv->Mat(1, 6) = 0.0;
  dsysv->Mat(1, 7) = -1.0;
  dsysv->Mat(1, 8) = 0.0;
  dsysv->Mat(1, 9) = 1.0;

  dsysv->Mat(2, 0) = -1.0;
  dsysv->Mat(2, 1) = -1.0;
  dsysv->Mat(2, 2) = 4.0;
  dsysv->Mat(2, 3) = 0.0;
  dsysv->Mat(2, 4) = 0.0;
  dsysv->Mat(2, 5) = -1.0;
  dsysv->Mat(2, 6) = 0.0;
  dsysv->Mat(2, 7) = 0.0;
  dsysv->Mat(2, 8) = -1.0;
  dsysv->Mat(2, 9) = 1.0;

  dsysv->Mat(3, 0) = -1.0;
  dsysv->Mat(3, 1) = 0.0;
  dsysv->Mat(3, 2) = 0.0;
  dsysv->Mat(3, 3) = 4.0;
  dsysv->Mat(3, 4) = -1.0;
  dsysv->Mat(3, 5) = -1.0;
  dsysv->Mat(3, 6) = -1.0;
  dsysv->Mat(3, 7) = 0.0;
  dsysv->Mat(3, 8) = 0.0;
  dsysv->Mat(3, 9) = 1.0;

  dsysv->Mat(4, 0) = 0.0;
  dsysv->Mat(4, 1) = -1.0;
  dsysv->Mat(4, 2) = 0.0;
  dsysv->Mat(4, 3) = -1.0;
  dsysv->Mat(4, 4) = 4.0;
  dsysv->Mat(4, 5) = -1.0;
  dsysv->Mat(4, 6) = 0.0;
  dsysv->Mat(4, 7) = -1.0;
  dsysv->Mat(4, 8) = 0.0;
  dsysv->Mat(4, 9) = 1.0;

  dsysv->Mat(5, 0) = 0.0;
  dsysv->Mat(5, 1) = 0.0;
  dsysv->Mat(5, 2) = -1.0;
  dsysv->Mat(5, 3) = -1.0;
  dsysv->Mat(5, 4) = -1.0;
  dsysv->Mat(5, 5) = 4.0;
  dsysv->Mat(5, 6) = 0.0;
  dsysv->Mat(5, 7) = 0.0;
  dsysv->Mat(5, 8) = -1.0;
  dsysv->Mat(5, 9) = 1.0;

  dsysv->Mat(6, 0) = -1.0;
  dsysv->Mat(6, 1) = 0.0;
  dsysv->Mat(6, 2) = 0.0;
  dsysv->Mat(6, 3) = -1.0;
  dsysv->Mat(6, 4) = 0.0;
  dsysv->Mat(6, 5) = 0.0;
  dsysv->Mat(6, 6) = 4.0;
  dsysv->Mat(6, 7) = -1.0;
  dsysv->Mat(6, 8) = -1.0;
  dsysv->Mat(6, 9) = 1.0;

  dsysv->Mat(7, 0) = 0.0;
  dsysv->Mat(7, 1) = -1.0;
  dsysv->Mat(7, 2) = 0.0;
  dsysv->Mat(7, 3) = 0.0;
  dsysv->Mat(7, 4) = -1.0;
  dsysv->Mat(7, 5) = 0.0;
  dsysv->Mat(7, 6) = -1.0;
  dsysv->Mat(7, 7) = 4.0;
  dsysv->Mat(7, 8) = -1.0;
  dsysv->Mat(7, 9) = 1.0;

  dsysv->Mat(8, 0) = 0.0;
  dsysv->Mat(8, 1) = 0.0;
  dsysv->Mat(8, 2) = -1.0;
  dsysv->Mat(8, 3) = 0.0;
  dsysv->Mat(8, 4) = 0.0;
  dsysv->Mat(8, 5) = -1.0;
  dsysv->Mat(8, 6) = -1.0;
  dsysv->Mat(8, 7) = -1.0;
  dsysv->Mat(8, 8) = 4.0;
  dsysv->Mat(8, 9) = 1.0;

  dsysv->Mat(9, 0) = 1.0;
  dsysv->Mat(9, 1) = 1.0;
  dsysv->Mat(9, 2) = 1.0;
  dsysv->Mat(9, 3) = 1.0;
  dsysv->Mat(9, 4) = 1.0;
  dsysv->Mat(9, 5) = 1.0;
  dsysv->Mat(9, 6) = 1.0;
  dsysv->Mat(9, 7) = 1.0;
  dsysv->Mat(9, 8) = 1.0;
  dsysv->Mat(9, 9) = 0.0;


  dsysv->Rhs(0) = -0.125;
  dsysv->Rhs(1) = -0.125;
  dsysv->Rhs(2) = -0.125;
  dsysv->Rhs(3) = -0.125;
  dsysv->Rhs(4) = 1.0;
  dsysv->Rhs(5) = -0.125;
  dsysv->Rhs(6) = -0.125;
  dsysv->Rhs(7) = -0.125;
  dsysv->Rhs(8) = -0.125;
  dsysv->Rhs(9) = 0.0;

  dsysv->Compute();

  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->Sol(0), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(    0.0, dsysv->Sol(1), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->Sol(2), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(    0.0, dsysv->Sol(3), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(   0.25, dsysv->Sol(4), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(    0.0, dsysv->Sol(5), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->Sol(6), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(    0.0, dsysv->Sol(7), std::numeric_limits<double>::epsilon());
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->Sol(8), std::numeric_limits<double>::epsilon());
#endif
}

void SolverTestSuite::GivensTest()
{
  givens->Mat(0, 0) = 1.0;
  givens->Mat(0, 1) = -4.0;
  givens->Mat(0, 2) = 0.0;
  givens->Mat(0, 3) = -5.0;
  givens->Mat(0, 4) = 0.0;

  givens->Mat(1, 0) = 0.0;
  givens->Mat(1, 1) = 2.0;
  givens->Mat(1, 2) = 2.0;
  givens->Mat(1, 3) = 5.0;
  givens->Mat(1, 4) = 4.0;

  givens->Mat(2, 0) = 4.0;
  givens->Mat(2, 1) = -16.0;
  givens->Mat(2, 2) = 0.0;
  givens->Mat(2, 3) = -25.0;
  givens->Mat(2, 4) = 0.0;

  givens->Mat(3, 0) = -1.0;
  givens->Mat(3, 1) = -2.0;
  givens->Mat(3, 2) = -4.0;
  givens->Mat(3, 3) = -5.0;
  givens->Mat(3, 4) = -8.0;

  givens->Mat(4, 0) = 0.0;
  givens->Mat(4, 1) = 0.0;
  givens->Mat(4, 2) = 2.0;
  givens->Mat(4, 3) = 0.0;
  givens->Mat(4, 4) = 5.0;

  givens->Rhs(0) = 12.5;
  givens->Rhs(1) = 8.0;
  givens->Rhs(2) = -10.0;
  givens->Rhs(3) = 3.5;
  givens->Rhs(4) = -16.0;

  givens->Compute();

  CPPUNIT_ASSERT_DOUBLES_EQUAL(  8.5, givens->Sol(0), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-16.0, givens->Sol(1), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(-18.0, givens->Sol(2), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( 12.0, givens->Sol(3), 1.0e-12);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(  4.0, givens->Sol(4), 1.0e-12);
}