Feature 1926 gridstat openmp

met/configure.ac

@@ -1053,6 +1053,7 @@ AC_SUBST([PYTHON_LIBS])
 # Check for the MET_DEVELOPMENT environment variable
 
 AM_CONDITIONAL([ENABLE_DEVELOPMENT], [test -n "$MET_DEVELOPMENT"])
+AM_CONDITIONAL([ENABLE_OPENMP],      [test -n "$MET_OPENMP"])
 
 # Add MET_BASE to the CPPFLAGS
 
@@ -1067,9 +1068,17 @@ AC_SUBST(FC_LIBS, [-lgfortran])
 # For Intel compilers, "-g" slows down runtimes considerably (MET #1778).
 # For development, retain the "-g" option. Otherwise, discard it.
 
-AM_COND_IF([ENABLE_DEVELOPMENT], [: ${CFLAGS="-g -O2"}], [: ${CFLAGS="-O2"}])
-AM_COND_IF([ENABLE_DEVELOPMENT], [: ${CXXFLAGS="-g"}  ], [: ${CXXFLAGS=""}] )
-AM_COND_IF([ENABLE_DEVELOPMENT], [: ${FFLAGS="-g -O2"}], [: ${FFLAGS="-O2"}])
+AM_COND_IF([ENABLE_DEVELOPMENT], [: ${CFLAGS_OPT="-g -O2"}], [: ${CFLAGS_OPT="-O2"}])
+AM_COND_IF([ENABLE_DEVELOPMENT], [: ${CXXFLAGS_OPT="-g"}  ], [: ${CXXFLAGS_OPT="-O3"}] )
+AM_COND_IF([ENABLE_DEVELOPMENT], [: ${FFLAGS_OPT="-g -O2"}], [: ${FFLAGS_OPT="-O2"}])
+
+AM_COND_IF([ENABLE_OPENMP], [: ${CFLAGS_OMP="-fopenmp"}],   [: ${CFLAGS_OMP=""}])
+AM_COND_IF([ENABLE_OPENMP], [: ${CXXFLAGS_OMP="-fopenmp"}], [: ${CXXFLAGS_OMP=""}])
+AM_COND_IF([ENABLE_OPENMP], [: ${FFLAGS_OMP="-fopenmp"}],   [: ${FFLAGS_OMP=""}])
+
+CFLAGS="$CFLAGS_OPT $CFLAGS_OMP"
+CXXFLAGS="$CXXFLAGS_OPT $CXXFLAGS_OMP"
+FFLAGS="$FFLAGS_OPT $FFLAGS_OMP"
 
 # Checks for programs.
 

met/src/basic/vx_util/data_plane_util.cc

@@ -16,6 +16,8 @@ using namespace std;
 #include <string.h>
 #include <unistd.h>
 
+#include "omp.h"
+
 #include "data_plane_util.h"
 #include "interp_util.h"
 #include "two_to_one.h"
@@ -247,83 +249,92 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp,
       }
    }
 
