Write a custom CUDA kernel to implement the PairNorm operator supporting all four modes: 'None', 'PN', 'PN-SI', and 'PN-SCS'.

Modes Definitions:
1. 'None': Identity.
2. 'PN': x = s * (x - col_mean) / global_row_norm_avg
3. 'PN-SI': x = s * (x - col_mean) / row_norm(x - col_mean)
4. 'PN-SCS': x = s * x / row_norm(x) - col_mean

Optimization Strategy:
The implementation uses a two-stage approach to handle the dependency on global column means and the different normalization logic.

1. Stage 1: Column Statistics Kernel (Common for PN, PN-SI, PN-SCS)
   - A generic reduction kernel computes the Sum and Sum-of-Squares for each column.
   - Using `atomicAdd` allows this to scale to large Batch sizes (N).
   - Results are used to compute `col_mean` on the CPU.

2. Stage 2: Application Kernels (Mode-Specific)
   - **For 'PN'**: The normalization factor is a scalar derived from global stats. A fully vectorized (float4) element-wise kernel applies the centering and scaling.
   - **For 'PN-SI' and 'PN-SCS'**: These require per-row norms. We implement a **Row-wise Reduction Kernel**:
     - One Thread Block is assigned to one Row (Sample).
     - Threads efficiently load the row and the pre-computed column means into registers/shared memory.
     - A parallel block reduction computes the row's L2 norm (either of `x` or `x - col_mean` depending on mode).
     - Threads then apply the formula and write back.

3. Precision:
   - Accumulators use `double` precision to ensure stability when calculating variance and norms.
   - Inputs/Outputs are `float`.
  
Here's an example to show you the syntax of inline embedding custom CUDA operators in torch: The example given architecture is:   
  
```python
import torch
import torch.nn as nn

BATCH_SIZE = 4096
DIM = 4096
SHAPE = (BATCH_SIZE, DIM)
SCALE = 1
MODE = 'PN'

class PairNorm(nn.Module):
    '''
    PairNorm
    https://openreview.net/pdf?id=rkecl1rtwB
    '''
    def __init__(self, mode='PN', scale=1):
        """
            mode:
              'None' : No normalization 
              'PN'   : Original version
              'PN-SI'  : Scale-Individually version
              'PN-SCS' : Scale-and-Center-Simultaneously version
        """
        assert mode in ['None', 'PN',  'PN-SI', 'PN-SCS']
        super(PairNorm, self).__init__()
        self.mode = mode
        self.scale = scale
                
    def forward(self, x):
        if self.mode == 'None':
            return x
        
        col_mean = x.mean(dim=0)      
        if self.mode == 'PN':
            x = x - col_mean
            rownorm_mean = (1e-6 + x.pow(2).sum(dim=1).mean()).sqrt() 
            x = self.scale * x / rownorm_mean

        if self.mode == 'PN-SI':
            x = x - col_mean
            rownorm_individual = (1e-6 + x.pow(2).sum(dim=1, keepdim=True)).sqrt()
            x = self.scale * x / rownorm_individual

        if self.mode == 'PN-SCS':
            rownorm_individual = (1e-6 + x.pow(2).sum(dim=1, keepdim=True)).sqrt()
            x = self.scale * x / rownorm_individual - col_mean

        return x

class Model(nn.Module):
    def __init__(self, mode, scale):
        super(Model, self).__init__()
        self.norm = PairNorm(mode=mode, scale=scale)
    
    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return self.norm(x)

def get_inputs():
    x = torch.randn(SHAPE, dtype=torch.float32)
    return [x.contiguous()]

def get_init_inputs():
    return [MODE, SCALE]