mma_kernel.h

/*
 Copyright (c) 2021-2025, 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.
*/
 
#ifndef MMA_CUDA_KERNEL_H
#define MMA_CUDA_KERNEL_H
 
template <typename T>
__global__ void mma_Ljjxinv_kernel(T* __restrict__ Ljjxinv,
     const T* __restrict__ pjlambda, const T* __restrict__ qjlambda,
     const T* __restrict__ x, const T* __restrict__ low, const T* __restrict__ upp,
     const T* __restrict__ alpha, const T* __restrict__ beta,
     const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    const T xt = x[tj];
    T val = -1.0 / (2.0 * pjlambda[tj] / pow(upp[tj] - xt, 3) +
                    2.0 * qjlambda[tj] / pow(xt - low[tj], 3));
    // Remove the sensitivity for the active primal constraints
    bool is_alpha = xt == alpha[tj];
    bool is_beta  = xt == beta[tj];
    Ljjxinv[tj] = (is_alpha || is_beta) ? T(0.0) : val;
  }
}
 
 
template <typename T>
__global__ void mma_dipsolvesub1_kernel(T* __restrict__ x,
     const T* __restrict__ pjlambda, const T* __restrict__ qjlambda,
     const T* __restrict__ low, const T* __restrict__ upp,
     const T* __restrict__ alpha, const T* __restrict__ beta,
     const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    T pj = sqrt(pjlambda[tj]);
    T qj = sqrt(qjlambda[tj]);
    T denom = pj + qj;
    T val = (pj * low[tj] + qj * upp[tj]) / denom;
 
    // Clamp x between alpha and beta using branchless min/max
    x[tj] = fmax(fmin(val, beta[tj]), alpha[tj]);
  }
 
}
 
template <typename T>
__global__ void mattrans_v_mul_kernel(T* __restrict__ output,
     const T* __restrict__ pij, const T* __restrict__ lambda,
     const int m, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    output[tj] = 0.0;
    for (int i = 0; i < m; i++) {
      // output[tj] = output[tj] + pij[tj + i * n] * lambda[i];
      output[tj] = output[tj] + pij[i + tj * m] * lambda[i];
    }
  }
}
 
template <typename T>
__global__ void mma_sub1_kernel(T* __restrict__ xlow, T* __restrict__ xupp,
     const T* __restrict__ x, const T* __restrict__ xmin,
     const T* __restrict__ xmax, const T asyinit, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     T xgap = xmax[tj] - xmin[tj];
     xlow[tj] = x[tj] - asyinit * xgap;
     xupp[tj] = x[tj] + asyinit * xgap;
  }
}
 
 
template< typename T >
__global__ void mma_sub2_kernel(T* __restrict__ low, T* __restrict__ upp,
     const T* __restrict__ x, const T* __restrict__ xold1,
     const T* __restrict__ xold2, const T* __restrict__ xdiff,
     const T asydecr, const T asyincr, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj >= n) return;
 
  // Load data into registers for faster accessing compare to global memory 
  // when accessing repeatedly)
  const T xval     = x[tj];
  const T xold1val = xold1[tj];
  const T xold2val = xold2[tj];
  const T lowval   = low[tj];
  const T uppval   = upp[tj];
  const T xdiffval = xdiff[tj];
 
  // Compute the product
  const T prod = (xval - xold1val) * (xold1val - xold2val);
 
  // Compute asy_factor without branching
  T asy_factor = (prod < T(0)) ? asydecr :
                 (prod > T(0)) ? asyincr : T(1);
 
  // Update low and upp using fma (fused multiply-add) for numerical stability
  T new_low = fma(-asy_factor, (xold1val - lowval), xval);
  T new_upp = fma(asy_factor,  (uppval - xold1val), xval);
 
  // Apply bounds
  new_low = max(new_low, xval - T(10.0) * xdiffval);
  new_low = min(new_low, xval - T(0.01) * xdiffval);
 
  new_upp = min(new_upp, xval + T(10.0) * xdiffval);
  new_upp = max(new_upp, xval + T(0.01) * xdiffval);
 
  // Write results back
  low[tj] = new_low;
  upp[tj] = new_upp;
}
 
