Repository: pytorch
Issue: Numerical-reproducibility issue in torch.matmul
Labels: module: cuda,triaged,module: numerical-reproducibility

## 🐛 Bug

<!-- A clear and concise description of what the bug is. -->
I tried to utilize the boardcast feature in torch.matmul, but here are some numerical errors when involves certain types of sparse matrix. To be clear, I did not use any sparse matrix format but most elements in the matrix ```lapa_c``` are zero (see the sample code I provide below).  

## To Reproduce
<!-- If you have a code sample, error messages, stack traces, please provide it here as well -->
```Python
import torch

def Build_CoefficientMatrix(c:list,meshx:int,periodic:bool=False):
    '''
    c is list of FD coefficient
    e.g. for 1st derivative with 2nd accuracy central difference:
    c=[-0.5,0] 
    '''
    if 2*len(c)-1>=meshx: raise ValueError
    acc = len(c)   
    
    tmp=[]
    c.reverse()
    for i in range(acc):
        x = torch.cat((torch.cat((torch.zeros((i,meshx-i)),
                                    c[i]*torch.eye(meshx-i)),dim=0),
                                    torch.zeros((meshx,i))
                                    ),dim=1)
        tmp.append(x)
    re=tmp[0]
    for k in tmp[1:]:
        re+=k+k.T

    if periodic:
        re[:acc,-acc:]=re[acc:2*acc,:acc]
        re[-acc:,:acc]=re[:acc,acc:2*acc]
    return re

lapa_c=Build_CoefficientMatrix([1/90,-3/20,3/2,-49/18],64,periodic=True)

x=torch.rand(2,64,1)
out1=torch.matmul(lapa_c,x)
out2=torch.matmul(lapa_c,x[0,:,0])
print(out1[0,0,0].item())
print(out2[0].item())
```
The ```Build_CoefficientMatrix``` is used to build a [band matrix](https://en.wikipedia.org/wiki/Band_matrix) (with some non-zero elements in the bottom left and upper right corner when ```periodic=True```)


## Expected behavior
The two printed outputs are supposed to be the same, but actually they are not. 

##
if you change: 
* the batch size of ```x``` back to ```1```(```x=torch.rand(1,64,1)```)
* the ```periodic``` parameter in ```Build_CoefficientMatrix``` to ```False```
* the ```c``` parameter in ```Build_CoefficientMatrix``` to a list with length no greater than 3

the two results match each other.


## Environment

PyTorch version: 1.8.1
Is debug build: False
CUDA used to build PyTorch: 11.1
ROCM used to build PyTorch: N/A

OS: Ubuntu 20.04.2 LTS (x86_64)
GCC version: (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Clang version: Could not collect
CMake version: version 3.16.3
Libc version: glibc-2.31

Python version: 3.8 (64-bit runtime)
Python platform: Linux-5.8.0-55-generic-x86_64-with-glibc2.10
Is CUDA available: True
CUDA runtime version: Could not collect
GPU models and configuration: GPU 0: GeForce RTX 3070
Nvidia driver version: 460.80
cuDNN version: Could not collect
HIP runtime version: N/A
MIOpen runtime version: N/A

Versions of relevant libraries:
[pip3] numpy==1.19.2
[pip3] numpydoc==1.1.0
[pip3] torch==1.8.1
[pip3] torch-cluster==1.5.8
[pip3] torch-geometric==1.6.3
[pip3] torch-scatter==2.0.5
[pip3] torch-sparse==0.6.8
[pip3] torch-spline-conv==1.2.0
[pip3] torchaudio==0.8.0a0+e4e171a
[pip3] torchvision==0.9.1
[conda] blas                      1.0                         mkl  
[conda] cudatoolkit               11.1.1               h6406543_8    conda-forge
[conda] ffmpeg                    4.3                  hf484d3e_0    pytorch
[conda] mkl                       2020.4             h726a3e6_304    conda-forge
[conda] mkl-service               2.3.0            py38h1e0a361_2    conda-forge
[conda] mkl_fft                   1.3.0            py38h5c078b8_1    conda-forge
[conda] mkl_random                1.2.0            py38hc5bc63f_1    conda-forge
[conda] numpy                     1.19.2           py38h54aff64_0  
[conda] numpy-base                1.19.2           py38hfa32c7d_0  
[conda] numpydoc                  1.1.0                      py_1    conda-forge
[conda] pytorch                   1.8.1           py3.8_cuda11.1_cudnn8.0.5_0    pytorch
[conda] torch-cluster             1.5.8                    pypi_0    pypi
[conda] torch-geometric           1.6.3                    pypi_0    pypi
[conda] torch-scatter             2.0.5                    pypi_0    pypi
[conda] torch-sparse              0.6.8                    pypi_0    pypi
[conda] torch-spline-conv         1.2.0                    pypi_0    pypi
[conda] torchaudio                0.8.1                      py38    pytorch
[conda] torchvision               0.9.1                py38_cu111    pytorch

## Additional context

<!-- Add any other context about the problem here. -->


cc @VitalyFedyunin @ngimel @zasdfgbnm @ptrblck