math_kernel.h

Go to the documentation of this file.

#ifndef __MATH_MATH_KERNEL_H__
#define __MATH_MATH_KERNEL_H__
/*
 Copyright (c) 2021-2026, The Neko Authors
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions
 are met:
 
   * Redistributions of source code must retain the above copyright
     notice, this list of conditions and the following disclaimer.
 
   * Redistributions in binary form must reproduce the above
     copyright notice, this list of conditions and the following
     disclaimer in the documentation and/or other materials provided
     with the distribution.
 
   * Neither the name of the authors nor the names of its
     contributors may be used to endorse or promote products derived
     from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.
*/
 
#include <device/device_config.h>
 
template< typename T >
__global__ void cmult_kernel(T * __restrict__ a,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c * a[i];
  }
}
 
template< typename T >
__global__ void masked_gather_copy_kernel(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[i] = b[mask[i+1]-1];
  }
}
 
 
template< typename T >
__global__ void masked_gather_copy_aligned_kernel(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[i] = b[mask[i]];
  }
}
 
__device__ __forceinline__
void face_gather_nonlinear_index(int *index, const int idx, const int lx,
                                 const int ly, const int lz) {
  const int idx2 = idx - 1;
  index[3] = idx2 / (lx * ly * lz);
  index[2] = (idx2 - (lx * ly * lz) * index[3]) / (lx * ly);
  index[1] = (idx2 - (lx * ly * lz) * index[3] - (lx * ly) * index[2]) / lx;
  index[0] = (idx2 - (lx * ly * lz) * index[3] - (lx * ly) * index[2]) -
    lx * index[1];
  index[0]++;
  index[1]++;
  index[2]++;
  index[3]++;
}
 
__device__ __forceinline__
int face_gather_idx(const int i, const int j, const int k, const int l,
                    const int n1, const int n2, const int nf) {
  return ((i) + (n1) * (((j) - 1) + (n2) * (((k) - 1) + (nf) * (((l) - 1))))) - 1;
}
 
template< typename T >
__global__ void face_masked_gather_copy_kernel(T * __restrict__ a,
                                               const T * __restrict__ b,
                                               const int * __restrict__ mask,
                                               const int * __restrict__ facet,
                                               const int n1,
                                               const int n2,
                                               const int lx,
                                               const int ly,
                                               const int lz,
                                               const int n_mask) {
  int index[4];
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int m = idx; m < n_mask; m += str) {
    const int f = facet[m + 1];
    face_gather_nonlinear_index(index, mask[m + 1], lx, ly, lz);
 
    switch (f) {
    case 1:
    case 2:
      a[m] = b[face_gather_idx(index[1], index[2], f, index[3], n1, n2, 6)];
      break;
    case 3:
    case 4:
      a[m] = b[face_gather_idx(index[0], index[2], f, index[3], n1, n2, 6)];
      break;
    case 5:
    case 6:
      a[m] = b[face_gather_idx(index[0], index[1], f, index[3], n1, n2, 6)];
      break;
    }
  }
}
 
 
template< typename T >
__global__ void masked_scatter_copy_kernel(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[mask[i+1]-1] = b[i];
  }
}
 
template< typename T >
__global__ void masked_scatter_copy_aligned_kernel(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[mask[i]] = b[i];
  }
}
 
#if __CUDA_ARCH__ < 600
#include <cassert>
#endif
 
#define MASKED_ATOMIC_RED_BLOCK 256
 