template< typename T >
__global__ void mma_sub3_kernel(const T* __restrict__ x,
     const T* __restrict__ df0dx, const T* __restrict__ dfdx,
     T* __restrict__ low, T* __restrict__ upp, const T* __restrict__ xmin,
   const T* __restrict__ xmax, T* __restrict__ alpha, T* __restrict__ beta,
   T* __restrict__ p0j, T* __restrict__ q0j, T* __restrict__ pij,
   T* __restrict__ qij, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     T xgap = xmax[tj] - xmin[tj];
     alpha[tj] = max(max(xmin[tj], low[tj] +
        0.1 * (x[tj] - low[tj])), x[tj] - 0.5 * xgap);
     beta[tj] = min(min(xmax[tj], upp[tj] - 0.1 * (upp[tj] - x[tj])), x[tj] +
        0.5 * xgap);
 
     p0j[tj] = pow(upp[tj] - x[tj], 2) * (1.001 * max(df0dx[tj], 0.0) +
        0.001 * max(-df0dx[tj], 0.0) + 0.00001 / max(0.00001, xgap));
     q0j[tj] = pow(x[tj] - low[tj], 2) * (0.001 * max(df0dx[tj], 0.0) +
        1.001 * max(-df0dx[tj], 0.0) + 0.00001 / max(0.00001, xgap));
        
     for (int i = 0; i < m; i++) {
        pij[i + tj*m] = pow(upp[tj] - x[tj], 2) *
         (1.001 * max(dfdx[i + tj*m], 0.0) + 0.001 *
         max(-dfdx[i + tj*m], 0.0) + 0.00001 / max(0.00001, xgap));
        qij[i + tj*m] = pow(x[tj] - low[tj], 2) *
         (0.001 * max(dfdx[i + tj*m], 0.0) + 1.001 *
       max(-dfdx[i + tj*m], 0.0) + 0.00001 / max(0.00001, xgap));
     }
  }
}
 
 
template< typename T >
__global__ void mma_sub4_kernel(const T* __restrict__ x, T* __restrict__ low,
     T* __restrict__ upp, T* __restrict__ pij, T* __restrict__ qij,
     T* __restrict__ temp, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     for (int i = 0; i < m; i++) {
        temp[i + tj*m] = pij[i + tj*m] / (upp[tj] - x[tj]) +
         qij[i + tj*m] / (x[tj] - low[tj]);
     }
  }
}
 
 
template <typename T>
__global__ void mma_max2_kernel(T* __restrict__ xsi, const T* __restrict__ x,
     T* __restrict__ alpha, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     xsi[tj] = max(1.0, 1.0 / (x[tj] - alpha[tj]));
  }
}
 
 
 
template <typename T>
__global__ void relambda_kernel(T* __restrict__ temp, const T* __restrict__ x,
     const T* __restrict__ xupp, const T* __restrict__ xlow,
     const T* __restrict__ pij, const T* __restrict__ qij, const int n,
   const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     for (int i = 0; i < m; i++) {
        temp[i + tj*m] = pij[i + tj*m] / (xupp[tj] - x[tj]) +
         qij[i + tj*m] / (x[tj] - xlow[tj]);
     }
  }
}
 
 
template <typename T>
__global__ void sub2cons2_kernel(T* __restrict__ a, const T* __restrict__ b,
     const T* __restrict__ c, const T* __restrict__ d,
     const T e, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
     a[tj] = b[tj]*(c[tj]-d[tj])-e;
  }
}
 
template< typename T>
__inline__ __device__ T max_reduce_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 >
__global__ void maxval_kernel(const T*  __restrict__ a, T *temp, 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 maxval = 0.0;
  for (int i = idx; i < n; i += str) {
     maxval = max(maxval, abs(a[i]));
  }
 
  maxval = max_reduce_warp<T>(maxval);
  if (lane == 0)
     shared[wid] = maxval;
  __syncthreads();
 
  maxval = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0.0;
  if (wid == 0)
     maxval = max_reduce_warp<T>(maxval);
 
  if (threadIdx.x == 0)
    temp[blockIdx.x] = maxval;
 
}
 
template <typename T>
__global__ void max_reduce_kernel(T*  __restrict__ bufred, const int n) {
 
  T maxval = 0.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)
  {
    maxval =max(maxval, bufred[i]);
  }
 
  __shared__ T shared[32];
  unsigned int lane = threadIdx.x % warpSize;
  unsigned int wid = threadIdx.x / warpSize;
 
  maxval = max_reduce_warp<T>(maxval);
  if (lane == 0)
    shared[wid] = maxval;
  __syncthreads();
 
  maxval = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0.0;
  if (wid == 0)
    maxval = max_reduce_warp<T>(maxval);
 
  if (threadIdx.x == 0)
    bufred[blockIdx.x] = maxval;
}
 
