You write custom CUDA kernels to replace the pytorch operators in the given architecture to get speedups.  

You have complete freedom to choose the set of operators you want to replace. You may make the decision to replace some operators with custom CUDA kernels and leave others unchanged. You may replace multiple operators with custom implementations, consider operator fusion opportunities (combining multiple operators into a single kernel, for example, combining matmul+relu), or algorithmic changes (such as online softmax). You are only limited by your imagination.  

**SPECIAL INSTRUCTIONS FOR REDUCTION-INTENSIVE OPERATORS (Loss Functions):**

When the target operator involves reduction operations that sum across tensor dimensions (like loss functions), you MUST implement the following optimized strategy:

1. **TWO-LEVEL REDUCTION APPROACH**: Instead of using atomic operations on the final result, implement:
   - **Level 1**: Block-level reduction using shared memory
   - **Level 2**: Final reduction on GPU using torch.sum() on partial results
   - This eliminates atomic operation contention for large tensors

2. **SHARED MEMORY REDUCTION PATTERN**:
   cpp
// Each thread accumulates locally
float local_loss = 0.0f;
for (int element_idx = idx; element_idx < total_elements; element_idx += stride) {
// Compute loss for each element
local_loss += computed_loss;
}

// Block-level reduction using shared memory
extern shared float shared_mem[];
shared_mem[threadIdx.x] = local_loss;
__syncthreads();

// Reduction within block
for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
if (threadIdx.x < stride) {
shared_mem[threadIdx.x] += shared_mem[threadIdx.x + stride];
}
__syncthreads();
}

// Each block writes partial sum
if (threadIdx.x == 0) {
partial_sums[blockIdx.x] = shared_mem[0];
}



3. **DYNAMIC BLOCK ALLOCATION**: Based on total_elements:
cpp
int num_blocks;
if (total_elements <= 2048) num_blocks = 4;
else if (total_elements <= 8192) num_blocks = 8;
else if (total_elements <= 32768) num_blocks = 16;
else if (total_elements <= 131072) num_blocks = 32;
else num_blocks = 64;



4. **ELEMENT PROCESSING STRATEGY**: Each thread processes multiple elements:
cpp
int stride = blockDim.x * gridDim.x;
for (int element_idx = idx; element_idx < total_elements; element_idx += stride) {
// Process element at element_idx
}



5. **FINAL REDUCTION ON GPU**: Use torch.sum() on partial results:
cpp
auto partial_sums = torch::zeros({num_blocks}, input.options());
// … kernel execution …
auto total_loss = torch::sum(partial_sums);
return total_loss;



6. **COMPILATION FLAGS**: Use balanced optimization:
python
extra_cuda_cflags=[
“-O3”,
“–use_fast_math”,
“-gencode=arch=compute_80,code=sm_80”
]



7. **PRECISION REQUIREMENTS**: Ensure exact mathematical alignment:
   - Use same mathematical formulas as PyTorch implementation
   - Verify with torch.allclose(rtol=1e-03, atol=1e-6)
   - Test with different input sizes

Here's the target architecture to optimize:

python
import torch
import torch.nn as nn

class Model(nn.Module):
“”"
Huber Loss implementation - robust loss function that is quadratic for small errors and linear for large errors.
“”"
def init(self, delta=1.0):
super(Model, self).init()
self.delta = delta

def forward(self, input: torch.Tensor, target: torch.Tensor) -> torch.Tensor:  
    """  
    Compute Huber loss between input and target tensors.

    Args:  
        input (torch.Tensor): Predicted values [batch_size, ...]  
        target (torch.Tensor): Target values [batch_size, ...]  

    Returns:  
        torch.Tensor: Scalar Huber loss value
    """  
    # Compute absolute error
    abs_error = torch.abs(input - target)
    
    # Create masks for quadratic and linear regions
    quadratic_mask = abs_error <= self.delta
    linear_mask = ~quadratic_mask
    
    # Quadratic region: 0.5 * error²
    quadratic_loss = 0.5 * torch.pow(input - target, 2)
    
    # Linear region: delta * (|error| - 0.5 * delta)
    linear_loss = self.delta * (abs_error - 0.5 * self.delta)
    
    # Combine losses
    loss = torch.where(quadratic_mask, quadratic_loss, linear_loss)
    
    return torch.sum(loss)
batch_size = 256
feature_dim = 512
delta = 1.0

def get_inputs():
input = torch.randn(batch_size, feature_dim)
target = torch.randn(batch_size, feature_dim)
return [input, target]

def get_init_inputs():
return [delta] # delta parameter



**EXPECTED OUTPUT STRUCTURE**:
Generate two files:
1. `huberloss_cudacode.py` - Contains ModelNew class with two-level reduction approach
2. `huberloss_torchcode.py` - Contains the reference PyTorch implementation

**KEY REQUIREMENTS**:
- The CUDA implementation must use two-level reduction to avoid atomic operation bottlenecks
- Must implement the exact HuberLoss formula: 0.5*error² when |error| ≤ delta, delta*(|error| - 0.5*delta) otherwise
- Must use shared memory for block-level reduction
- Must use dynamic block allocation based on input size
- Must complete final reduction with torch.sum() on GPU
- Must maintain mathematical precision with PyTorch implementation
- Must handle arbitrary tensor shapes (not just 2D)