# HG changeset patch
# User Shan Kang <shan.kang@intel.com>
# Date 1612267147 -28800
#      Tue Feb 02 19:59:07 2021 +0800
# Node ID 25ad95899eb049ba4e1b4df62b99d51a617a4cf2
# Parent  df3a18e894881c6cf9b37ddbe0689e8782273d02
Partial specialization on do_bsxfun_op for double and float data

* liboctave/array/dNDArray.cc (do_bsxfun_op): optimized do_bsxfun_op
for double data through VML.
* liboctave/array/fNDArray.cc (do_bsxfun_op): optimized do_bsxfun_op
for float data through VML.

diff -r df3a18e89488 -r 25ad95899eb0 liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc	Tue Feb 02 19:52:17 2021 +0800
+++ b/liboctave/array/dNDArray.cc	Tue Feb 02 19:59:07 2021 +0800
@@ -27,6 +27,13 @@
 #  include "config.h"
 #endif
 
+#if defined (USE_VML)
+#if defined (NO_IFUNC_SUPPORT)
+#define GCC_IFUN_UNAVAILABLE 1
+#endif
+#  include "vml_math.h"
+#endif
+
 #include <istream>
 #include <limits>
 #include <ostream>
@@ -44,6 +51,175 @@
 
 #include "bsxfun-defs.cc"
 