template <typename T>
__global__ void delx_kernel(T* __restrict__ delx, const T* __restrict__ x,
     const T* __restrict__ xlow, const T* __restrict__ xupp,
     const T* __restrict__ pij, const T* __restrict__ qij,
   const T* __restrict__ p0j, const T* __restrict__ q0j,
     const T* __restrict__ alpha, const T* __restrict__ beta,
   const T* __restrict__ lambda, const T epsi, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    delx[tj]=0;
    for (int i = 0; i < m; i++) {
      delx[tj] = delx[tj] + pij[i + tj * m] *
         lambda[i] / pow(xupp[tj] - x[tj], 2) -
           qij[i + tj * m] * lambda[i] / pow(x[tj] - xlow[tj], 2);
    }
    delx[tj] = delx[tj] + p0j[tj] / pow(xupp[tj] - x[tj], 2) -
       q0j[tj] / pow(x[tj] - xlow[tj], 2) - epsi / (x[tj] - alpha[tj])
         + epsi / (beta[tj] - x[tj]);
  }
}
 
 
template <typename T>
__global__ void GG_kernel(T* __restrict__ GG, const T* __restrict__ x,
     const T* __restrict__ xlow, const T* __restrict__ xupp,
     const T* __restrict__ pij, const T* __restrict__ qij, const int n,
     const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    for (int ggdumiter = 0; ggdumiter < m; ggdumiter++) {
      GG[ggdumiter  + m*tj] = pij[ggdumiter + m*tj] / pow(xupp[tj] - x[tj], 2) -
         qij[ggdumiter + m*tj] / pow(x[tj] - xlow[tj], 2);
    }
  }
}
template <typename T>
__global__ void diagx_kernel(T* __restrict__ diagx, const T* __restrict__ x,
     const T* __restrict__ xsi, const T* __restrict__ xlow,
     const T* __restrict__ xupp, const T* __restrict__ p0j,
   const T* __restrict__ q0j, const T* __restrict__ pij,
   const T* __restrict__ qij, const T* alpha, const T*  beta,
   const T*  eta, const T* lambda, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    T sum = 0;
    T sum1 = 0;
    for (int i = 0; i < m; i++) {
      sum = sum + pij[tj *m+ i] * lambda[i];
      sum1 = sum1 + qij[tj*m + i] * lambda[i];
    }
    diagx[tj] = (p0j[tj] + sum) / pow(xupp[tj] - x[tj], 3) +
       (q0j[tj] + sum1) / pow(x[tj] - xlow[tj], 3);
    diagx[tj] = 2.0 * diagx[tj] + xsi[tj] / (x[tj] - alpha[tj]) +
       eta[tj] / (beta[tj] - x[tj]);
  }
}
 
 
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>
__global__ void mmareduce_kernel(T* __restrict__ 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 mmasum_kernel(const T*  __restrict__ a, T*  __restrict__ buf_h,
     const int n, const int m, const int k) {
 
  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;
  for (int i = idx; i < n; i += str)
  {
    sum += a[m * i + k ];
  }
 
  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 >
__global__ void mmasumbb_kernel(const T*  __restrict__ GG,
     const T*  __restrict__ delx, const T*  __restrict__ diagx,
     T*  __restrict__ buf_h, const int n, const int m, const int k) {
 
  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;
  for (int i = idx; i < n; i += str)
  {
    sum += GG[ k + i * m] * delx[i] / diagx[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 >
__global__ void mmasumHess_kernel(const T*  __restrict__ hijx,
     const T*  __restrict__ Ljjxinv, T*  __restrict__ buf_h, const int n,
   const int m, const int k0, const int k1) {
 
  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;
  // this is similar to mmasumAA_kernel but with Ljjxinv_d = 1/diagx
  __shared__ T shared[32];
  T sum = 0;
  for (int i = idx; i < n; i += str)
  {
    sum += hijx[ k0 + i * m] * Ljjxinv[i]  * hijx[ k1 + i * m];
  }
 
  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 >
__global__ void mmasumAA_kernel(const T*  __restrict__ GG,
     const T*  __restrict__ diagx, T*  __restrict__ buf_h, const int n,
   const int m, const int k0, const int k1) {
 
  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;
  for (int i = idx; i < n; i += str)
  {
    sum += GG[ k0 + i * m] /diagx[i]  * GG[ k1 + i * m];
  }
 
  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>
__global__ void mma_copy_kernel(T* __restrict__ a, const T* __restrict__ b,
     const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
    if(tj<n)
      a[tj+m]=b[tj];
}
 
 
 
template <typename T>
__global__ void AA_kernel(T* __restrict__ temp, const T* __restrict__ GG,
     const T* __restrict__ diagx, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    for (int i0 = 0; i0 < m; i0++) {
      for (int i1 = 0; i1 < m; i1++) {
        temp[tj + i0 * (n + 1) + i1 * (m + 1) * (n + 1)] = GG[i0 * n + tj] *
         (1.0 / diagx[tj]) * GG[i1 * n + tj];
      }
    }
  }
}
 
 
template <typename T>
__global__ void dx_kernel(T* __restrict__ dx, const T* __restrict__ delx,
     const T* __restrict__ diagx, const T* __restrict__ GG,
     const T* __restrict__ dlambda, const int n, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    dx[tj] = -delx[tj]/diagx[tj];
    for(int i=0;i<m;i++){
      dx[tj] =dx[tj] - GG[tj*m+i]*dlambda[i]/diagx[tj];
    }
  }
}
 
 
 
template <typename T>
__global__ void dxsi_kernel(T* __restrict__ dxsi, const T* __restrict__ xsi,
     const T* __restrict__ dx, const T* __restrict__ x,
     const T* __restrict__ alpha, const T epsi, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    dxsi[tj]= -xsi[tj] + (epsi-dx[tj]*xsi[tj])/(x[tj] - alpha[tj]);
  }
}
 
 
template <typename T>
__global__ void deta_kernel(T* __restrict__ deta, const T* __restrict__ eta,
     const T* __restrict__ dx, const T* __restrict__ x,
     const T* __restrict__ beta, const T epsi, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    deta[tj] = -eta[tj] + (epsi + dx[tj] * eta[tj]) / (beta[tj] - x[tj]);
  }
}
 
 
 
template <typename T>
__global__ void RexCalculation_kernel(T* __restrict__ rex,
     const T* __restrict__ x, const T* __restrict__ xlow,
   const T* __restrict__ xupp, const T* __restrict__ pij,
   const T* __restrict__ p0j, const T* __restrict__ qij,
   const T* __restrict__ q0j, const T* __restrict__ lambda,
   const T* __restrict__ xsi, const T* __restrict__ eta, const int n,
   const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    rex[tj] = 0.0;
    for (int i = 0; i < m; i++) {
      rex[tj] = rex[tj] + pij[i +tj*m] * lambda[i] / pow(xupp[tj] - x[tj], 2) -
         qij[i +tj*m] * lambda[i] / pow(x[tj] - xlow[tj], 2);
    }
    rex[tj] = rex[tj] + p0j[tj] / pow(xupp[tj] - x[tj], 2) -
       q0j[tj] / pow(x[tj] - xlow[tj], 2) - xsi[tj] + eta[tj];
  }
}
 
template <typename T>
__global__ void rey_calculation_kernel(T* __restrict__ rey,
     const T* __restrict__ c, const T* __restrict__ d, const T* __restrict__ y,
     const T* __restrict__ lambda, const T* __restrict__ mu, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    rey[tj] = c[tj] + d[tj] * y[tj] - lambda[tj] - mu[tj];
  }
}
 
 
template< typename T >
__global__ void norm_kernel(const T*  __restrict__ a, T*  __restrict__ 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;
  for (int i = idx; i < n; i += str)
  {
    sum += pow(a[i], 2);
  }
 
  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>
__global__ void sub2cons_kernel(T* __restrict__ a, const T* __restrict__ b,
     const T* __restrict__ c,
  const T d, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    a[tj] = b[tj]*c[tj]-d;
  }
}
 
 
template <typename T>
__global__ void dely_kernel(T* __restrict__ dely, const T* __restrict__ c,
     const T* __restrict__ d, const T* __restrict__ y,
     const T* __restrict__ lambda, const T epsi, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    dely[tj] = c[tj] + d[tj]*y[tj] - lambda[tj] - epsi/y[tj];
  }
}
 
 
 
