text: // This file is part of Eigen, a lightweight C++ template library
  // for linear algebra.
  //
  // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
  //
  // This Source Code Form is subject to the terms of the Mozilla
  // Public License 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/.
  
  #ifndef EIGEN_PARALLELIZER_H
  #define EIGEN_PARALLELIZER_H
  
  namespace Eigen {
  
  namespace internal {
  
  /** \internal */
  inline void manage_multi_threading(Action action, int *v) {
      static EIGEN_UNUSED int m_maxThreads = -1;
  
      if (action == SetAction) {
          eigen_internal_assert(v != 0);
          m_maxThreads = *v;
      } else if (action == GetAction) {
          eigen_internal_assert(v != 0);
  #ifdef EIGEN_HAS_OPENMP
          if (m_maxThreads > 0)
              *v = m_maxThreads;
          else
              *v = omp_get_max_threads();
  #else
          *v = 1;
  #endif
      } else {
          eigen_internal_assert(false);
      }
  }
  
  } // namespace internal
  
  /** Must be call first when calling Eigen from multiple threads */
  inline void initParallel() {
      int nbt;
      internal::manage_multi_threading(GetAction, &nbt);
      std::ptrdiff_t l1, l2;
      internal::manage_caching_sizes(GetAction, &l1, &l2);
  }
  
  /** \returns the max number of threads reserved for Eigen
   * \sa setNbThreads */
  inline int nbThreads() {
      int ret;
      internal::manage_multi_threading(GetAction, &ret);
      return ret;
  }
  
  /** Sets the max number of threads reserved for Eigen
   * \sa nbThreads */
  inline void setNbThreads(int v) {
      internal::manage_multi_threading(SetAction, &v);
  }
  
  namespace internal {
  
  template <typename Index> struct GemmParallelInfo {
      GemmParallelInfo() : sync(-1), users(0), rhs_start(0), rhs_length(0) {}
  
      int volatile sync;
      int volatile users;
  
      Index rhs_start;
      Index rhs_length;
  };
  
  template <bool Condition, typename Functor, typename Index>
  void parallelize_gemm(const Functor &func, Index rows, Index cols,
                        bool transpose) {
      // TODO when EIGEN_USE_BLAS is defined,
      // we should still enable OMP for other scalar types
  #if !(defined(EIGEN_HAS_OPENMP)) || defined(EIGEN_USE_BLAS)
      // FIXME the transpose variable is only needed to properly split
      // the matrix product when multithreading is enabled. This is a temporary
      // fix to support row-major destination matrices. This whole
      // parallelizer mechanism has to be redisigned anyway.
      EIGEN_UNUSED_VARIABLE(transpose);
      func(0, rows, 0, cols);
  #else
  
      // Dynamically check whether we should enable or disable OpenMP.
      // The conditions are:
      // - the max number of threads we can create is greater than 1
      // - we are not already in a parallel code
      // - the sizes are large enough
  
      // 1- are we already in a parallel session?
      // FIXME omp_get_num_threads()>1 only works for openmp, what if the user
      // does not use openmp?
      if ((!Condition) || (omp_get_num_threads() > 1))
          return func(0, rows, 0, cols);
  
      Index size = transpose ? cols : rows;
  
      // 2- compute the maximal number of threads from the size of the product:
      // FIXME this has to be fine tuned
      Index max_threads = std::max<Index>(1, size / 32);
  
      // 3 - compute the number of threads we are going to use
      Index threads = std::min<Index>(nbThreads(), max_threads);
  
      if (threads == 1)
          return func(0, rows, 0, cols);
  
      Eigen::initParallel();
      func.initParallelSession();
  
      if (transpose)
          std::swap(rows, cols);
  
      Index blockCols = (cols / threads) & ~Index(0x3);
      Index blockRows = (rows / threads) & ~Index(0x7);
  
      GemmParallelInfo<Index> *info = new GemmParallelInfo<Index>[threads];
  
  #pragma omp parallel for schedule(static, 1) num_threads(threads)
      for (Index i = 0; i < threads; ++i) {
          Index r0 = i * blockRows;
          Index actualBlockRows = (i + 1 == threads) ? rows - r0 : blockRows;
  
          Index c0 = i * blockCols;
          Index actualBlockCols = (i + 1 == threads) ? cols - c0 : blockCols;
  
          info[i].rhs_start = c0;
          info[i].rhs_length = actualBlockCols;
  
          if (transpose)
              func(0, cols, r0, actualBlockRows, info);
          else
              func(r0, actualBlockRows, 0, cols, info);
      }
  
      delete[] info;
  #endif
  }
  
  } // end namespace internal
  
  } // end namespace Eigen
  
  #endif // EIGEN_PARALLELIZER_H
  