template< typename T >
__global__ void masked_atomic_reduction_kernel(T * __restrict__ a,
                                               T * __restrict__ b,
                                               int * __restrict__ mask,
                                               const int n,
                                               const int m) {
 
#if __CUDA_ARCH__ >= 600
  __shared__ T bins[MASKED_ATOMIC_RED_BLOCK];
  __shared__ int bins_base;
 
  const int lane = threadIdx.x % 32;
  const int str = MASKED_ATOMIC_RED_BLOCK * gridDim.x;
 
  bins[threadIdx.x] = 0.0;
 
  /* Stride over whole blocks, so every thread reaches every barrier and
     every lane every shuffle; lanes past the end carry no target */
  for (int base = blockIdx.x * MASKED_ATOMIC_RED_BLOCK; base < m;
       base += str) {
 
    if (threadIdx.x == 0)
      bins_base = mask[base + 1] - 1 - MASKED_ATOMIC_RED_BLOCK / 2;
    __syncthreads();
    const int kbase = bins_base;
 
    const int i = base + threadIdx.x;
    int key = -1;
    T val = 0.0;
    if (i < m) {
      key = mask[i + 1] - 1;
      val = b[i];
    }
 
    /* A run starts where the target changes and ends in the lane before
       the next run's head */
    const int key_prev = __shfl_up_sync(0xffffffff, key, 1);
    const bool head = (lane == 0) || (key != key_prev);
    const unsigned int heads = __ballot_sync(0xffffffff, head);
 
    if (heads == 0xffffffff) {
      /* No two neighbouring lanes share a target, so there is nothing to
         combine: scatter directly, as a plain atomic kernel would. The
         ballot is the same in every lane, so the branch is warp-uniform. */
      if (key >= 0)
        atomicAdd(&a[key], val);
    }
    else {
      const unsigned int next = heads & ~((2u << lane) - 1u);
      const int run_end = next ? __ffs(next) - 2 : 31;
 
      /* Segmented suffix sum, the head ends up with the sum of its run */
#pragma unroll
      for (int off = 1; off < 32; off *= 2) {
        const T tmp = __shfl_down_sync(0xffffffff, val, off);
        if (lane + off <= run_end)
          val += tmp;
      }
 
      if (head && key >= 0) {
        const int s = key - kbase;
        if (s >= 0 && s < MASKED_ATOMIC_RED_BLOCK)
          atomicAdd(&bins[s], val);
        else
          atomicAdd(&a[key], val);
      }
    }
    __syncthreads();
 
    /* Flush the bins and clear them for the next stride; the barrier at
       the top of the loop orders this against the next accumulation */
    const T sum = bins[threadIdx.x];
    if (sum != (T) 0.0) {
      atomicAdd(&a[kbase + threadIdx.x], sum);
      bins[threadIdx.x] = 0.0;
    }
  }
#else
  if (blockIdx.x * blockDim.x + threadIdx.x == 0)
    assert(0 && "masked_atomic_reduction_kernel requires compute capability 6.0 or higher.");
#endif
}
 
template< typename T >
__global__ void masked_copy_kernel_0(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[mask[i+1]-1] = b[mask[i+1]-1];
  }
}
 
template< typename T >
__global__ void masked_copy_kernel_aligned(T * __restrict__ a,
                                   T * __restrict__ b,
                                   int * __restrict__ mask,
                                   const int n,
                                   const int n_mask) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n_mask; i += str) {
    a[mask[i]] = b[mask[i]];
  }
}
 
template <typename T>
__global__ void cfill_mask_kernel(T* __restrict__ a,
                                  const T c,
                                  const int size,
                                  int* __restrict__ mask,
                                  const int mask_size) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < mask_size; i += str) { a[mask[i]] = c; }
}
 
template< typename T >
__global__ void cmult2_kernel(T * __restrict__ a,
                 T * __restrict__ b,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c * b[i];
  }
}
 
template< typename T >
__global__ void cdiv_kernel(T * __restrict__ a,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c / a[i];
  }
}
 
template< typename T >
__global__ void cdiv2_kernel(T * __restrict__ a,
                 T * __restrict__ b,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c / b[i];
  }
}
 
template< typename T >
__global__ void cadd_kernel(T * __restrict__ a,
                            const T c,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + c;
  }
}
 
template< typename T >
__global__ void cadd2_kernel(T * __restrict__ a,
                             T * __restrict__ b,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] + c;
  }
}
 
template< typename T >
__global__ void cwrap_kernel(T * __restrict__ a,
                             const T min_val,
                             const T max_val,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
  const T l = max_val - min_val;
 
  for (int i = idx; i < n; i += str) {
    a[i] = min_val + fmod(fmod(a[i] - min_val, l) + l, l);
  }
}
 
template< typename T >
__global__ void sqrt_inplace_kernel(T * __restrict__ a,
                                    const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = sqrt(a[i]);
  }
}
 
