Remove omp pragmas (#2242) [ci skip]

* remove omp pragmas

* continue -> return

include/igl/bfs_orient.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "bfs_orient.h"
 #include "orientable_patches.h"
+#include "parallel_for.h"
 #include <Eigen/Sparse>
 #include <queue>
 
@@ -35,8 +36,7 @@ IGL_INLINE void igl::bfs_orient(
     FF = F;
   }
   // loop over patches
-#pragma omp parallel for
-  for(int c = 0;c<num_cc;c++)
+  parallel_for(num_cc,[&](const int c)
   {
     queue<typename DerivedF::Scalar> Q;
     // find first member of patch c
@@ -89,7 +89,7 @@ IGL_INLINE void igl::bfs_orient(
         }
       }
     }
-  }
+  },1000);
 
   // make sure flip is OK if &FF = &F
 }

include/igl/ceil.cpp

@@ -14,15 +14,6 @@ IGL_INLINE void igl::ceil(
   Eigen::PlainObjectBase<DerivedY>& Y)
 {
   using namespace std;
-  //Y = DerivedY::Zero(m,n);
-//#pragma omp parallel for
-  //for(int i = 0;i<m;i++)
-  //{
-  //  for(int j = 0;j<n;j++)
-  //  {
-  //    Y(i,j) = std::ceil(X(i,j));
-  //  }
-  //}
   typedef typename DerivedX::Scalar Scalar;
   Y = X.unaryExpr([](const Scalar &x)->Scalar{return std::ceil(x);}).template cast<typename DerivedY::Scalar >();
 }

include/igl/cumprod.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cumprod.h"
+#include "parallel_for.h"
 #include <numeric>
 #include <iostream>
 
@@ -26,8 +27,7 @@ IGL_INLINE void igl::cumprod(
   // (Optimizations assume ColMajor order)
   if(dim == 1)
   {
-#pragma omp parallel for
-    for(int o = 0;o<num_outer;o++)
+    parallel_for(num_outer,[&](const int o)
     {
       typename DerivedX::Scalar prod = 1;
       for(int i = 0;i<num_inner;i++)
@@ -35,15 +35,14 @@ IGL_INLINE void igl::cumprod(
         prod *= X(i,o);
         Y(i,o) = prod;
       }
-    }
+    },1000);
   }else
   {
     for(int i = 0;i<num_inner;i++)
     {
-      // Notice that it is *not* OK to put this above the inner loop
+      // Notice that it is *not* OK to put parallel_for this above the inner loop
       // Though here it doesn't seem to pay off...
-//#pragma omp parallel for
-      for(int o = 0;o<num_outer;o++)
+      parallel_for(num_outer,[&](const int o)
       {
         if(i == 0)
         {
@@ -52,7 +51,7 @@ IGL_INLINE void igl::cumprod(
         {
           Y(o,i) = Y(o,i-1) * X(o,i);
         }
-      }
+      },1000);
     }
   }
 }

include/igl/cumsum.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cumsum.h"
+#include "parallel_for.h"
 #include <numeric>
 #include <iostream>
 
@@ -42,8 +43,7 @@ IGL_INLINE void igl::cumsum(
     {
       Y.row(0).setConstant(0);
     }
-#pragma omp parallel for
-    for(Eigen::Index o = 0;o<num_outer;o++)
+    parallel_for(num_outer,[&](const int o)
     {
       typename DerivedX::Scalar sum = 0;
       for(Eigen::Index i = 0;i<num_inner;i++)
@@ -52,7 +52,7 @@ IGL_INLINE void igl::cumsum(
         const Eigen::Index yi = zero_prefix?i+1:i;
         Y(yi,o) = sum;
       }
-    }
+    },1000);
   }else
   {
     if(zero_prefix)
@@ -62,10 +62,7 @@ IGL_INLINE void igl::cumsum(
     for(Eigen::Index i = 0;i<num_inner;i++)
     {
       const Eigen::Index yi = zero_prefix?i+1:i;
-      // Notice that it is *not* OK to put this above the inner loop
-      // Though here it doesn't seem to pay off...
-//#pragma omp parallel for
-      for(Eigen::Index o = 0;o<num_outer;o++)
+      parallel_for(num_outer,[&](const int o)
       {
         if(i == 0)
         {
@@ -74,7 +71,7 @@ IGL_INLINE void igl::cumsum(
         {
           Y(o,yi) = Y(o,yi-1) + X(o,i);
         }
-      }
+      },1000);
     }
   }
 }
@@ -98,4 +95,4 @@ template void igl::cumsum<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>
 template void igl::cumsum<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>, class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> &);
 template void igl::cumsum<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> > const&, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> >&);
 #endif
