Repository: pytorch
Issue: Large cumulative sums appear to be nondeterministic. 
Labels: high priority,module: cuda,triaged,module: determinism

### 🐛 Describe the bug

@timothytmzhou and I found that cumulative sums over large input arrays do not seem to deterministically produce the same output array:

```python
import torch
import math
import random
import numpy as np
import os

# make ourselves as deterministic as possible.
os.environ["CUBLAS_WORKSPACE_CONFIG"] = ":4096:8"
np.random.seed(0)
torch.manual_seed(0)
random.seed(0)
torch.use_deterministic_algorithms(True)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False

# (deterministic) messy number generator
v = torch.Tensor([(x*math.pi)%1 for x in range(1000000)]).to(torch.device("cuda"))

# init a and b to the same cumulative sum.
a = b = v.cumsum(0)
counter = 0

# this loop should never terminate.
while ((a==b).all()):
	b = v.cumsum(0)
	counter+=1

print(f"non-determinism after {counter} runs.")
print("errors:",(a-b).abs()[a-b!=0])
print("first error index:",torch.where((a-b)!=0)[0][0])
```
yields, on my machine:
```
non-determinism after 1 runs.
errors: tensor([0.0020, 0.0020, 0.0020,  ..., 0.2188, 0.2188, 0.2188], device='cuda:0')
first error index: tensor(38400, device='cuda:0')
```
Multiple runs of this file exhibit unique non-determinism. Here's another run:
```
non-determinism after 1 runs.
errors: tensor([0.0039, 0.0039, 0.0039,  ..., 0.0938, 0.0938, 0.0938], device='cuda:0')
first error index: tensor(76800, device='cuda:0')
```
We decided to test this after reading [this paper](https://www.cs.purdue.edu/homes/lintan/publications/variance-ase20.pdf). Of note is the following paragraph:

> Scheduling and floating-point precision: GPU programming uses warp as a unit of computation. A warp consists of 32 parallel threads with concurrent memory access. Due to the limited precision (32-bits) of DL models, rounding errors are introduced at every step of floating-point calculation. In the GPU model, concurrent accesses to a single memory address must be serialized to prevent race conditions with no guaranteed order of access. Therefore, the rounding error introduced in each warp may vary across fixed-seed identical training runs due to the different accessing orders.

### Versions

Collecting environment information...
PyTorch version: 1.11.0
Is debug build: False
CUDA used to build PyTorch: 11.6
ROCM used to build PyTorch: N/A

OS: Arch Linux (x86_64)
GCC version: (GCC) 11.2.0
Clang version: 13.0.1
CMake version: version 3.22.3
Libc version: glibc-2.35

Python version: 3.10.4 (main, Mar 23 2022, 23:05:40) [GCC 11.2.0] (64-bit runtime)
Python platform: Linux-5.16.16-arch1-1-x86_64-with-glibc2.35
Is CUDA available: True
CUDA runtime version: 11.6.112
GPU models and configuration: GPU 0: NVIDIA GeForce GTX 1050 Ti
Nvidia driver version: 510.54
cuDNN version: Probably one of the following:
/usr/lib/libcudnn.so.8.3.2
/usr/lib/libcudnn_adv_infer.so.8.3.2
/usr/lib/libcudnn_adv_train.so.8.3.2
/usr/lib/libcudnn_cnn_infer.so.8.3.2
/usr/lib/libcudnn_cnn_train.so.8.3.2
/usr/lib/libcudnn_ops_infer.so.8.3.2
/usr/lib/libcudnn_ops_train.so.8.3.2
HIP runtime version: N/A
MIOpen runtime version: N/A
Is XNNPACK available: True

Versions of relevant libraries:
[pip3] mypy-extensions==0.4.3
[pip3] numpy==1.22.3
[pip3] torch==1.11.0
[conda] Could not collect

cc @ezyang @gchanan @zou3519 @ngimel @mruberry @kurtamohler