-   // Build the grid template
-   GridTemplateFactory gtf;
-   GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon);
-
-   mlog << Debug(3)
-        << "Computing fractional coverage field using the "
-        << t.get_str() << " threshold and the "
-        << interpmthd_to_string(InterpMthd_Nbrhd) << "(" << gt->size()
-        << ") " << gt->getClassName() << " interpolation method.\n";
-
-   // Initialize the fractional coverage field
-   frac_dp = dp;
-   frac_dp.set_constant(bad_data_double);
-
-   // Compute the fractional coverage meeting the threshold criteria
-   for(x=0; x<dp.nx(); x++) {
-      for(y=0; y<dp.ny(); y++) {
-
-         // For a new column, reset the grid template and counts.
-         if(y == 0) {
-
-            // Initialize counts
-            n_vld = n_thr = 0;
-
-            // Sum all the points
-            for(gp  = gt->getFirstInGrid(x, y, dp.nx(), dp.ny());
-                gp != NULL;
-                gp  = gt->getNextInGrid()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld++;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
-            }
-         }
-         // Subtract off the bottom edge, shift up, and add the top.
-         else {
-
-            // Subtract points from the the bottom edge
-            for(gp  = gt->getFirstInBotEdge();
-                gp != NULL;
-                gp  = gt->getNextInBotEdge()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld--;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr--;
-            }
-
-            // Increment Y
-            gt->incBaseY(1);
-
-            // Add points from the the top edge
-            for(gp  = gt->getFirstInTopEdge();
-                gp != NULL;
-                gp  = gt->getNextInTopEdge()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld++;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
-            }
-         }
-
-         // Check for enough valid data and compute fractional coverage
-         if((double)(n_vld)/gt->size() >= vld_t && n_vld != 0) {
-            frac_dp.set((double) n_thr/n_vld, x, y);
-         }
-
-      } // end for y
-
-      // Increment X
-      if(x < (dp.nx() - 1)) gt->incBaseX(1);
-
-   } // end for x
-
-   delete gt;
+#pragma omp parallel default(shared) private(x,y,n_vld,n_thr,gp,v)
+   {
+
+     // Build the grid template
+     GridTemplateFactory gtf;
+     GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon);
+
+#pragma omp single
+     {
+       mlog << Debug(3)
+            << "Computing fractional coverage field using the "
+            << t.get_str() << " threshold and the "
+            << interpmthd_to_string(InterpMthd_Nbrhd) << "(" << gt->size()
+            << ") " << gt->getClassName() << " interpolation method.\n";
+
+       // Initialize the fractional coverage field
+       frac_dp = dp;
+       frac_dp.set_constant(bad_data_double);
+     }
+
+     // Compute the fractional coverage meeting the threshold criteria
+#pragma omp for schedule (static)
+     for(x=0; x<dp.nx(); x++) {
+        for(y=0; y<dp.ny(); y++) {
+
+           // For a new column, reset the grid template and counts.
+           if(y == 0) {
+
+              // Initialize counts
+              n_vld = n_thr = 0;
+
+              // Sum all the points
+              for(gp  = gt->getFirstInGrid(x, y, dp.nx(), dp.ny());
+                  gp != NULL;
+                  gp  = gt->getNextInGrid()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld++;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
+              }
+           }
+           // Subtract off the bottom edge, shift up, and add the top.
+           else {
+
+              // Subtract points from the the bottom edge
+              for(gp  = gt->getFirstInBotEdge();
+                  gp != NULL;
+                  gp  = gt->getNextInBotEdge()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld--;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr--;
+              }
+
+              // Increment Y
+              gt->incBaseY(1);
+
+              // Add points from the the top edge
+              for(gp  = gt->getFirstInTopEdge();
+                  gp != NULL;
+                  gp  = gt->getNextInTopEdge()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld++;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
+              }
+           }
+
+           // Check for enough valid data and compute fractional coverage
+           if((double)(n_vld)/gt->size() >= vld_t && n_vld != 0) {
+              frac_dp.set((double) n_thr/n_vld, x, y);
+           }
+
+        } // end for y
+
+        // Increment X
+        if(x < (dp.nx() - 1)) gt->incBaseX(1);
+
+     } // end for x
+
+     delete gt;
+
+   } // End of omp parallel
 
    return;
 }

met/src/tools/core/grid_stat/grid_stat.cc

@@ -127,6 +127,8 @@ using namespace std;
 #include <sys/types.h>
 #include <unistd.h>
 
+#include "omp.h"
+
 #include "grid_stat.h"
 
 #include "vx_statistics.h"
@@ -185,6 +187,15 @@ static bool read_data_plane(VarInfo* info, DataPlane& dp, Met2dDataFile* mtddf,
 
 int main(int argc, char *argv[]) {
 
+#ifdef _OPENMP
+#pragma omp parallel
+#pragma omp single
+  {
+    // Report number of threads if compiled with OpenMP
+    std::cout << "Running on " << omp_get_num_threads() << " threads.\n";
+  }
+#endif
+
    // Set handler to be called for memory allocation error
    set_new_handler(oom);