-#endif
+#endif

include/igl/dqs.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "dqs.h"
+#include "parallel_for.h"
 #include <Eigen/Geometry>
 template <
   typename DerivedV,
@@ -41,8 +42,7 @@ IGL_INLINE void igl::dqs(
 
   // Loop over vertices
   const int nv = V.rows();
-#pragma omp parallel for if (nv>10000)
-  for(int i = 0;i<nv;i++)
+  parallel_for(nv,[&](const int i)
   {
     Q b0(0,0,0,0);
     Q be(0,0,0,0);
@@ -64,7 +64,7 @@ IGL_INLINE void igl::dqs(
     typename Q::Scalar a0 = c0.w();
     typename Q::Scalar ae = ce.w();
     U.row(i) =  v + 2*d0.cross(d0.cross(v) + a0*v) + 2*(a0*de - ae*d0 + d0.cross(de));
-  }
+  },1000);
 
 }
 

include/igl/embree/bone_visible.cpp

@@ -9,6 +9,7 @@
 #include "../project_to_line.h"
 #include "../EPS.h"
 #include "../Hit.h"
+#include "../parallel_for.h"
 #include "../Timer.h"
 #include <iostream>
 
@@ -52,9 +53,8 @@ IGL_INLINE void igl::embree::bone_visible(
   flag.resize(V.rows());
   const double sd_norm = (s-d).norm();
   // Embree seems to be parallel when constructing but not when tracing rays
-#pragma omp parallel for
   // loop over mesh vertices
-  for(int v = 0;v<V.rows();v++)
+  parallel_for(V.rows(),[&](const int v)
   {
     const Vector3d Vv = V.row(v);
     // Project vertex v onto line segment sd
@@ -135,7 +135,7 @@ IGL_INLINE void igl::embree::bone_visible(
       // no hit so vectex v is visible
       flag(v) = true;
     }
-  }
+  },10000);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/embree/reorient_facets_raycast.cpp

@@ -154,7 +154,7 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
   vector<pair<int  , int  >> C_vote_parity(num_cc, make_pair(0, 0));        // sum of parity count for each ray
 
   if (is_verbose) cout << "shooting rays... ";