template< typename T >
__global__ void maxval2_kernel(const T* __restrict__ a, const T* __restrict__ b,
     T* __restrict__ temp, const T cons, 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 maxval = cons * a[0] / b[0];
  for (int i = idx; i < n; i += str)
  {
    maxval = max(maxval, cons * a[i] / b[i]);
  }
 
  maxval = max_reduce_warp<T>(maxval);
  if (lane == 0)
    shared[wid] = maxval;
  __syncthreads();
 
  maxval = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0.0;
  if (wid == 0)
    maxval = max_reduce_warp<T>(maxval);
 
  if (threadIdx.x == 0)
    temp[blockIdx.x] = maxval;
 
}
 
 
template< typename T >
__global__ void maxval3_kernel(const T* __restrict__ a, const T* __restrict__ b,
     const T* __restrict__ c, T* __restrict__ temp, const T cons, 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 maxval = cons * a[0] / b[0];
  for (int i = idx; i < n; i += str)
  {
    maxval = max(maxval, cons * a[i] / (b[i] - c[i]));
  }
 
  maxval = max_reduce_warp<T>(maxval);
  if (lane == 0)
    shared[wid] = maxval;
  __syncthreads();
 
  maxval = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
  if (wid == 0)
    maxval = max_reduce_warp<T>(maxval);
 
  if (threadIdx.x == 0)
    temp[blockIdx.x] = maxval;
 
}
 
 
template <typename T>
__global__ void kkt_rex_kernel(T* __restrict__ rex, const T* __restrict__ df0dx,
     const T* __restrict__ dfdx, const T* __restrict__ xsi,
     const T* __restrict__ eta, const T* __restrict__ lambda, const int n,
   const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    rex[tj] = 0.0;
    for (int i = 0; i < m; i++) {
      rex[tj] = rex[tj] + dfdx[i + tj*m] * lambda[i];
    }
    rex[tj] += df0dx[tj] - xsi[tj] + eta[tj];
  }
}
 
 
template <typename T>
__global__ void maxcons_kernel(T* __restrict__ a, const T b,
     const T c, const T* __restrict__ d, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    a[tj] = max(b, c * d[tj]);
  }
}
 
 template< typename T >
 __global__ void glsum_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 sum = 0;
   for (int i = idx; i<n ; i += str)
   {
     sum += a[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 >
__global__ void glsc2_kernel(const T * a,
                              const T * b,
                              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];
   }
 
   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>
__global__ void add2inv2_kernel(T* __restrict__ a, const T* __restrict__ b,
     const T c, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    a[tj] = a[tj]+c/b[tj];
  }
}
 
