Line data    Source code

       1             : /*
       2             :  * File:   TestCGSolver.cpp
       3             :  * Author: chung
       4             :  *
       5             :  * Created on Nov 16, 2015, 6:41:56 PM
       6             :  */
       7             : 
       8             : #include "TestCGSolver.h"
       9             : #include "CGSolver.h"
      10             : #include <iostream>
      11             : 
      12           1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestCGSolver);
      13             : 
      14           3 : TestCGSolver::TestCGSolver() {
      15           3 : }
      16             : 
      17           6 : TestCGSolver::~TestCGSolver() {
      18           6 : }
      19             : 
      20           3 : void TestCGSolver::setUp() {
      21           3 : }
      22             : 
      23           3 : void TestCGSolver::tearDown() {
      24           3 : }
      25             : 
      26           1 : void TestCGSolver::testSolve() {
      27           1 :     size_t n = 800;
      28             : 
      29           1 :     Matrix b = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
      30           2 :     Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_SYMMETRIC);
      31           2 :     Matrix Y = MatrixFactory::MakeIdentity(n, 50.0);
      32           1 :     A += Y;
      33             : 
      34           2 :     Matrix ID(n, n, Matrix::MATRIX_DIAGONAL);
      35         801 :     for (size_t j = 0; j < n; ++j) {
      36         800 :         ID.set(j, j, 1 / A.get(j, j));
      37             :     }
      38           2 :     MatrixOperator Aop(A);
      39           2 :     MatrixOperator M(ID);
      40             : 
      41           1 :     size_t max_iter = n;
      42           1 :     const double tolerance = 1e-4;
      43           2 :     CGSolver solver(Aop, M, tolerance, max_iter);
      44           2 :     Matrix sol(n, 1);
      45           1 :     int status = solver.solve(b, sol);
      46           1 :     _ASSERT(ForBESUtils::is_status_ok(status));
      47           1 :     _ASSERT(solver.last_num_iter() < max_iter);
      48           1 :     _ASSERT(solver.last_error() < tolerance);
      49           1 :     _ASSERT(solver.last_num_iter() > 0);
      50             : 
      51           2 :     Matrix Asol = Aop.call(sol);
      52         801 :     for (size_t i = 0; i < n; i++) {
      53         800 :         _ASSERT_NUM_EQ(Asol[i], b[i], tolerance);
      54           1 :     }
      55             : 
      56           1 : }
      57             : 
      58           1 : void TestCGSolver::testSolve2() {
      59           1 :     size_t n = 500;
      60             : 
      61           1 :     Matrix b = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
      62           2 :     Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_SYMMETRIC);
      63             : 
      64           2 :     Matrix ID = MatrixFactory::MakeIdentity(n, 1.0);
      65           2 :     MatrixOperator Aop(A);
      66           2 :     MatrixOperator M(ID);
      67             : 
      68           1 :     size_t max_iter = 5;
      69           2 :     CGSolver solver(Aop, M, 1e-4, max_iter);
      70           2 :     Matrix sol(n, 1);
      71           1 :     int status = solver.solve(b, sol);
      72           1 :     _ASSERT_EQ(ForBESUtils::STATUS_MAX_ITERATIONS_REACHED, status);
      73           2 :     _ASSERT_EQ(max_iter, solver.last_num_iter());
      74           1 : }
      75             : 
      76           1 : void TestCGSolver::testSolveNoPredcond() {
      77           1 :     size_t n = 500;
      78             : 
      79           1 :     Matrix b = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
      80           2 :     Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_SYMMETRIC);
      81             : 
      82           2 :     MatrixOperator Aop(A);
      83             : 
      84           2 :     CGSolver solver(Aop);
      85           2 :     Matrix sol(n, 1);
      86           1 :     int status = solver.solve(b, sol);
      87           1 :     _ASSERT(ForBESUtils::is_status_ok(status));
      88           1 :     const double default_tolerance = 1e-4;
      89           2 :     _ASSERT(solver.last_error() < default_tolerance);
      90           4 : }
      91             :