template< typename T >
__global__ void power_kernel(T * __restrict__ ap,
                             const T * __restrict__ a,
                             const T p,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    ap[i] = pow(a[i], p);
  }
}
 
template< typename T >
__global__ void cfill_kernel(T * __restrict__ a,
                             const T c,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c;
  }
}
 
template< typename T >
__global__ void add2_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + b[i];
  }
}
 
template< typename T >
__global__ void add3_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const T * __restrict__ c,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] + c[i];
  }
}
 
template< typename T >
__global__ void add4_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const T * __restrict__ c,
                            const T * __restrict__ d,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] + c[i] + d[i];
  }
}
 
template< typename T >
__global__ void add2s1_kernel(T * __restrict__ a,
                              const T * __restrict__ b,
                              const T c1,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c1 * a[i] + b[i];
  }
}
 
template< typename T >
__global__ void add2s2_many_kernel(T  * __restrict__  x,
                                   const T ** p,
                                   const T * alpha,
                                   const int p_cur,
                                   const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
 
  for (int i = idx; i < n; i+= str) {
    T tmp = 0.0;
    for (int j = 0; j < p_cur; j ++) {
      tmp += p[j][i]*alpha[j];
    }
    x[i] += tmp;
  }
}
 
template< typename T >
__global__ void add2s2_kernel(T * __restrict__ a,
                              const T * __restrict__ b,
                              const T c1,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + c1 * b[i];
  }
}
 
template< typename T >
__global__ void addsqr2s2_kernel(T * __restrict__ a,
                                 const T * __restrict__ b,
                                 const T c1,
                                 const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + c1 * (b[i] * b[i]);
  }
}
 
template< typename T >
__global__ void add3s2_kernel(T * __restrict__ a,
                              const T * __restrict__ b,
                              const T * __restrict__ c,
                              const T c1,
                              const T c2,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c1 * b[i] + c2 * c[i];
  }
}
 
template< typename T >
__global__ void add4s3_kernel(T * __restrict__ a,
                              const T * __restrict__ b,
                              const T * __restrict__ c,
                              const T * __restrict__ d,
                              const T c1,
                              const T c2,
                              const T c3,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = c1 * b[i] + c2 * c[i] + c3 * d[i];
  }
}
 
template< typename T >
__global__ void add5s4_kernel(T * __restrict__ a,
                              const T * __restrict__ b,
                              const T * __restrict__ c,
                              const T * __restrict__ d,
                              const T * __restrict__ e,
                              const T c1,
                              const T c2,
                              const T c3,
                              const T c4,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + c1 * b[i] + c2 * c[i] + c3 * d[i] + c4 * e[i];
  }
}
 
template< typename T >
__global__ void invcol1_kernel(T * __restrict__ a,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
  const T one = 1.0;
 
  for (int i = idx; i < n; i += str) {
    a[i] = one / a[i];
  }
}
 
template< typename T >
__global__ void invcol2_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] / b[i];
  }
}
 
template< typename T >
__global__ void invcol3_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const T * __restrict__ c,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] / c[i];
  }
}
 
template< typename T >
__global__ void col2_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] * b[i];
  }
}
 
template< typename T >
__global__ void col3_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const T * __restrict__ c,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] * c[i];
  }
}
 
template< typename T >
__global__ void subcol3_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const T * __restrict__ c,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] - b[i] * c[i];
  }
}
 
template< typename T >
__global__ void sub2_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] - b[i];
  }
}
 
template< typename T >
__global__ void sub3_kernel(T * __restrict__ a,
                            const T * __restrict__ b,
                            const T * __restrict__ c,
                            const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = b[i] - c[i];
  }
}
 
template< typename T >
__global__ void addcol3_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const T * __restrict__ c,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + b[i] * c[i];
  }
 
}
 
template< typename T >
__global__ void addcol4_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const T * __restrict__ c,
                               const T * __restrict__ d,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + b[i] * c[i] * d[i];
  }
 
}
 
template< typename T >
__global__ void addcol3s2_kernel(T * __restrict__ a,
                               const T * __restrict__ b,
                               const T * __restrict__ c,
                               const T s,
                               const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = a[i] + s * b[i] * c[i];
  }
 
}
 
