Multigrid Linear System Solver (#104)

This revision includes:
* Successive Over-Relaxation (SOR) method implementation
* Red-black Gauss-Seidel method implementation
* Multigrid (MG) method implementation (Issue #17)
* Multigrid Preconditioned Conjugate Gradient (MGPCG) method implementation (Issue #17)
* Python bindings for all the new APIs for the new solvers above

Author: doyubkim

include/jet/detail/cg-inl.h

@@ -34,8 +34,6 @@ void pcg(
     BlasType::set(0, q);
     BlasType::set(0, s);
 
-    M->build(A);
-
     // r = b - Ax
     BlasType::residual(A, *x, b, r);
 
@@ -90,7 +88,9 @@ void pcg(
     }
 
     *lastNumberOfIterations = iter;
-    *lastResidualNorm = std::sqrt(sigmaNew);
+
+    // std::fabs(sigmaNew) - Workaround for negative zero
+    *lastResidualNorm = std::sqrt(std::fabs(sigmaNew));
 }
 
 template <typename BlasType>

include/jet/detail/fdm_mg_linear_system2-inl.h

@@ -0,0 +1,52 @@
+// Copyright (c) 2017 Doyub Kim
+//
+// I am making my contributions/submissions to this project solely in my
+// personal capacity and am not conveying any rights to any intellectual
+// property of any third parties.
+
+#ifndef INCLUDE_JET_FDM_MG_LINEAR_SYSTEM2_INL_H_
+#define INCLUDE_JET_FDM_MG_LINEAR_SYSTEM2_INL_H_
+
+#include <jet/fdm_mg_linear_system2.h>
+
+namespace jet {
+
+template <typename T>
+void FdmMgUtils2::resizeArrayWithCoarsest(const Size2& coarsestResolution,
+                                          size_t numberOfLevels,
+                                          std::vector<Array2<T>>* levels) {
+    numberOfLevels = std::max(numberOfLevels, kOneSize);
+
+    levels->resize(numberOfLevels);
+
+    // Level 0 is the finest level, thus takes coarsestResolution ^
+    // numberOfLevels.
+    // Level numberOfLevels - 1 is the coarsest, taking coarsestResolution.
+    Size2 res = coarsestResolution;
+    for (size_t level = 0; level < numberOfLevels; ++level) {
+        (*levels)[numberOfLevels - level - 1].resize(res);
+        res.x = res.x << 1;
+        res.y = res.y << 1;
+    }
+}
+
+template <typename T>
+void FdmMgUtils2::resizeArrayWithFinest(const Size2& finestResolution,
+                                        size_t maxNumberOfLevels,
+                                        std::vector<Array2<T>>* levels) {
+    Size2 res = finestResolution;
+    size_t i = 1;
+    for (; i < maxNumberOfLevels; ++i) {
+        if (res.x % 2 == 0 && res.y % 2 == 0) {
+            res.x = res.x >> 1;
+            res.y = res.y >> 1;
+        } else {
+            break;
+        }
+    }
+    resizeArrayWithCoarsest(res, i, levels);
+}
+
+}  // namespace jet
+
+#endif  // INCLUDE_JET_FDM_MG_LINEAR_SYSTEM2_INL_H_

include/jet/detail/fdm_mg_linear_system3-inl.h

@@ -0,0 +1,54 @@
+// Copyright (c) 2017 Doyub Kim
+//
+// I am making my contributions/submissions to this project solely in my
+// personal capacity and am not conveying any rights to any intellectual
+// property of any third parties.
+
+#ifndef INCLUDE_JET_FDM_MG_LINEAR_SYSTEM3_INL_H_
+#define INCLUDE_JET_FDM_MG_LINEAR_SYSTEM3_INL_H_
+
+#include <jet/fdm_mg_linear_system3.h>
+
+namespace jet {
+
+template <typename T>
+void FdmMgUtils3::resizeArrayWithCoarsest(const Size3& coarsestResolution,
+                                          size_t numberOfLevels,
+                                          std::vector<Array3<T>>* levels) {
+    numberOfLevels = std::max(numberOfLevels, kOneSize);
+
+    levels->resize(numberOfLevels);
+
+    // Level 0 is the finest level, thus takes coarsestResolution ^
+    // numberOfLevels.
+    // Level numberOfLevels - 1 is the coarsest, taking coarsestResolution.
+    Size3 res = coarsestResolution;
+    for (size_t level = 0; level < numberOfLevels; ++level) {
+        (*levels)[numberOfLevels - level - 1].resize(res);
+        res.x = res.x << 1;
+        res.y = res.y << 1;
+        res.z = res.z << 1;
+    }
+}
+
+template <typename T>
+void FdmMgUtils3::resizeArrayWithFinest(const Size3& finestResolution,
+                                        size_t maxNumberOfLevels,
+                                        std::vector<Array3<T>>* levels) {
+    Size3 res = finestResolution;
+    size_t i = 1;
+    for (; i < maxNumberOfLevels; ++i) {
+        if (res.x % 2 == 0 && res.y % 2 == 0 && res.z % 2 == 0) {
+            res.x = res.x >> 1;
+            res.y = res.y >> 1;
+            res.z = res.z >> 1;
+        } else {
+            break;
+        }
+    }
+    resizeArrayWithCoarsest(res, i, levels);
+}
+
+}  // namespace jet
+
+#endif  // INCLUDE_JET_FDM_MG_LINEAR_SYSTEM3_INL_H_

include/jet/detail/math_utils-inl.h

@@ -38,7 +38,6 @@ inline T max3(T x, T y, T z) {
     return std::max(std::max(x, y), z);
 }
 
-//! Returns minimum among n-elements.
 template <typename T>
 inline T minn(const T* x, size_t n) {
     T m = x[0];
@@ -48,7 +47,6 @@ inline T minn(const T* x, size_t n) {
     return m;
 }
 
-//! Returns maximum among n-elements.
 template <typename T>
 inline T maxn(const T* x, size_t n) {
     T m = x[0];
@@ -68,7 +66,6 @@ inline T absmax(T x, T y) {
     return (x*x > y*y) ? x : y;
 }
 
-//! Returns absolute minimum among n-elements.
 template <typename T>
 inline T absminn(const T* x, size_t n) {
     T m = x[0];
@@ -78,7 +75,6 @@ inline T absminn(const T* x, size_t n) {
     return m;
 }
 
-//! Returns absolute maximum among n-elements.
 template <typename T>
 inline T absmaxn(const T* x, size_t n) {
     T m = x[0];
@@ -88,6 +84,34 @@ inline T absmaxn(const T* x, size_t n) {
     return m;
 }
 
+template <typename T>
+inline size_t argmin2(T x, T y) {
+    return (x < y) ? 0 : 1;
+}
+
+template <typename T>
+inline size_t argmax2(T x, T y) {
+    return (x > y) ? 0 : 1;
+}
+
+template <typename T>
+inline size_t argmin3(T x, T y, T z) {
+    if (x < y) {
+        return (x < z) ? 0 : 2;
+    } else {
+        return (y < z) ? 1 : 2;
+    }
+}
+
+template <typename T>
+inline size_t argmax3(T x, T y, T z) {
+    if (x > y) {
+        return (x > z) ? 0 : 2;
+    } else {
+        return (y > z) ? 1 : 2;
+    }
+}
+
 template <typename T>
 inline T square(T x) {
     return x * x;

include/jet/detail/mg-inl.h

@@ -23,7 +23,7 @@ MgResult mgVCycle(const MgMatrix<BlasType>& A, MgParameters<BlasType> params,
                      &((*x)[currentLevel]), &((*buffer)[currentLevel]));
 
     // 2) if currentLevel is the coarsest grid, goto 5)
-    if (currentLevel < params.maxNumberOfLevels - 1) {
+    if (currentLevel < A.levels.size() - 1) {
         auto r = buffer;
         BlasType::residual(A[currentLevel], (*x)[currentLevel],
                            (*b)[currentLevel], &(*r)[currentLevel]);
@@ -80,6 +80,16 @@ typename BlasType::MatrixType& MgMatrix<BlasType>::operator[](size_t i) {
     return levels[i];
 }
 
+template <typename BlasType>
+const typename BlasType::MatrixType& MgMatrix<BlasType>::finest() const {
+    return levels.front();
+}
+
+template <typename BlasType>
+typename BlasType::MatrixType& MgMatrix<BlasType>::finest() {
+    return levels.front();
+}
+
 template <typename BlasType>
 const typename BlasType::VectorType& MgVector<BlasType>::operator[](
     size_t i) const {
@@ -91,6 +101,16 @@ typename BlasType::VectorType& MgVector<BlasType>::operator[](size_t i) {
     return levels[i];
 }
 
+template <typename BlasType>
+const typename BlasType::VectorType& MgVector<BlasType>::finest() const {
+    return levels.front();
+}
+
+template <typename BlasType>
+typename BlasType::VectorType& MgVector<BlasType>::finest() {
+    return levels.front();
+}
+
 template <typename BlasType>
 MgResult mgVCycle(const MgMatrix<BlasType>& A, MgParameters<BlasType> params,
                   MgVector<BlasType>* x, MgVector<BlasType>* b,

include/jet/detail/parallel-inl.h

@@ -158,26 +158,80 @@ void parallelFor(IndexType start, IndexType end, const Function& func,
     }
 }
 
+template <typename IndexType, typename Function>
+void parallelRangeFor(IndexType start, IndexType end, const Function& func,
+                      ExecutionPolicy policy) {
+    if (start > end) {
+        return;
+    }
+
+    // Estimate number of threads in the pool
+    static const unsigned int numThreadsHint =
+        std::thread::hardware_concurrency();
+    const unsigned int numThreads =
+        (policy == ExecutionPolicy::kParallel)
+            ? (numThreadsHint == 0u ? 8u : numThreadsHint)
+            : 1;
+
+    // Size of a slice for the range functions
+    IndexType n = end - start + 1;
+    IndexType slice =
+        (IndexType)std::round(n / static_cast<double>(numThreads));
+    slice = std::max(slice, IndexType(1));
+
+    // Create pool and launch jobs
+    std::vector<std::thread> pool;
+    pool.reserve(numThreads);
+    IndexType i1 = start;
+    IndexType i2 = std::min(start + slice, end);
+    for (unsigned int i = 0; i + 1 < numThreads && i1 < end; ++i) {
+        pool.emplace_back(func, i1, i2);
+        i1 = i2;
+        i2 = std::min(i2 + slice, end);
+    }
+    if (i1 < end) {
+        pool.emplace_back(func, i1, end);
+    }
+
+    // Wait for jobs to finish
+    for (std::thread& t : pool) {
+        if (t.joinable()) {
+            t.join();
+        }
+    }
+}
+
 template <typename IndexType, typename Function>
 void parallelFor(IndexType beginIndexX, IndexType endIndexX,
                  IndexType beginIndexY, IndexType endIndexY,
                  const Function& function, ExecutionPolicy policy) {
     parallelFor(beginIndexY, endIndexY,
-                [&](size_t j) {
+                [&](IndexType j) {
                     for (IndexType i = beginIndexX; i < endIndexX; ++i) {
                         function(i, j);
                     }
                 },
                 policy);
 }
 
+template <typename IndexType, typename Function>
+void parallelRangeFor(IndexType beginIndexX, IndexType endIndexX,
+                      IndexType beginIndexY, IndexType endIndexY,
+                      const Function& function, ExecutionPolicy policy) {
+    parallelRangeFor(beginIndexY, endIndexY,
+                     [&](IndexType jBegin, IndexType jEnd) {
+                         function(beginIndexX, endIndexX, jBegin, jEnd);
+                     },
+                     policy);
+}
+
 template <typename IndexType, typename Function>
 void parallelFor(IndexType beginIndexX, IndexType endIndexX,
                  IndexType beginIndexY, IndexType endIndexY,
                  IndexType beginIndexZ, IndexType endIndexZ,
                  const Function& function, ExecutionPolicy policy) {
     parallelFor(beginIndexZ, endIndexZ,
-                [&](size_t k) {
+                [&](IndexType k) {
                     for (IndexType j = beginIndexY; j < endIndexY; ++j) {
                         for (IndexType i = beginIndexX; i < endIndexX; ++i) {
                             function(i, j, k);
@@ -187,6 +241,19 @@ void parallelFor(IndexType beginIndexX, IndexType endIndexX,
                 policy);
 }
 
+template <typename IndexType, typename Function>
+void parallelRangeFor(IndexType beginIndexX, IndexType endIndexX,
+                      IndexType beginIndexY, IndexType endIndexY,
+                      IndexType beginIndexZ, IndexType endIndexZ,
+                      const Function& function, ExecutionPolicy policy) {
+    parallelRangeFor(beginIndexZ, endIndexZ,
+                     [&](IndexType kBegin, IndexType kEnd) {
+                         function(beginIndexX, endIndexX, beginIndexY,
+                                  endIndexY, kBegin, kEnd);
+                     },
+                     policy);
+}
+
 template <typename IndexType, typename Value, typename Function,
           typename Reduce>
 Value parallelReduce(IndexType start, IndexType end, const Value& identity,

include/jet/fdm_gauss_seidel_solver2.h

@@ -16,10 +16,10 @@ namespace jet {
 class FdmGaussSeidelSolver2 final : public FdmLinearSystemSolver2 {
  public:
     //! Constructs the solver with given parameters.
-    FdmGaussSeidelSolver2(
-        unsigned int maxNumberOfIterations,
-        unsigned int residualCheckInterval,
-        double tolerance);
+    FdmGaussSeidelSolver2(unsigned int maxNumberOfIterations,
+                          unsigned int residualCheckInterval, double tolerance,
+                          double sorFactor = 1.0,
+                          bool useRedBlackOrdering = false);
 
     //! Solves the given linear system.
     bool solve(FdmLinearSystem2* system) override;
@@ -36,16 +36,30 @@ class FdmGaussSeidelSolver2 final : public FdmLinearSystemSolver2 {
     //! Returns the last residual after the Gauss-Seidel iterations.
     double lastResidual() const;
 
+    //! Returns the SOR (Successive Over Relaxation) factor.
+    double sorFactor() const;
+
+    //! Returns true if red-black ordering is enabled.
+    bool useRedBlackOrdering() const;
+
+    //! Performs single natural Gauss-Seidel relaxation step.
+    static void relax(const FdmMatrix2& A, const FdmVector2& b,
+                      double sorFactor, FdmVector2* x);
+
+    //! Performs single Red-Black Gauss-Seidel relaxation step.
+    static void relaxRedBlack(const FdmMatrix2& A, const FdmVector2& b,
+                              double sorFactor, FdmVector2* x);
+
  private:
     unsigned int _maxNumberOfIterations;
     unsigned int _lastNumberOfIterations;
     unsigned int _residualCheckInterval;
     double _tolerance;
     double _lastResidual;
+    double _sorFactor;
+    bool _useRedBlackOrdering;
 
     FdmVector2 _residual;
-
-    void relax(FdmLinearSystem2* system);
 };
 
 //! Shared pointer type for the FdmGaussSeidelSolver2.