template <typename T>
__global__ void max2_kernel(T* __restrict__ a, const T b,
     const T* __restrict__ c, const T d, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if (tj < n) {
    a[tj]=max(b, d*c[tj]);
  }
}
 
template <typename T>
__global__ void updatebb_kernel(T* __restrict__ bb,
     const T* __restrict__ dellambda, const T* __restrict__ dely,
   const T* __restrict__ d, const T* __restrict__ mu,
   const T* __restrict__ y, const T delz, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if(tj<m)
    bb[tj]=dellambda[tj] + dely[tj]/(d[tj] + mu[tj]/y[tj]) - bb[tj];
  else if(tj<m+1)
    bb[tj]=delz;
}
 
 
 
template <typename T>
__global__ void updateAA_kernel(T* __restrict__ AA,
     const T* __restrict__ globaltmp_mm, const T* __restrict__ s,
   const T* __restrict__ lambda, const T* __restrict__ d,
     const T* __restrict__ mu, const T* __restrict__ y, const T* __restrict__ a,
   const T zeta, const T z, const int m) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if(tj<m)
    {
      AA[tj+tj*(m+1)]=globaltmp_mm[tj+tj*m] + (s[tj] / lambda[tj] +
         1.0/ (d[tj] + mu[tj] / y[tj]));
      AA[tj+m*(m+1)]=a[tj];
      AA[m+tj*(m+1)]=a[tj];
    }
  else if(tj<m+1)
    AA[tj+tj*(m+1)]= -zeta/z;
}
 
template <typename T>
__global__ void dy_kernel(T* __restrict__ dy, const T* __restrict__ dely,
     const T* __restrict__ dlambda, const T* __restrict__ d,
     const T* __restrict__ mu, const T* __restrict__ y, const int n) {
  int tj = blockIdx.x * blockDim.x + threadIdx.x;
  if(tj<n)
    dy[tj] = (-dely[tj]+dlambda[tj])/(d[tj] + mu[tj]/y[tj]);
}
 
#endif