template< typename T >
__global__ void vdot3_kernel(T * __restrict__ dot,
                             const T * __restrict__ u1,
                             const T * __restrict__ u2,
                             const T * __restrict__ u3,
                             const T * __restrict__ v1,
                             const T * __restrict__ v2,
                             const T * __restrict__ v3,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    dot[i] = u1[i] * v1[i]  + u2[i] * v2[i] + u3[i] * v3[i];
  }
 
}
 
template< typename T >
__global__ void vcross_kernel(T * __restrict__ u1,
                             T * __restrict__ u2,
                             T * __restrict__ u3,
                             const T * __restrict__ v1,
                             const T * __restrict__ v2,
                             const T * __restrict__ v3,
                             const T * __restrict__ w1,
                             const T * __restrict__ w2,
                             const T * __restrict__ w3,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    u1[i] = v2[i]*w3[i] - v3[i]*w2[i];
    u2[i] = v3[i]*w1[i] - v1[i]*w3[i];
    u3[i] = v1[i]*w2[i] - v2[i]*w1[i];
  }
 
}
 
 
template< typename T>
__inline__ __device__ T reduce_warp(T val) {
  val += __shfl_down_sync(0xffffffff, val, 16);
  val += __shfl_down_sync(0xffffffff, val, 8);
  val += __shfl_down_sync(0xffffffff, val, 4);
  val += __shfl_down_sync(0xffffffff, val, 2);
  val += __shfl_down_sync(0xffffffff, val, 1);
  return val;
}
 
template< typename T>
__inline__ __device__ T reduce_max_warp(T val) {
  val = max(val, __shfl_down_sync(0xffffffff, val, 16));
  val = max(val, __shfl_down_sync(0xffffffff, val, 8));
  val = max(val, __shfl_down_sync(0xffffffff, val, 4));
  val = max(val, __shfl_down_sync(0xffffffff, val, 2));
  val = max(val, __shfl_down_sync(0xffffffff, val, 1));
  return val;
}
 
template< typename T>
__inline__ __device__ T reduce_min_warp(T val) {
  val = min(val, __shfl_down_sync(0xffffffff, val, 16));
  val = min(val, __shfl_down_sync(0xffffffff, val, 8));
  val = min(val, __shfl_down_sync(0xffffffff, val, 4));
  val = min(val, __shfl_down_sync(0xffffffff, val, 2));
  val = min(val, __shfl_down_sync(0xffffffff, val, 1));
  return val;
}
 
template< typename T >
__global__ void reduce_kernel(T * bufred, const int n) {
 
  T sum = 0;
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
  for (int i = idx; i<n ; i += str)
  {
    sum += bufred[i];
  }
 
  __shared__ T shared[32];
  unsigned int lane = threadIdx.x % warpSize;
  unsigned int wid = threadIdx.x / warpSize;
 
  sum = reduce_warp<T>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T>(sum);
 
  if (threadIdx.x == 0)
    bufred[blockIdx.x] = sum;
}
 
template< typename T >
__global__ void reduce_max_kernel(T * bufred, const T ninf, const int n) {
 
  T max_val = ninf;
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
  for (int i = idx; i<n ; i += str)
  {
    max_val = max(max_val, bufred[i]);
  }
 
  __shared__ T shared[32];
  unsigned int lane = threadIdx.x % warpSize;
  unsigned int wid = threadIdx.x / warpSize;
 
  max_val = reduce_max_warp<T>(max_val);
  if (lane == 0)
    shared[wid] = max_val;
  __syncthreads();
 
  max_val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : ninf;
  if (wid == 0)
    max_val = reduce_max_warp<T>(max_val);
 
  if (threadIdx.x == 0)
    bufred[blockIdx.x] = max_val;
}
 