+#if defined (USE_VML)
+void *
+octave_quit_wrapper (void *input)
+{
+  octave_quit ();
+  return NULL;
+}
+
+template <>
+inline void
+mx_inline_pow<> (size_t n, double *r, const double *x, const double *y)
+{
+  vml_vpow_hook (x, y, r, n, octave_quit_wrapper, NULL, NULL, NULL);
+}
+
+template <>
+Array<double> do_bsxfun_op <> (const Array<double>& x, const Array<double>& y,
+              void (*op_vv) (size_t, double *, const double *, const double *),
+              void (*op_sv) (size_t, double *, double, const double *),
+              void (*op_vs) (size_t, double *, const double *, double))
+{
+  int nd = std::max (x.ndims (), y.ndims ());
+  dim_vector dvx = x.dims ().redim (nd);
+  dim_vector dvy = y.dims ().redim (nd);
+  void (*fp_vv) (size_t, double *, const double *, const double *);
+  void (*fp_sv) (size_t, double *, double, const double *);
+  void (*fp_vs) (size_t, double *, const double *, double);
+  fp_vv = mx_inline_pow;
+  fp_sv = mx_inline_pow;
+  fp_vs = mx_inline_pow;
+
+  // Construct the result dimensions.
+  dim_vector dvr;
+  dvr.resize (nd);
+  for (int i = 0; i < nd; i++)
+    {
+      octave_idx_type xk = dvx(i);
+      octave_idx_type yk = dvy(i);
+      if (xk == 1)
+        dvr(i) = yk;
+      else if (yk == 1 || xk == yk)
+        dvr(i) = xk;
+      else
+        (*current_liboctave_error_handler)
+          ("bsxfun: nonconformant dimensions: %s and %s",
+           x.dims ().str ().c_str (), y.dims ().str ().c_str ());
+    }
+
+  Array<double> retval (dvr);
+  Array<double> tmpval1 (dvr);
+  const double *xvec = x.fortran_vec ();
+  const double *yvec = y.fortran_vec ();
+  double *rvec = retval.fortran_vec ();
+  double *tmpvec1 = tmpval1.fortran_vec ();
+
+  // Fold the common leading dimensions.
+  octave_idx_type start, ldr = 1;
+  for (start = 0; start < nd; start++)
+    {
+      if (dvx(start) != dvy(start))
+        break;
+      ldr *= dvr(start);
+    }
+
+  if (retval.isempty ())
+    ; // do nothing
+  else if (start == nd){
+    op_vv (retval.numel (), rvec, xvec, yvec);
+  }
+  else
+    {
+      // Determine the type of the low-level loop.
+      bool xsing = false;
+      bool ysing = false;
+      if (ldr == 1)
+        {
+          xsing = dvx(start) == 1;
+          ysing = dvy(start) == 1;
+          if (xsing || ysing)
+            {
+              ldr *= dvx(start) * dvy(start);
+              start++;
+            }
+        }
+      dim_vector cdvx = dvx.cumulative ();
+      dim_vector cdvy = dvy.cumulative ();
+      // Nullify singleton dims to achieve a spread effect.
+      for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
+        {
+          if (dvx(i) == 1)
+            cdvx(i-1) = 0;
+          if (dvy(i) == 1)
+            cdvy(i-1) = 0;
+        }
+
+        octave_idx_type niter = dvr.numel (start);
+        // The index array.
+        OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
+        size_t v_idx1 = 0;
+        size_t v_idx2 = 0;
+
+        if(op_vv==fp_vv || op_sv==fp_sv || op_vs==fp_vs)
+        {
+          for (octave_idx_type iter = 0; iter < niter; iter++)
+            {
+              octave_quit ();
+
+              // Compute indices.
+              octave_idx_type xidx = cdvx.cum_compute_index (idx);
+              octave_idx_type yidx = cdvy.cum_compute_index (idx);
+              octave_idx_type ridx = dvr.compute_index (idx);
+
+              // Apply the low-level loop.
+              if (xsing){
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = xvec[xidx];
+                  rvec[v_idx2++] = *(yvec + yidx + i);
+                }
+              }
+              else if (ysing){
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = *(xvec + xidx + i);
+                  rvec[v_idx2++] = yvec[yidx];
+                }
+              }
+              else{
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = *(xvec + xidx + i);
+                  rvec[v_idx2++] = *(yvec + yidx + i);
+                }
+              }
+              dvr.increment_index (idx + start, start);
+          }
+          op_vv(retval.numel (), rvec, tmpvec1, rvec);
+        } else
+        {
+          for (octave_idx_type iter = 0; iter < niter; iter++)
+          {
+            octave_quit ();
+
+            // Compute indices.
+            // FIXME: performance impact noticeable?
+            octave_idx_type xidx = cdvx.cum_compute_index (idx);
+            octave_idx_type yidx = cdvy.cum_compute_index (idx);
+            octave_idx_type ridx = dvr.compute_index (idx);
+
+            // Apply the low-level loop.
+            if (xsing){
+              op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx);
+            }
+            else if (ysing){
+              op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx]);
+            }
+            else{
+              op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx);
+            }
+
+            dvr.increment_index (idx + start, start);
+          }
+        }
+    }
+
+  return retval;
+}
+#endif
+
 NDArray::NDArray (const Array<octave_idx_type>& a, bool zero_based,
                   bool negative_to_nan)
 {
diff -r df3a18e89488 -r 25ad95899eb0 liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc	Tue Feb 02 19:52:17 2021 +0800
+++ b/liboctave/array/fNDArray.cc	Tue Feb 02 19:59:07 2021 +0800
@@ -44,6 +44,176 @@
 
 #include "bsxfun-defs.cc"
 
+#if defined (USE_VML)
+extern "C" {
+   typedef void *(*vml_hook_func_t1)(void *);
+   int vml_vpowf_hook(const float *input_array, const float *input_array1,
+   float *result_array, unsigned int size, vml_hook_func_t1 pre_hook,
+   void *pre_args, vml_hook_func_t1 post_hook, void *post_args);
+}
+extern void *octave_quit_wrapper (void *input);
+
+template <>
+inline void
+mx_inline_pow<> (size_t n, float *r, const float *x, const float *y)
+{
+  vml_vpowf_hook (x, y, r, n, octave_quit_wrapper, NULL, NULL, NULL);
+}
+
+template <>
+Array<float> do_bsxfun_op <> (const Array<float>& x, const Array<float>& y,
+              void (*op_vv) (size_t, float *, const float *, const float *),
+              void (*op_sv) (size_t, float *, float, const float *),
+              void (*op_vs) (size_t, float *, const float *, float))
+{
+  int nd = std::max (x.ndims (), y.ndims ());
+  dim_vector dvx = x.dims ().redim (nd);
+  dim_vector dvy = y.dims ().redim (nd);
+  void (*fp_vv) (size_t, float *, const float *, const float *);
+  void (*fp_sv) (size_t, float *, float, const float *);
+  void (*fp_vs) (size_t, float *, const float *, float);
+  fp_vv = mx_inline_pow;
+  fp_sv = mx_inline_pow;
+  fp_vs = mx_inline_pow;
+
+  // Construct the result dimensions.
+  dim_vector dvr;
+  dvr.resize (nd);
+  for (int i = 0; i < nd; i++)
+    {
+      octave_idx_type xk = dvx(i);
+      octave_idx_type yk = dvy(i);
+      if (xk == 1)
+        dvr(i) = yk;
+      else if (yk == 1 || xk == yk)
+        dvr(i) = xk;
+      else
+        (*current_liboctave_error_handler)
+          ("bsxfun: nonconformant dimensions: %s and %s",
+           x.dims ().str ().c_str (), y.dims ().str ().c_str ());
+    }
+
+  Array<float> retval (dvr);
+  Array<float> tmpval1 (dvr);
+  const float *xvec = x.fortran_vec ();
+  const float *yvec = y.fortran_vec ();
+  float *rvec = retval.fortran_vec ();
+  float *tmpvec1 = tmpval1.fortran_vec ();
+
+  // Fold the common leading dimensions.
+  octave_idx_type start, ldr = 1;
+  for (start = 0; start < nd; start++)
+    {
+      if (dvx(start) != dvy(start))
+        break;
+      ldr *= dvr(start);
+    }
+
+  if (retval.isempty ())
+    ; // do nothing
+  else if (start == nd){
+    op_vv (retval.numel (), rvec, xvec, yvec);
+  }
+  else
+    {
+      // Determine the type of the low-level loop.
+      bool xsing = false;
+      bool ysing = false;
+      if (ldr == 1)
+        {
+          xsing = dvx(start) == 1;
+          ysing = dvy(start) == 1;
+          if (xsing || ysing)
+            {
+              ldr *= dvx(start) * dvy(start);
+              start++;
+            }
+        }
+      dim_vector cdvx = dvx.cumulative ();
+      dim_vector cdvy = dvy.cumulative ();
+      // Nullify singleton dims to achieve a spread effect.
+      for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
+        {
+          if (dvx(i) == 1)
+            cdvx(i-1) = 0;
+          if (dvy(i) == 1)
+            cdvy(i-1) = 0;
+        }
+
+        octave_idx_type niter = dvr.numel (start);
+        // The index array.
+        OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
+        size_t v_idx1 = 0;
+        size_t v_idx2 = 0;
+
+        if(op_vv==fp_vv || op_sv==fp_sv || op_vs==fp_vs)
+        {
+          for (octave_idx_type iter = 0; iter < niter; iter++)
+            {
+              octave_quit ();
+
+              // Compute indices.
+              octave_idx_type xidx = cdvx.cum_compute_index (idx);
+              octave_idx_type yidx = cdvy.cum_compute_index (idx);
+              octave_idx_type ridx = dvr.compute_index (idx);
+
+              // Apply the low-level loop.
+              if (xsing){
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = xvec[xidx];
+                  rvec[v_idx2++] = *(yvec + yidx + i);
+                }
+              }
+              else if (ysing){
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = *(xvec + xidx + i);
+                  rvec[v_idx2++] = yvec[yidx];
+                }
+              }
+              else{
+                for(int i = 0; i < ldr; i++)
+                {
+                  tmpvec1[v_idx1++] = *(xvec + xidx + i);
+                  rvec[v_idx2++] = *(yvec + yidx + i);
+                }
+              }
+              dvr.increment_index (idx + start, start);
+          }
+          op_vv(retval.numel (), rvec, tmpvec1, rvec);
+        } else
+        {
+          for (octave_idx_type iter = 0; iter < niter; iter++)
+          {
+            octave_quit ();
+
+            // Compute indices.
+            // FIXME: performance impact noticeable?
+            octave_idx_type xidx = cdvx.cum_compute_index (idx);
+            octave_idx_type yidx = cdvy.cum_compute_index (idx);
+            octave_idx_type ridx = dvr.compute_index (idx);
+
+            // Apply the low-level loop.
+            if (xsing){
+              op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx);
+            }
+            else if (ysing){
+              op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx]);
+            }
+            else{
+              op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx);
+            }
+
+            dvr.increment_index (idx + start, start);
+          }
+        }
+    }
+
+  return retval;
+}
+#endif
+
 FloatNDArray::FloatNDArray (const charNDArray& a)
   : MArray<float> (a.dims ())
 {
