lebedov / scikit-cuda

Python interface to GPU-powered libraries
http://scikit-cuda.readthedocs.org/
Other
986 stars 179 forks source link

linalg.eig fails for symmetric matrices larger than 2 #286

Open willw1 opened 4 years ago

willw1 commented 4 years ago

Problem

import pycuda.gpuarray as gpuarray import pycuda.autoinit import numpy as np from skcuda import linalg linalg.init()

Compute right eigenvectors of a symmetric matrix A and verify Avr = vrw

a = np.array(([1,3],[3,5]), np.float32, order='F') #--- Works

a = np.array(([1,2,3],[2,4,5],[3,5,6]), np.float32, order='F') #--- Fails - np.allclose() below gives False a_gpu = gpuarray.to_gpu(a) vr_gpu, w_gpu = linalg.eig(a_gpu, 'N', 'V') np.allclose(np.dot(a, vr_gpu.get()), np.dot(vr_gpu.get(), np.diag(w_gpu.get())), 1e-4)

Environment

ubuntu16.04 CUDA 8.0 Python 3.6.0 PyCUDA VERSION = (2019, 1, 2) installed with: pip install pycuda scikit-cuda version 0.5.3 installed with: pip install scikit-cuda

--EDIT-- Also fails with CUDA 10.1 on Ubuntu 18.04

--EDIT 2--

The function eig() should be returning the transpose of vr_gpu to agree with the documentation. As it stands, the eigenvectors are returned as rows in vr_gpu. The documentation (and test code) say/assume columns.

np.allclose(np.dot(a, vr_gpu.get().T), np.dot(vr_gpu.get().T, np.diag(w_gpu.get())), 1e-4,1e-4)

gives True in the above 3x3 case.