Merge pull request #215 from intel/develop

Develop

include/yask_common_api.hpp

@@ -33,6 +33,7 @@ IN THE SOFTWARE.
 #define YASK_COMMON_API
 
 #include <string>
+#include <vector>
 #include <iostream>
 #include <ostream>
 #include <memory>
@@ -204,6 +205,55 @@ namespace yask {
         virtual ~yask_null_output() {}
     };
 
+    /// Create finite-difference (FD) coefficients for the standard center form.
+    /**
+       Find FD coefficients with `radius` sample points to both the left and right
+       of the center sample and evaluation point on a uniformly-spaced grid. 
+       The FD has `radius * 2`-order accuracy.
+       @returns `radius * 2 + 1` FD coefficients.
+    */
+    std::vector<double>
+    get_center_fd_coefficients(int derivative_order
+                               /**< [in] `1` for 1st derivative, `2` for 2nd, etc. */,
+                               int radius
+                               /**< [in] Number of points to either side of the center point. */ );
+
+    /// Create finite-difference (FD) coefficients for the standard forward form.
+    /**
+       Find FD coefficients with `accuracy_order` sample points to the right
+       of the center sample and evaluation point on a uniformly-spaced grid. 
+       @returns `accuracy_order + 1` FD coefficients.
+    */
+    std::vector<double>
+    get_forward_fd_coefficients(int derivative_order
+                                /**< [in] `1` for 1st derivative, `2` for 2nd, etc. */,
+                                int accuracy_order
+                                /**< [in] Number of points to the right of the center point. */ );
+
+    /// Create finite-difference (FD) coefficients for the standard backward form.
+    /**
+       Find FD coefficients with `accuracy_order` sample points to the left
+       of the center sample and evaluation point on a uniformly-spaced grid. 
+       @returns `accuracy_order + 1` FD coefficients.
+    */
+    std::vector<double>
+    get_backward_fd_coefficients(int derivative_order
+                                 /**< [in] `1` for 1st derivative, `2` for 2nd, etc. */,
+                                 int accuracy_order
+                                 /**< [in] Number of points to the left of the center point. */ );
+
+    /// Create finite-difference (FD) coefficients at arbitrary evaluation and sample points.
+    /**
+       @returns `sample_points` FD coefficients.
+    */
+    std::vector<double>
+    get_arbitrary_fd_coefficients(int derivative_order
+                                  /**< [in] `1` for 1st derivative, `2` for 2nd, etc. */,
+                                  double eval_point
+                                  /**< [in] Location of evaluation point. */,
+                                  const std::vector<double> sample_points
+                                  /**< [in] Locations of sampled points. Must have at least 2. */ );
+
     /** @}*/
 
 } // namespace yask.

src/common/common.mk

@@ -55,7 +55,8 @@ endif
 
 # Common source.
 COMM_DIR	:=	$(SRC_DIR)/common
-COMM_SRC_NAMES	:=	output common_utils tuple combo
+COMM_SRC_NAMES	:=	output common_utils tuple combo fd_coeff fd_coeff2
+COEFF_DIR	:=	$(SRC_DIR)/contrib/coefficients
 
 # YASK stencil compiler.
 # This is here because both the compiler and kernel

src/common/common_utils.cpp

