Merge branch 'fix/collections-autotuning' into 'master'

Fix/collections autotuning

Contrary to its name, this MR actually sets the code to use Collections by default. Current policy is to use StdMat approach when P<5 and IterPerExp otherwise. This can still be overridden in the XML file and autotuning can optionally be turned on if desired.

See merge request !476

Author: mdave

library/Collections/Collection.cpp

@@ -57,19 +57,20 @@ Collection::Collection(
         OperatorType opType = (OperatorType)i;
         ImplementationType impType;
 
-        it = impTypes.find(opType);
-        impType = it == impTypes.end() ? eIterPerExp : it->second;
+        if ((it = impTypes.find(opType)) != impTypes.end())
+        {
+            impType = it->second;
+            OperatorKey opKey(pCollExp[0]->DetShapeType(), opType, impType,
+                              pCollExp[0]->IsNodalNonTensorialExp());
 
-        OperatorKey opKey(pCollExp[0]->DetShapeType(), opType, impType,
-                          pCollExp[0]->IsNodalNonTensorialExp());
-
-        stringstream ss;
-        ss << opKey;
-        ASSERTL0(GetOperatorFactory().ModuleExists(opKey),
+            stringstream ss;
+            ss << opKey;
+            ASSERTL0(GetOperatorFactory().ModuleExists(opKey),
                  "Requested unknown operator "+ss.str());
 
-        m_ops[opType] = GetOperatorFactory().CreateInstance(
+            m_ops[opType] = GetOperatorFactory().CreateInstance(
                                                 opKey, pCollExp, m_geomData);
+        }
     }
 }
 

library/Collections/Collection.h

@@ -77,6 +77,8 @@ class Collection
                       Array<OneD,       NekDouble>           &output1,
                       Array<OneD,       NekDouble>           &output2);
 
+        inline bool HasOperator(const OperatorType &op);
+
     protected:
         StdRegions::StdExpansionSharedPtr                     m_stdExp;
         vector<SpatialDomains::GeometrySharedPtr>             m_geom;
@@ -131,6 +133,11 @@ inline void Collection::ApplyOperator(
     (*m_ops[op])(inarray, output0, output1, output2, wsp);
 }
 
+inline bool Collection::HasOperator(const OperatorType &op)
+{
+    return (m_ops.find(op) != m_ops.end());
+}
+
 }
 }
 

library/Collections/CollectionOptimisation.cpp

