test implementation of DWF inverter

programs/Hadrons/Application.cc

@@ -156,6 +156,7 @@ void Application::configLoop(void)
     {
         LOG(Message) << "Starting measurement for trajectory " << t
                      << std::endl;
+        env_.loadUnitGauge();
         execute(program_);
         env_.freeAll();
     }

programs/Hadrons/Environment.cc

@@ -29,9 +29,15 @@ using namespace Hadrons;
 // constructor /////////////////////////////////////////////////////////////////
 Environment::Environment(void)
 {
+    std::vector<int> seed4d({1,2,3,4});
+    
     grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
         GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
         GridDefaultMpi()));
+    gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
+    rng4d_.reset(new GridParallelRNG(grid4d_.get()));
+    rng4d_->SeedFixedIntegers(seed4d);
+    gauge_.reset(new LatticeGaugeField(grid4d_.get()));
 }
 
 // dry run /////////////////////////////////////////////////////////////////////
@@ -45,30 +51,67 @@ bool Environment::isDryRun(void)
     return dryRun_;
 }
 
+// grids ///////////////////////////////////////////////////////////////////////
+GridCartesian * Environment::get4dGrid(void)
+{
+    return grid4d_.get();
+}
+
+GridRedBlackCartesian * Environment::getRb4dGrid(void)
+{
+    return gridRb4d_.get();
+}
+
+GridCartesian * Environment::get5dGrid(const unsigned int Ls)
+{
+    try
+    {
+        return grid5d_.at(Ls).get();
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR("no 5D grid with Ls= " << Ls);
+    }
+}
+
+GridRedBlackCartesian * Environment::getRb5dGrid(const unsigned int Ls)
+{
+    try
+    {
+        return gridRb5d_.at(Ls).get();
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR("no red-black 5D grid with Ls= " << Ls);
+    }
+}
+
 // quark propagators ///////////////////////////////////////////////////////////
 void Environment::addProp(const std::string name, const unsigned int Ls)
 {
+    GridCartesian *p4 = grid4d_.get();
+
     if (propExists(name))
     {
         HADRON_ERROR("propagator '" + name + "' already exists");
     }
     if (Ls > 1)
     {
-        GridCartesian *pt;
+        GridCartesian *p;
         
         try
         {
-            pt = grid5d_.at(Ls).get();
+            p = grid5d_.at(Ls).get();
         }
         catch(std::out_of_range &)
         {
-            grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls,
+            grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, p4));
-                                                             grid4d_.get()));
+            gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, p4));
-            pt = grid5d_[Ls].get();
+            p = grid5d_[Ls].get();
         }
         if (!isDryRun())
         {
-            prop_[name].reset(new LatticePropagator(pt));
+            prop_[name].reset(new LatticePropagator(p));
         }
         else
         {
@@ -80,7 +123,7 @@ void Environment::addProp(const std::string name, const unsigned int Ls)
     {
         if (!isDryRun())
         {
-            prop_[name].reset(new LatticePropagator(grid4d_.get()));
+            prop_[name].reset(new LatticePropagator(p4));
         }
         else
         {
@@ -143,6 +186,22 @@ unsigned int Environment::nProp(void)
     return size;
 }
 
+// gauge configuration /////////////////////////////////////////////////////////
+LatticeGaugeField * Environment::getGauge(void)
+{
+    return gauge_.get();
+}
+
+void Environment::loadUnitGauge(void)
+{
+    SU3::ColdConfiguration(*rng4d_, *gauge_);
+}
+
+void Environment::loadRandomGauge(void)
+{
+    SU3::HotConfiguration(*rng4d_, *gauge_);
+}
+
 // general free ////////////////////////////////////////////////////////////////
 void Environment::free(const std::string name)
 {

programs/Hadrons/Environment.hpp

@@ -31,28 +31,45 @@ class Environment
 {
     SINGLETON(Environment);
 public:
-    typedef std::unique_ptr<GridCartesian>     GridPt;
+    typedef std::unique_ptr<GridCartesian>         GridPt;
-    typedef std::unique_ptr<LatticePropagator> PropPt;
+    typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
+    typedef std::unique_ptr<GridParallelRNG>       RngPt;
+    typedef std::unique_ptr<LatticePropagator>     PropPt;
+    typedef std::unique_ptr<LatticeGaugeField>     GaugePt;
 public:
     // dry run
-    void                dryRun(const bool isDry);
+    void                    dryRun(const bool isDry);
-    bool                isDryRun(void);
+    bool                    isDryRun(void);
+    // grids
+    GridCartesian *         get4dGrid(void);
+    GridRedBlackCartesian * getRb4dGrid(void);
+    GridCartesian *         get5dGrid(const unsigned int Ls);
+    GridRedBlackCartesian * getRb5dGrid(const unsigned int Ls);
     // quark propagators
-    void                addProp(const std::string name,
+    void                    addProp(const std::string name,
-                                const unsigned int Ls = 1);
+                                    const unsigned int Ls = 1);
-    void                freeProp(const std::string name);
+    void                    freeProp(const std::string name);
-    LatticePropagator * getProp(const std::string name);
+    LatticePropagator *     getProp(const std::string name);
-    bool                propExists(const std::string name);
+    bool                    propExists(const std::string name);
-    unsigned int        nProp(void);
+    unsigned int            nProp(void);
+    // gauge configuration
+    LatticeGaugeField *     getGauge(void);
+    void                    loadUnitGauge(void);
+    void                    loadRandomGauge(void);
     // general free
-    void                free(const std::string name);
+    void                    free(const std::string name);
-    void                freeAll(void);
+    void                    freeAll(void);
 private:
     bool                                dryRun_{false};
     GridPt                              grid4d_;
     std::map<unsigned int, GridPt>      grid5d_;
+    GridRbPt                            gridRb4d_;
+    std::map<unsigned int, GridRbPt>    gridRb5d_;
+    RngPt                               rng4d_;
     std::map<std::string, PropPt>       prop_;
     std::map<std::string, unsigned int> propSize_;
+    GaugePt                             gauge_;
+    
 };
 
 END_HADRONS_NAMESPACE

programs/Hadrons/MQuark.cc

@@ -66,4 +66,20 @@ void MQuark::execute(Environment &env)
 {
     LOG(Message) << "computing quark propagator '" << getName() << "'"
                  << std::endl;
+    
+    GridCartesian         *g4d   = env.get4dGrid(),
+                          *g5d   = env.get5dGrid(par_.Ls);
+    GridRedBlackCartesian *gRb4d = env.getRb4dGrid(),
+                          *gRb5d = env.getRb5dGrid(par_.Ls);
+    LatticeGaugeField     &Umu   = *env.getGauge();
+    LatticeFermion        src(g5d); src=zero;
+    LatticeFermion        result(g5d); result=zero;
+    
+    RealD mass=0.1;
+    RealD M5=1.8;
+    DomainWallFermionR Ddwf(Umu, *g5d, *gRb5d, *g4d, *gRb4d, mass, M5);
+    
+    ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+    SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
+    SchurSolver(Ddwf,src,result);
 }