-#pragma omp parallel for
+// #pragma omp parallel for
   for (int i = 0; i < (int)ray_face.size(); ++i)
   {
     int      f = ray_face[i];
@@ -173,27 +173,27 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
     if (!hits_back .empty() && hits_back [0].id == f) hits_back .erase(hits_back .begin());
 
     if (use_parity) {
-#pragma omp atomic
+// #pragma omp atomic
       C_vote_parity[c].first  += hits_front.size() % 2;
-#pragma omp atomic
+// #pragma omp atomic
       C_vote_parity[c].second += hits_back .size() % 2;
 
     } else {
       if (hits_front.empty())
       {
-#pragma omp atomic
+// #pragma omp atomic
         C_vote_infinity[c].first++;
       } else {
-#pragma omp atomic
+// #pragma omp atomic
         C_vote_distance[c].first += hits_front[0].t;
       }
 
       if (hits_back.empty())
       {
-#pragma omp atomic
+// #pragma omp atomic
         C_vote_infinity[c].second++;
       } else {
-#pragma omp atomic
+// #pragma omp atomic
         C_vote_distance[c].second += hits_back[0].t;
       }
     }

include/igl/floor.cpp

@@ -15,15 +15,6 @@ IGL_INLINE void igl::floor(
   Eigen::PlainObjectBase<DerivedY>& Y)
 {
   using namespace std;
-  //Y = DerivedY::Zero(m,n);
-//#pragma omp parallel for
-  //for(int i = 0;i<m;i++)
-  //{
-  //  for(int j = 0;j<n;j++)
-  //  {
-  //    Y(i,j) = std::floor(X(i,j));
-  //  }
-  //}
   typedef typename DerivedX::Scalar Scalar;
   Y = X.unaryExpr([](const Scalar &x)->Scalar{return std::floor(x);}).template cast<typename DerivedY::Scalar >();
 }

include/igl/gaussian_curvature.cpp

@@ -32,18 +32,11 @@ IGL_INLINE void igl::gaussian_curvature(
   assert(K.cols() == 1);
   const int Frows = F.rows();
   //K_G(x_i) = (2π - ∑θj)
-//#ifndef IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE
-//#  define IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE 1000
-//#endif
-//#pragma omp parallel for if (Frows>IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE)
   for(int f = 0;f<Frows;f++)
   {
     // throw normal at each corner
     for(int j = 0; j < 3;j++)
     {
-      // Q: Does this need to be critical?
-      // H: I think so, sadly. Maybe there's a way to use reduction
-//#pragma omp critical
       K(F(f,j),0) -=  A(f,j);
     }
   }

include/igl/histc.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "histc.h"
+#include "parallel_for.h"
 #include <cassert>
 #include <iostream>
 
@@ -22,12 +23,10 @@ IGL_INLINE void igl::histc(
   assert(m == B.size());
   N.resize(n,1);
   N.setConstant(0);
-#pragma omp parallel for
   for(int j = 0;j<m;j++)
   {
     if(B(j) >= 0)
     {
-#pragma omp atomic
       N(int(B(j)))++;
     }
   }
@@ -46,15 +45,14 @@ IGL_INLINE void igl::histc(
       E.topLeftCorner(E.size()-1,1)).maxCoeff() >= 0 && 
     "E should be monotonically increasing");
   B.resize(m,1);
-#pragma omp parallel for
-  for(int j = 0;j<m;j++)
+  parallel_for(m,[&](const int j)
   {
     const double x = X(j);
     // Boring one-offs
     if(x < E(0) || x > E(E.size()-1))
     {
       B(j) = -1;
-      continue;
+      return;
     }
     // Find x in E
     int l = 0;
@@ -81,7 +79,7 @@ IGL_INLINE void igl::histc(
       k = l;
     }
     B(j) = k;
-  }
+  },1000);
 }
 
 template <typename DerivedE>

include/igl/in_element.cpp