template< typename T >
__global__ void reduce_min_kernel(T * bufred, const T pinf, const int n) {
 
  T min_val = pinf;
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
  for (int i = idx; i<n ; i += str)
  {
    min_val = min(min_val, bufred[i]);
  }
 
  __shared__ T shared[32];
  unsigned int lane = threadIdx.x % warpSize;
  unsigned int wid = threadIdx.x / warpSize;
 
  min_val = reduce_min_warp<T>(min_val);
  if (lane == 0)
    shared[wid] = min_val;
  __syncthreads();
 
  min_val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : pinf;
  if (wid == 0)
    min_val = reduce_min_warp<T>(min_val);
 
  if (threadIdx.x == 0)
    bufred[blockIdx.x] = min_val;
}
 
template< typename T_acc >
__global__ void glsc3_reduce_kernel( T_acc * bufred,
                                     const int n,
                                     const int j
                                   ) {
   __shared__ T_acc buf[1024] ;
   const int idx = threadIdx.x;
   const int y= blockIdx.x;
   const int step = blockDim.x;
 
   buf[idx]=0;
   for (int i=idx ; i<n ; i+=step)
   {
     buf[idx] += bufred[i*j + y];
   }
   __syncthreads();
 
   int i = 512;
   while (i != 0)
   {
     if(threadIdx.x < i && (threadIdx.x + i) < n )
     {
        buf[threadIdx.x] += buf[threadIdx.x + i] ;
     }
     i = i>>1;
     __syncthreads();
   }
 
   bufred[y] = buf[0];
}
 
 
template< typename T >
__global__ void vlsc3_kernel(const T * a,
                             const T * b,
                             const T * c,
                             T * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T shared[32];
  T sum = 0.0;
  for (int i = idx; i < n; i+= str) {
    sum += a[i] * b[i] * c[i];
  }
 
  sum = reduce_warp<T>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T>(sum);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = sum;
}
 
template< typename T, typename T_acc >
__global__ void glsc3_kernel(const T * a,
                             const T * b,
                             const T * c,
                             T_acc * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T_acc shared[32];
  T_acc sum = 0.0;
  for (int i = idx; i < n; i+= str) {
    sum += static_cast<T_acc>(a[i] * b[i] * c[i]);
  }
 
  sum = reduce_warp<T_acc>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T_acc>(sum);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = sum;
}
 
template< typename T, typename T_acc >
__global__ void glsc3_many_kernel(const T * a,
                                  const T ** b,
                                  const T * c,
                                  T_acc * buf_h,
                                  const int j,
                                  const int n) {
 
  const int idx = blockIdx.x * blockDim.y + threadIdx.y;
  const int str = blockDim.y * gridDim.x;
  const int y = threadIdx.x;
 
  __shared__ T_acc buf[1024];
  T_acc tmp = 0;
  if(y < j){
    for (int i = idx; i < n; i+= str) {
      tmp += static_cast<T_acc>(a[i] * b[threadIdx.x][i] * c[i]);
    }
  }
 
  buf[threadIdx.y*blockDim.x+y] = tmp;
  __syncthreads();
 
  int i = blockDim.y>>1;
  while (i != 0) {
    if (threadIdx.y < i) {
      buf[threadIdx.y*blockDim.x +y] += buf[(threadIdx.y + i)*blockDim.x+y];
    }
    __syncthreads();
    i = i>>1;
  }
  if (threadIdx.y == 0) {
    if( y < j) {
      buf_h[j*blockIdx.x+y] = buf[y];
    }
  }
}
 
template< typename T, typename T_acc >
__global__ void glsc2_kernel(const T * a,
                             const T * b,
                             T_acc * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T_acc shared[32];
  T_acc sum = 0.0;
  for (int i = idx; i < n; i+= str) {
    sum += static_cast<T_acc>(a[i] * b[i]);
  }
 
  sum = reduce_warp<T_acc>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T_acc>(sum);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = sum;
 
}
 
template< typename T, typename T_acc >
__global__ void glsubnorm2_kernel(const T * a,
                                  const T * b,
                                  T_acc * buf_h,
                                  const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T_acc shared[32];
  T_acc sum = 0.0;
  for (int i = idx; i < n; i+= str) {
    sum += static_cast<T_acc>(pow(a[i] - b[i], 2.0));
  }
 
  sum = reduce_warp<T_acc>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T_acc>(sum);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = sum;
 
}
 