@@ -50,6 +50,7 @@ CollectionOptimisation::CollectionOptimisation(
 {
     int i;
     map<ElmtOrder, ImplementationType> defaults;
+    map<ElmtOrder, ImplementationType> defaultsPhysDeriv;
     map<ElmtOrder, ImplementationType>::iterator it;
     bool verbose  = (pSession.get()) &&
                     (pSession->DefinesCmdLineArgument("verbose")) &&
@@ -58,7 +59,7 @@ CollectionOptimisation::CollectionOptimisation(
     m_setByXml    = false;
     m_autotune    = false;
     m_maxCollSize = 0;
-    m_defaultType = defaultType == eNoImpType ? eNoCollection : defaultType;
+    m_defaultType = defaultType == eNoImpType ? eIterPerExp : defaultType;
 
     map<string, LibUtilities::ShapeType> elTypes;
     map<string, LibUtilities::ShapeType>::iterator it2;
@@ -73,14 +74,38 @@ CollectionOptimisation::CollectionOptimisation(
     // Set defaults for all element types.
     for (it2 = elTypes.begin(); it2 != elTypes.end(); ++it2)
     {
-        defaults[ElmtOrder(it2->second, -1)] = m_defaultType;
+        defaults          [ElmtOrder(it2->second, -1)] = m_defaultType;
+        defaultsPhysDeriv [ElmtOrder(it2->second, -1)] = m_defaultType;
+    }
+
+    if (defaultType == eNoImpType)
+    {
+        for (it2 = elTypes.begin(); it2 != elTypes.end(); ++it2)
+        {
+            defaultsPhysDeriv [ElmtOrder(it2->second, -1)] = eNoCollection;
+            for (int i = 1; i < 5; ++i)
+            {
+                defaults[ElmtOrder(it2->second, i)] = eStdMat;
+            }
+            for (int i = 1; i < 3; ++i)
+            {
+                defaultsPhysDeriv[ElmtOrder(it2->second, i)] = eSumFac;
+            }
+        }
     }
 
     map<string, OperatorType> opTypes;
     for (i = 0; i < SIZE_OperatorType; ++i)
     {
         opTypes[OperatorTypeMap[i]] = (OperatorType)i;
-        m_global[(OperatorType)i] = defaults;
+        switch ((OperatorType)i)
+        {
+            case ePhysDeriv:
+                m_global[(OperatorType)i] = defaultsPhysDeriv;
+                break;
+            default:
+                m_global[(OperatorType)i] = defaults;
+        }
     }
 
     map<string, ImplementationType> impTypes;
@@ -98,17 +123,23 @@ CollectionOptimisation::CollectionOptimisation(
 
         TiXmlElement *xmlCol = master->FirstChildElement("COLLECTIONS");
 
+        // Check if user has specified some options
         if (xmlCol)
         {
+            // Set the maxsize and default implementation type if provided
             const char *maxSize = xmlCol->Attribute("MAXSIZE");
             m_maxCollSize = (maxSize ? atoi(maxSize) : 0);
 
             const char *defaultImpl = xmlCol->Attribute("DEFAULT");
             m_defaultType = defaultType;
+
+            // If user has specified a default impl type, autotuning
+            // and set this default across all operators.
             if (defaultType == eNoImpType && defaultImpl)
             {
                 const std::string collinfo = string(defaultImpl);
                 m_autotune = boost::iequals(collinfo, "auto");
+
                 if (!m_autotune)
                 {
                     for(i = 1; i < Collections::SIZE_ImplementationType; ++i)
@@ -122,7 +153,7 @@ CollectionOptimisation::CollectionOptimisation(
                     }
 
                     ASSERTL0(i != Collections::SIZE_ImplementationType,
-                             "Unknown default collection scheme: "+collinfo);
+                         "Unknown default collection scheme: "+collinfo);
 
                     // Override default types
                     for (it2 = elTypes.begin(); it2 != elTypes.end(); ++it2)
@@ -137,6 +168,7 @@ CollectionOptimisation::CollectionOptimisation(
                 }
             }
 
+            // Now process operator-specific implementation selections
             TiXmlElement *elmt = xmlCol->FirstChildElement();
             while (elmt)
             {
@@ -323,7 +355,23 @@ OperatorImpMap CollectionOptimisation::SetWithTimings(
     CollectionVector coll;
     for(int imp = 1; imp < SIZE_ImplementationType; ++imp)
     {
-        OperatorImpMap impTypes = SetFixedImpType((ImplementationType) imp);
+        ImplementationType impType = (ImplementationType)imp;
+        OperatorImpMap impTypes;
+        for (int i = 0; i < SIZE_OperatorType; ++i)
+        {
+            OperatorType opType = (OperatorType)i;
+            OperatorKey opKey(pCollExp[0]->DetShapeType(), opType, impType,
+                              pCollExp[0]->IsNodalNonTensorialExp());
+
+            if (GetOperatorFactory().ModuleExists(opKey))
+            {
+                impTypes[opType] = impType;
+            }
+            else
+            {
+                cout << "Note: Implementation does not exist: " << opKey << endl;
+            }
+        }
 
         Collection collloc(pCollExp,impTypes);
         coll.push_back(collloc);
@@ -357,17 +405,24 @@ OperatorImpMap CollectionOptimisation::SetWithTimings(
         // call collection implementation in thorugh ExpList.
         for (int imp = 0; imp < coll.size(); ++imp)
         {
-            t.Start();
-            for(int n = 0; n < Ntest[i]; ++n)
+            if (coll[imp].HasOperator(OpType))
             {
-                coll[imp].ApplyOperator(OpType,
+                t.Start();
+                for(int n = 0; n < Ntest[i]; ++n)
+                {
+                    coll[imp].ApplyOperator(OpType,
                                       inarray,
                                       outarray1,
                                       outarray2,
                                       outarray3);
+                }
+                t.Stop();
+                timing[imp] = t.TimePerTest(Ntest[i]);
+            }
+            else
+            {
+                timing[imp] = 1000.0;
             }
-            t.Stop();
-            timing[imp] = t.TimePerTest(Ntest[i]);
         }
         // determine optimal implementation. Note +1 to
         // remove NoImplementationType flag
@@ -379,7 +434,14 @@ OperatorImpMap CollectionOptimisation::SetWithTimings(
                  << ImplementationTypeMap[minImp] << "\t (";
             for(int j = 0; j < coll.size(); ++j)
             {
-                cout << timing[j] ;
+                if (timing[j] > 999.0)
+                {
+                    cout << "-";
+                }
+                else
+                {
+                    cout << timing[j] ;
+                }
                 if(j != coll.size()-1)
                 {
                     cout <<", ";

library/Collections/Operator.cpp

@@ -88,10 +88,10 @@ bool operator< (OperatorKey const &p1, OperatorKey const &p2)
  */
 std::ostream &operator<<(std::ostream &os, OperatorKey const &p)
 {
-    os <<                       boost::get<0>(p)  << " "
-       << OperatorTypeMap      [boost::get<1>(p)] << " "
-       << ImplementationTypeMap[boost::get<2>(p)] << " "
-       << ImplementationTypeMap[boost::get<3>(p)];
+    os << LibUtilities::ShapeTypeMap[boost::get<0>(p)] << ", "
+       << OperatorTypeMap           [boost::get<1>(p)] << ", "
+       << ImplementationTypeMap     [boost::get<2>(p)] << ", "
+       << (boost::get<3>(p) ? "Nodal" : "Modal");
     return os;
 }