 }
 
 template <typename T>
+__global__ void normalize_sub_kernel(size_t n, size_t sub_n, T* x, size_t incx) {
+    auto index        = threadIdx.x + blockIdx.x * blockDim.x;
+    const auto stride = blockDim.x * gridDim.x;
+
+    for (; index < n; index += stride) {
+        // 1. compute the mean
+
+        T m = 0;
+
+        for(size_t i = 0; i < sub_n; ++i){
+            m += x[(index * sub_n + i) * incx];
+        }
+
+        m /= T(sub_n);
+
+        // 2. Compute the standard deviation
+
+        T s = 0;
+
+        for (size_t i = 0; i < sub_n; ++i) {
+            s += (x[(index * sub_n + i) * incx] - m) * (x[(index * sub_n + i) * incx] - m);
+        }
+
+        s = sqrtf(s / T(sub_n));
+
+        // 2. Normalize
+
+        for (size_t i = 0; i < sub_n; ++i) {
+            x[(index * sub_n + i) * incx] = (x[(index * sub_n + i) * incx] - m) / s;
+        }
+    }
+}
+
+template <typename T>
+__global__ void normalize_sub_kernel1(size_t n, size_t sub_n, T* x) {
+    auto index        = threadIdx.x + blockIdx.x * blockDim.x;
+    const auto stride = blockDim.x * gridDim.x;
+
+    for (; index < n; index += stride) {
+        // 1. compute the mean
+
+        T m = 0;
+
+        for(size_t i = 0; i < sub_n; ++i){
+            m += x[index * sub_n + i];
+        }
+
+        m /= T(sub_n);
+
+        // 2. Compute the standard deviation
+
+        T s = 0;
+
+        for (size_t i = 0; i < sub_n; ++i) {
+            s += (x[index * sub_n + i] - m) * (x[index * sub_n + i] - m);
+        }
+
+        s = sqrtf(s / T(sub_n));
+
+        // 2. Normalize
+
+        for (size_t i = 0; i < sub_n; ++i) {
+            x[index * sub_n + i] = (x[index * sub_n + i] - m) / s;
+        }
+    }
+}
+
+template <typename T>
 void normalize_flat_kernel1_run(size_t n, T* x, size_t incx, T mean) {
     static int blockSize   = 0;
     static int minGridSize = 0;

}

void egblas_snormalize_sub(size_t n, float* x, size_t sub_n, size_t incx) {
    for(size_t i = 0; i < n; ++i){
        egblas_snormalize_flat(sub_n, x + i * sub_n * incx, incx);
    }
}

void egblas_dnormalize_sub(size_t n, double* x, size_t sub_n, size_t incx) {
    for(size_t i = 0; i < n; ++i){
        egblas_dnormalize_flat(sub_n, x + i * sub_n * incx, incx);
    }
}

    cudaDeviceSynchronize();
}

template <typename T>
void normalize_sub_kernel_run(size_t n, size_t sub_n, T* x, size_t incx) {
    static int blockSize   = 0;
    static int minGridSize = 0;

    if (!blockSize) {
        cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, normalize_sub_kernel<T>, 0, 0);
    }

    int gridSize = ((n / incx) + blockSize - 1) / blockSize;
    normalize_sub_kernel<T><<<gridSize, blockSize>>>(n, sub_n, x, incx);

    cudaDeviceSynchronize();
}

template <typename T>
void normalize_sub_kernel1_run(size_t n, size_t sub_n, T* x) {
    static int blockSize   = 0;
    static int minGridSize = 0;

    if (!blockSize) {
        cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, normalize_sub_kernel1<T>, 0, 0);
    }

    int gridSize = (n + blockSize - 1) / blockSize;
    normalize_sub_kernel1<T><<<gridSize, blockSize>>>(n, sub_n, x);

    cudaDeviceSynchronize();
}

void egblas_snormalize_flat(size_t n, float* x, size_t incx) {
    auto m = egblas_smean(x, n, incx);
    normalize_flat_kernel1_run(n, x, incx, m);