template< typename T, typename T_acc >
__global__ void glsum_kernel(const T * a,
                             T_acc * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T_acc shared[32];
  T_acc sum = 0;
  for (int i = idx; i<n ; i += str)
  {
    sum += static_cast<T_acc>(a[i]);
  }
 
  sum = reduce_warp<T_acc>(sum);
  if (lane == 0)
    shared[wid] = sum;
  __syncthreads();
 
  sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    sum = reduce_warp<T_acc>(sum);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = sum;
 
}
 
template< typename T >
__global__ void glmax_kernel(const T * a,
                             const T ninf,
                             T * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T shared[32];
  T max_val = ninf;
  for (int i = idx; i<n ; i += str)
  {
    max_val = max(max_val, a[i]);
  }
 
  max_val = reduce_max_warp<T>(max_val);
  if (lane == 0)
    shared[wid] = max_val;
  __syncthreads();
 
  max_val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : ninf;
  if (wid == 0)
    max_val = reduce_max_warp<T>(max_val);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = max_val;
 
}
 
template< typename T >
__global__ void glamax_kernel(const T * a,
                              T * buf_h,
                              const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T shared[32];
  T max_val = T(0);
  for (int i = idx; i<n ; i += str)
  {
    max_val = max(max_val, fabs(a[i]));
  }
 
  max_val = reduce_max_warp<T>(max_val);
  if (lane == 0)
    shared[wid] = max_val;
  __syncthreads();
 
  max_val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : T(0);
  if (wid == 0)
    max_val = reduce_max_warp<T>(max_val);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = max_val;
 
}
 
template< typename T >
__global__ void glmin_kernel(const T * a,
                             const T pinf,
                             T * buf_h,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  const unsigned int lane = threadIdx.x % warpSize;
  const unsigned int wid = threadIdx.x / warpSize;
 
  __shared__ T shared[32];
  T min_val = pinf;
  for (int i = idx; i<n ; i += str)
  {
    min_val = min(min_val, a[i]);
  }
 
  min_val = reduce_min_warp<T>(min_val);
  if (lane == 0)
    shared[wid] = min_val;
  __syncthreads();
 
  min_val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : pinf;
  if (wid == 0)
    min_val = reduce_min_warp<T>(min_val);
 
  if (threadIdx.x == 0)
    buf_h[blockIdx.x] = min_val;
 
}
 
template< typename T >
__global__ void absval_kernel(T * __restrict__ a,
                             const int n) {
 
  const int idx = blockIdx.x * blockDim.x + threadIdx.x;
  const int str = blockDim.x * gridDim.x;
 
  for (int i = idx; i < n; i += str) {
    a[i] = fabs(a[i]);
  }
}
 
// ========================================================================== //
// Kernels for the point-wise operations
 
template <typename T>
__global__ void
    pwmax_vec2_kernel(T* __restrict__ a, const T* __restrict__ b, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = max(a[i], b[i]);
}
 
template <typename T>
__global__ void pwmax_vec3_kernel(
    T* __restrict__ a, const T* __restrict__ b, const T* __restrict__ c,
    const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = max(b[i], c[i]);
}
 
template <typename T>
__global__ void pwmax_sca2_kernel(T* __restrict__ a, const T c, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = max(a[i], c);
}
 
template <typename T>
__global__ void pwmax_sca3_kernel(
    T* __restrict__ a, const T* __restrict b, const T c, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = max(b[i], c);
}
 
template <typename T>
__global__ void
    pwmin_vec2_kernel(T* __restrict__ a, const T* __restrict__ b, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = min(a[i], b[i]);
}
 
template <typename T>
__global__ void pwmin_vec3_kernel(
    T* __restrict__ a, const T* __restrict__ b, const T* __restrict__ c,
    const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = min(b[i], c[i]);
}
 
template <typename T>
__global__ void pwmin_sca2_kernel(T* __restrict__ a, const T c, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = min(a[i], c);
}
 
template <typename T>
__global__ void pwmin_sca3_kernel(
    T* __restrict__ a, const T* __restrict b, const T c, const int n) {
 
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    const int str = blockDim.x * gridDim.x;
 
    for (int i = idx; i < n; i += str) a[i] = min(b[i], c);
}
 
#endif // __MATH_MATH_KERNEL_H__