/*
 * 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_dirichlet.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<SolverDirichlet>();
  dsysv = new DSYSV<SolverDirichlet>();

  dgesv->x = new double[5];
  dsysv->x = new double[10];

  dgesv->vec_size = dgesv->max_size = 5;
  dsysv->vec_size = dsysv->max_size = 10;

  dgesv->b = new double[5];
  dsysv->b = new double[10];

  dgesv->A = new double*[5];
  dsysv->A = new double*[10];

  for (int i=0; i<5; ++i) {
    dgesv->A[i] = new double[5];
  }

  for (int i=0; i<10; ++i)
    dsysv->A[i] = new double[10];
#endif

  givens = new Givens<SolverDirichlet>();
  givens->vec_size = givens->max_size = 5;
  givens->x = new double[5];
  givens->b = new double[5];
  givens->A = new double*[5];

  for (int i=0; i<5; ++i) {
    givens->A[i] = new double[5];
  }
}

void SolverTestSuite::tearDown()
{
#ifdef HAVE_LAPACK
  delete dgesv;
  delete dsysv;
#endif
  delete givens;
}

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

  dgesv->A[1][0] = 0.0;
  dgesv->A[1][1] = 2.0;
  dgesv->A[1][2] = 2.0;
  dgesv->A[1][3] = 5.0;
  dgesv->A[1][4] = 4.0;

  dgesv->A[2][0] = 4.0;
  dgesv->A[2][1] = -16.0;
  dgesv->A[2][2] = 0.0;
  dgesv->A[2][3] = -25.0;
  dgesv->A[2][4] = 0.0;

  dgesv->A[3][0] = -1.0;
  dgesv->A[3][1] = -2.0;
  dgesv->A[3][2] = -4.0;
  dgesv->A[3][3] = -5.0;
  dgesv->A[3][4] = -8.0;

  dgesv->A[4][0] = 0.0;
  dgesv->A[4][1] = 0.0;
  dgesv->A[4][2] = 2.0;
  dgesv->A[4][3] = 0.0;
  dgesv->A[4][4] = 5.0;

  dgesv->b[0] = 12.5;
  dgesv->b[1] = 8.0;
  dgesv->b[2] = -10.0;
  dgesv->b[3] = 3.5;
  dgesv->b[4] = -16.0;

  dgesv->Compute();

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


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

  dsysv->A[1][0] = -1.0;
  dsysv->A[1][1] = 4.0;
  dsysv->A[1][2] = -1.0;
  dsysv->A[1][3] = 0.0;
  dsysv->A[1][4] = -1.0;
  dsysv->A[1][5] = 0.0;
  dsysv->A[1][6] = 0.0;
  dsysv->A[1][7] = -1.0;
  dsysv->A[1][8] = 0.0;
  dsysv->A[1][9] = 1.0;

  dsysv->A[2][0] = -1.0;
  dsysv->A[2][1] = -1.0;
  dsysv->A[2][2] = 4.0;
  dsysv->A[2][3] = 0.0;
  dsysv->A[2][4] = 0.0;
  dsysv->A[2][5] = -1.0;
  dsysv->A[2][6] = 0.0;
  dsysv->A[2][7] = 0.0;
  dsysv->A[2][8] = -1.0;
  dsysv->A[2][9] = 1.0;

  dsysv->A[3][0] = -1.0;
  dsysv->A[3][1] = 0.0;
  dsysv->A[3][2] = 0.0;
  dsysv->A[3][3] = 4.0;
  dsysv->A[3][4] = -1.0;
  dsysv->A[3][5] = -1.0;
  dsysv->A[3][6] = -1.0;
  dsysv->A[3][7] = 0.0;
  dsysv->A[3][8] = 0.0;
  dsysv->A[3][9] = 1.0;

  dsysv->A[4][0] = 0.0;
  dsysv->A[4][1] = -1.0;
  dsysv->A[4][2] = 0.0;
  dsysv->A[4][3] = -1.0;
  dsysv->A[4][4] = 4.0;
  dsysv->A[4][5] = -1.0;
  dsysv->A[4][6] = 0.0;
  dsysv->A[4][7] = -1.0;
  dsysv->A[4][8] = 0.0;
  dsysv->A[4][9] = 1.0;

  dsysv->A[5][0] = 0.0;
  dsysv->A[5][1] = 0.0;
  dsysv->A[5][2] = -1.0;
  dsysv->A[5][3] = -1.0;
  dsysv->A[5][4] = -1.0;
  dsysv->A[5][5] = 4.0;
  dsysv->A[5][6] = 0.0;
  dsysv->A[5][7] = 0.0;
  dsysv->A[5][8] = -1.0;
  dsysv->A[5][9] = 1.0;

  dsysv->A[6][0] = -1.0;
  dsysv->A[6][1] = 0.0;
  dsysv->A[6][2] = 0.0;
  dsysv->A[6][3] = -1.0;
  dsysv->A[6][4] = 0.0;
  dsysv->A[6][5] = 0.0;
  dsysv->A[6][6] = 4.0;
  dsysv->A[6][7] = -1.0;
  dsysv->A[6][8] = -1.0;
  dsysv->A[6][9] = 1.0;

  dsysv->A[7][0] = 0.0;
  dsysv->A[7][1] = -1.0;
  dsysv->A[7][2] = 0.0;
  dsysv->A[7][3] = 0.0;
  dsysv->A[7][4] = -1.0;
  dsysv->A[7][5] = 0.0;
  dsysv->A[7][6] = -1.0;
  dsysv->A[7][7] = 4.0;
  dsysv->A[7][8] = -1.0;
  dsysv->A[7][9] = 1.0;

  dsysv->A[8][0] = 0.0;
  dsysv->A[8][1] = 0.0;
  dsysv->A[8][2] = -1.0;
  dsysv->A[8][3] = 0.0;
  dsysv->A[8][4] = 0.0;
  dsysv->A[8][5] = -1.0;
  dsysv->A[8][6] = -1.0;
  dsysv->A[8][7] = -1.0;
  dsysv->A[8][8] = 4.0;
  dsysv->A[8][9] = 1.0;

  dsysv->A[9][0] = 1.0;
  dsysv->A[9][1] = 1.0;
  dsysv->A[9][2] = 1.0;
  dsysv->A[9][3] = 1.0;
  dsysv->A[9][4] = 1.0;
  dsysv->A[9][5] = 1.0;
  dsysv->A[9][6] = 1.0;
  dsysv->A[9][7] = 1.0;
  dsysv->A[9][8] = 1.0;
  dsysv->A[9][9] = 0.0;


  dsysv->b[0] = -0.125;
  dsysv->b[1] = -0.125;
  dsysv->b[2] = -0.125;
  dsysv->b[3] = -0.125;
  dsysv->b[4] = 1.0;
  dsysv->b[5] = -0.125;
  dsysv->b[6] = -0.125;
  dsysv->b[7] = -0.125;
  dsysv->b[8] = -0.125;
  dsysv->b[9] = 0.0;

  dsysv->Compute();

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

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

  givens->A[1][0] = 0.0;
  givens->A[1][1] = 2.0;
  givens->A[1][2] = 2.0;
  givens->A[1][3] = 5.0;
  givens->A[1][4] = 4.0;

  givens->A[2][0] = 4.0;
  givens->A[2][1] = -16.0;
  givens->A[2][2] = 0.0;
  givens->A[2][3] = -25.0;
  givens->A[2][4] = 0.0;

  givens->A[3][0] = -1.0;
  givens->A[3][1] = -2.0;
  givens->A[3][2] = -4.0;
  givens->A[3][3] = -5.0;
  givens->A[3][4] = -8.0;

  givens->A[4][0] = 0.0;
  givens->A[4][1] = 0.0;
  givens->A[4][2] = 2.0;
  givens->A[4][3] = 0.0;
  givens->A[4][4] = 5.0;

  givens->b[0] = 12.5;
  givens->b[1] = 8.0;
  givens->b[2] = -10.0;
  givens->b[3] = 3.5;
  givens->b[4] = -16.0;

  givens->Compute();

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