test_Spectrum_CRSMatrix_CRSsolver.cpp

Go to the documentation of this file.

#include "commonParameters.h"
#include "communicator.h"
#include "bridgeIO.h"
using Bridge::vout;

#include "field.h"
#include "fopr_CRS.h"
#include "solver_CG.h"

#ifdef USE_TEST
#include "test.h"
#endif

//     // clover set
//     string solution = "solution4x8_clover.crs";
//     string matrix   = "matrix4x8_clover.crs";
//     string source   = "source4x8_clover.crs";

//     // 5d-overlap set
//     string solution = "solution4x4_ov5d.crs";
//     string matrix   = "matrix4x4_ov5d.crs";
//     string source   = "source4x4_ov5d.crs";

//     // domain-wall set
//     string solution = "solution4x4_dw.crs";
//     string matrix   = "matrix4x4_dw.crs";
//     string source   = "source4x4_dw.crs";


namespace Test_Spectrum_CRSMatrix {
  int CRSsolver(
    const string& solution,
    const string& matrix,
    const string& source,
    double& result       /* return value */
    )
  {
    Bridge::VerboseLevel vl = CommonParameters::Vlevel();

    // #####  Following part is common  #####

    // read source and solution vectors
    Field b, xq;

    b.read_text(source);
    xq.read_text(solution);

    // read CRS matrix
    Fopr_CRS *fopr = new Fopr_CRS(matrix);


    // setup of CG solver
    int    Niter     = 2000;
    double Stop_cond = 1.0e-28;
    Solver *solver   = new Solver_CG((Fopr *)fopr);
    solver->set_parameters(Niter, Stop_cond);

    Field b2(b);
    Field x(b);

    // setup of CGNE source
    fopr->Ddag(b2, b);

    // CGNE solver
    int    Nconv;
    double diff;
    fopr->set_mode("DdagD");
    solver->solve(x, b2, Nconv, diff);
    vout.general(vl, "solver(CG):  Nconv = %4d  diff = %12.6e\n", Nconv, diff);

    // check
    x -= xq;
    double xx = x.norm2();
    vout.general(vl, "standard norm2 = %.8e\n", xx);

    result = xx;

    delete solver;
    delete fopr;

    return EXIT_SUCCESS;
  }
} // namespace Test_Spectrum_CRSMatrix