@@ -43,7 +43,7 @@ namespace yask {
     // for numbers above 9 (at least up to 99).
 
     // Format: "major.minor.patch".
-    const string version = "2.20.00";
+    const string version = "2.21.00";
 
     string yask_get_version_string() {
         return version;

src/common/fd_coeff2.cpp

@@ -0,0 +1,75 @@
+/*****************************************************************************
+
+YASK: Yet Another Stencil Kernel
+Copyright (c) 2014-2019, Intel Corporation
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to
+deal in the Software without restriction, including without limitation the
+rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+sell copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+* The above copyright notice and this permission notice shall be included in
+  all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+IN THE SOFTWARE.
+
+*****************************************************************************/
+
+///////// FD coefficients implementation.
+
+#include "yask_common_api.hpp"
+#include "fd_coeff.hpp"
+#include "common_utils.hpp"
+
+using namespace std;
+
+namespace yask {
+
+    // C++-style interface.
+    vector<double> get_arbitrary_fd_coefficients(int derivative_order, double eval_point, const vector<double> sample_points) {
+        if (derivative_order < 1)
+            THROW_YASK_EXCEPTION("Error: get_fd_coefficients() called with derivative-order less than 1");
+        int n = int(sample_points.size());
+        if (n < 2)
+            THROW_YASK_EXCEPTION("Error: get_fd_coefficients() called with fewer than 2 sample points");
+        vector<double> coeffs(n);
+        fd_coeff(&coeffs[0], double(eval_point), derivative_order, &sample_points[0], n);
+        return coeffs;
+    }
+
+    // Common FD forms for uniform grid spacing.
+    vector<double> get_center_fd_coefficients(int derivative_order, int radius) {
+        if (radius < 1)
+            THROW_YASK_EXCEPTION("get_center_fd_coefficients() called with less than radius 1");
+        vector<double> pts;
+        for (int i = -radius; i <= radius; i++)
+            pts.push_back(i);
+        assert(sizeof(pts) == size_t(radius * 2 + 1));
+        return get_arbitrary_fd_coefficients(derivative_order, 0, pts);
+    }
+    vector<double> get_forward_fd_coefficients(int derivative_order, int accuracy_order) {
+        if (accuracy_order < 1)
+            THROW_YASK_EXCEPTION("get_forward_fd_coefficients() called with less than order-of-accuracy 1");
+        vector<double> pts;
+        for (int i = 0; i <= accuracy_order; i++)
+            pts.push_back(i);
+        return get_arbitrary_fd_coefficients(derivative_order, 0, pts);
+    }
+    vector<double> get_backward_fd_coefficients(int derivative_order, int accuracy_order) {
+        if (accuracy_order < 1)
+            THROW_YASK_EXCEPTION("get_backward_fd_coefficients() called with less than order-of-accuracy 1");
+        vector<double> pts;
+        for (int i = -accuracy_order; i <= 0; i++)
+            pts.push_back(i);
+        return get_arbitrary_fd_coefficients(derivative_order, 0, pts);
+    }
+
+} // yask namespace.

src/compiler/Makefile

@@ -50,7 +50,7 @@ YC_SRC_NAMES	:=	Expr ExprUtils Grid Eqs Print Vec Cpp CppIntrin YaskKernel Soln
 YC_STENCIL_NAMES:=	$(notdir $(patsubst %.cpp,%,$(wildcard $(YC_STENCIL_DIR)/*.cpp)))
 YC_OBJS		:=	$(addprefix $(YC_OBJ_DIR)/,$(addsuffix .o,$(YC_SRC_NAMES) $(COMM_SRC_NAMES)))
 YC_STENCIL_OBJS	:=	$(addprefix $(YC_OBJ_DIR)/,$(addsuffix .o,$(YC_STENCIL_NAMES)))
-YC_INC_DIRS	:=	$(INC_DIR) $(YC_LIB_SRC_DIR) $(COMM_DIR)
+YC_INC_DIRS	:=	$(INC_DIR) $(YC_LIB_SRC_DIR) $(COMM_DIR) $(COEFF_DIR)
 YC_INC_GLOB	:=	$(wildcard $(addsuffix /*.hpp,$(YC_INC_DIRS)))
 YC_STENCIL_INC_GLOB :=	$(wildcard $(YC_STENCIL_DIR)/*.hpp $(YC_STENCIL_DIR)/*/*.hpp)
 
@@ -87,6 +87,11 @@ $(YC_OBJ_DIR)/%.o: $(COMM_DIR)/%.cpp $(YC_INC_GLOB)
 	$(CXX_PREFIX) $(YC_CXX) $(YC_CXXFLAGS) -x c++ -fPIC -c -o $@ $<
 	@ls -l $@
 
+$(YC_OBJ_DIR)/%.o: $(COEFF_DIR)/%.cpp $(YC_INC_GLOB)
+	$(MKDIR) $(YC_OBJ_DIR)
+	$(CXX_PREFIX) $(YC_CXX) $(YC_CXXFLAGS) -x c++ -fPIC -c -o $@ $<
+	@ls -l $@
+
 ######## Primary targets.
 
 default: compiler
@@ -120,7 +125,7 @@ $(YC_SWIG_OUT_DIR)/yask_compiler_api_wrap.cpp: $(YC_SWIG_DIR)/yask*.i $(INC_DIR)
 	$(SWIG) -version
 	$(MKDIR) $(YC_SWIG_OUT_DIR) $(PY_OUT_DIR)
 	$(SWIG) -v -DYC_MODULE=$(YC_MODULE) -cppext cpp \
-	  -I$(INC_DIR) -I$(COMM_DIR) -I$(COMM_DIR)/swig \
+	  -I$(INC_DIR) -I$(COMM_DIR) -I$(COMM_DIR)/swig -I$(COEFF_DIR) \
 	  -c++ -python -o $@ -outdir $(PY_OUT_DIR) -builtin $<
 
 # Turn off asserts to work around known SWIG issue:

src/contrib/coefficients/fd_coeff.cpp

@@ -28,80 +28,60 @@ IN THE SOFTWARE.
 
 #include <iostream>
 #include <cstring>
+#include <vector>
 #include "fd_coeff.hpp"
-#define MIN(x, y) (((x) < (y)) ? (x): (y))
-#define MAX(x, y) (((x) > (y)) ? (x): (y))
 
 using namespace std;
-/*input:coeff=empty coefficient array (or one that can be overwritten)
-	eval_point=point at which the derivative is approximated
-	order=order of the derivative to approximate (e.g. f'' corresponds to order = 2)
-	points=array of points from which to construct the approximation of the derivative. usually an equi-spaced array with points [-radius*h, -(radius-1)*h,...0, ... radius*h]
-	num_points=number of elements in points[], e.g. the number of points used to approximate the derivative.
-	Note: if num_points < order+1, then the coefficients will all be 0
-
-
-  output:void, fills the coefficient array such that
-f^(m)[eval_point] ~~ sum of coeff[i]*f[point[i]] from i = 0 to num_points-1
-*/
-
-void fd_coeff(float *coeff, const float eval_point, const int order, float *points, const int num_points)
-{
-    float c1, c2, c3;
-    float x_0=eval_point;
-    float center=0;
-
-
-
-
-//  float* d = (float*) malloc((order+1)*num_points*num_points*sizeof(float));
-    float d[(order+1)*num_points*num_points];
-    int m_idx = (order+1)*num_points;
-    int n_idx = num_points;
-
-    //array initializer 1
-    /*
-    memset(d, 0.f, sizeof(d));
-    */
-
-    //array initializer 2
-    int sizeofd = (order+1)*(num_points)*(num_points)*sizeof(float);
-    memset(d, 0.f, sizeofd);
-
-
-    //array initializer 3
-    /*
-    for(int m=0; m <= order; ++m){
-	for(int n=0; n< num_points; ++n){
-	    for(int v=0; v<num_points;++v){
-	    d[m*m_idx+n*n_idx+v]=0.f;
-	    }}}
+namespace yask {
+
+    // C-style interface with arbitrary point locations.
+    /* input:coeff=empty coefficient array (or one that can be overwritten)
+       eval_point=point at which the derivative is approximated order=order
+       of the derivative to approximate (e.g. f'' corresponds to order = 2)
+       points=array of points from which to construct the approximation of
+       the derivative.  usually an equi-spaced array with points [-radius*h,
+       -(radius-1)*h,...0, ... radius*h] num_points=number of elements in
+       points[], e.g. the number of points used to approximate the
+       derivative.  Note: if num_points < order+1, then the coefficients will
+       all be 0
+       
+       output:void, fills the coefficient array such that
+       f^(m)[eval_point] ~~ sum of coeff[i]*f[point[i]] from i = 0 to num_points-1
     */
-
-
-    d[0]=1.f;
-    c1 = 1.f;
-
-    for(int n=1; n<=num_points-1;++n){
-        c2=1.f;
-	for(int v=0; v<=n-1; ++v){
-            c3 = points[n] - points[v];
-            c2 = c2*c3;
-            for(int m=0; m<=MIN(n, order); ++m){
-		d[m*m_idx+n*n_idx + v] = (points[n]-x_0)*d[m*m_idx + (n-1)*n_idx + v] - m*d[(m-1)*m_idx + (n-1)*n_idx + v];
-		d[m*m_idx + n*n_idx + v] *= 1.f/c3;
+    void fd_coeff(double *coeff, const double eval_point, const int order, const double *points, const int num_points)
+    {
+        double c1, c2, c3;
+        double x_0=eval_point;
+
+        int m_idx = (order+1)*num_points;
+        int n_idx = num_points;
+        int dsz = m_idx * n_idx;
+        double d[dsz];
+
+        for (int i = 1; i < dsz; i++)
+            d[i] = 0.0;
+        d[0]=1.0;
+        c1 = 1.0;
+
+        for(int n=1; n<=num_points-1;++n){
+            c2=1.0;
+            for(int v=0; v<=n-1; ++v){
+                c3 = points[n] - points[v];
+                c2 = c2*c3;
+                for(int m=0; m<=min(n, order); ++m){
+                    d[m*m_idx+n*n_idx + v] = (points[n]-x_0)*d[m*m_idx + (n-1)*n_idx + v] - m*d[(m-1)*m_idx + (n-1)*n_idx + v];
+                    d[m*m_idx + n*n_idx + v] *= 1.0/c3;
+                }
             }
-	}
-	for(int m=0; m<= MIN(n, order); ++m){
-            d[m*m_idx+n*n_idx+n] = m*d[(m-1)*m_idx+(n-1)*n_idx+(n-1)] - (points[n-1]-x_0)*d[m*m_idx+(n-1)*n_idx+n-1];
-            d[m*m_idx+n*n_idx+n] *= c1/c2;
-	}
-        c1=c2;
-    }
+            for(int m=0; m<= min(n, order); ++m){
+                d[m*m_idx+n*n_idx+n] = m*d[(m-1)*m_idx+(n-1)*n_idx+(n-1)] - (points[n-1]-x_0)*d[m*m_idx+(n-1)*n_idx+n-1];
+                d[m*m_idx+n*n_idx+n] *= c1/c2;
+            }
+            c1=c2;
+        }
 
-    for(int i=0; i<num_points; ++i){
-	coeff[i] = d[order*m_idx+(num_points-1)*n_idx + i];
+        for(int i=0; i<num_points; ++i){
+            coeff[i] = d[order*m_idx+(num_points-1)*n_idx + i];
+        }
     }
-
-//    free(d);
 }

src/contrib/coefficients/fd_coeff.hpp

@@ -26,7 +26,9 @@ IN THE SOFTWARE.
 // Finite-differences coefficients code.
 // Contributed by Jeremy Tillay.
 
-#ifndef HEADER_COEFF
-#define HEADER_COEFF
-void fd_coeff(float *coeff, const float eval_point, const int order, float *points, const int num_points);
-#endif
+#pragma once
+
+namespace yask {
+    void fd_coeff(double *coeff, const double eval_point, const int order, const double *points, const int num_points);
+}
+