@@ -6,7 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "in_element.h"
-
+#include "parallel_for.h"
 template <typename DerivedV, typename DerivedQ, int DIM>
 IGL_INLINE void igl::in_element(
   const Eigen::MatrixBase<DerivedV> & V,
@@ -19,8 +19,7 @@ IGL_INLINE void igl::in_element(
   using namespace Eigen;
   const int Qr = Q.rows();
   I.setConstant(Qr,1,-1);
-#pragma omp parallel for if (Qr>10000)
-  for(int e = 0;e<Qr;e++)
+  parallel_for(Qr,[&](const int e)
   {
     // find all
     const auto R = aabb.find(V,Ele,Q.row(e).eval(),true);
@@ -28,7 +27,7 @@ IGL_INLINE void igl::in_element(
     {
       I(e) = R[0];
     }
-  }
+  },10000);
 }
 
 template <typename DerivedV, typename DerivedQ, int DIM, typename Scalar>
@@ -44,14 +43,14 @@ IGL_INLINE void igl::in_element(
   const int Qr = Q.rows();
   std::vector<Triplet<Scalar> > IJV;
   IJV.reserve(Qr);
-#pragma omp parallel for if (Qr>10000)
+// #pragma omp parallel for if (Qr>10000)
   for(int e = 0;e<Qr;e++)
   {
     // find all
     const auto R = aabb.find(V,Ele,Q.row(e).eval(),false);
     for(const auto r : R)
     {
-#pragma omp critical
+// #pragma omp critical
       IJV.push_back(Triplet<Scalar>(e,r,1));
     }
   }

include/igl/per_face_normals.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "per_face_normals.h"
+#include "parallel_for.h"
 #include <Eigen/Geometry>
 
 #define SQRT_ONE_OVER_THREE 0.57735026918962573
@@ -19,8 +20,7 @@ IGL_INLINE void igl::per_face_normals(
   N.resize(F.rows(),3);
   // loop over faces
   int Frows = F.rows();
-#pragma omp parallel for if (Frows>10000)
-  for(int i = 0; i < Frows;i++)
+  parallel_for(Frows,[&](const int i)
   {
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v2 = V.row(F(i,2)) - V.row(F(i,0));
@@ -33,7 +33,7 @@ IGL_INLINE void igl::per_face_normals(
     {
       N.row(i) /= r;
     }
-  }
+  },10000);
 }
 
 template <typename DerivedV, typename DerivedF, typename DerivedN>
@@ -59,7 +59,7 @@ IGL_INLINE void igl::per_face_normals_stable(
 
   N.resize(F.rows(),3);
   // Grad all points
-  for(size_t f = 0;f<m;f++)
+  parallel_for(m,[&](const int f)
   {
     const RowVectorV3 p0 = V.row(F(f,0));
     const RowVectorV3 p1 = V.row(F(f,1));
@@ -97,7 +97,7 @@ IGL_INLINE void igl::per_face_normals_stable(
     }
     // sum better not be sure, or else NaN
     N.row(f) /= N.row(f).norm();
-  }
+  },10000);
 
 }
 

include/igl/project_to_line.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "project_to_line.h"
+#include "parallel_for.h"
 #include <cassert>
 #include <Eigen/Core>
 
@@ -40,8 +41,7 @@ IGL_INLINE void igl::project_to_line(
   t.resize(np,1);
   sqrD.resize(np,1);
   // loop over points
-#pragma omp parallel for if (np>10000)
-  for(int i = 0;i<np;i++)
+  parallel_for(np,[&](const int i)
   {
     const typename DerivedP::ConstRowXpr Pi = P.row(i);
     // vector from point i to source
@@ -50,7 +50,7 @@ IGL_INLINE void igl::project_to_line(
     // P projected onto line
     const DerivedD projP = (1-t(i))*S + t(i)*D;
     sqrD(i) = (Pi-projP).squaredNorm();
-  }
+  },10000);
 }
 
 template <typename Scalar>

include/igl/project_to_line_segment.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "project_to_line_segment.h"
 #include "project_to_line.h"
+#include "parallel_for.h"
 #include <Eigen/Core>
 
 template <
@@ -25,8 +26,7 @@ IGL_INLINE void igl::project_to_line_segment(
   project_to_line(P,S,D,t,sqrD);
   const int np = P.rows();
   // loop over points and fix those that projected beyond endpoints
-#pragma omp parallel for if (np>10000)
-  for(int p = 0;p<np;p++)
+  parallel_for(np,[&](const int p)
   {
     const DerivedP Pp = P.row(p);
     if(t(p)<0)
@@ -38,7 +38,7 @@ IGL_INLINE void igl::project_to_line_segment(
       sqrD(p) = (Pp-D).squaredNorm();
       t(p) = 1;
     }
-  }
+  },10000);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/signed_distance.cpp

@@ -316,8 +316,7 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   N.resize(np,3);
   C.resize(np,3);
   typedef typename AABB<DerivedV,3>::RowVectorDIMS RowVector3S;
-# pragma omp parallel for if(np>1000)
-  for(std::ptrdiff_t p = 0;p<np;p++)
+  parallel_for(np,[&](const int p)
   {
     typename DerivedV::Scalar s,sqrd;
     RowVector3S n,c;
@@ -328,7 +327,7 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
     I(p) = i;
     N.row(p) = n;
     C.row(p) = c;
-  }
+  },1000);
 }
 
 template <