MVP python BLAS example

[akochepasov]

Hi guys, I'm trying to put together a simple MVP BLAS example for a C-ordered array, according to https://docs.nvidia.com/cuda/nvmath-python/latest/bindings/index.html
So far the result is mixed.

import numpy as np
import nvmath

n, m = 503, 91
A = np.random.random((m, n))
C = np.zeros((n, n))

alpha, beta = 1.0, 0.0
uplo, trans = b'L', b'N'

k = A.shape[0]
lda = A.shape[1]
ldc = C.shape[1]

# Dealing with C-order
uplo_  = 1 if  uplo != b'U' else 0
trans_ = 0 if trans != b'N' else 1

a = np.array(alpha, dtype=np.float64)
b = np.array(beta, dtype=np.float64)

hndl = nvmath.bindings.cublas.create()
nvmath.bindings.cublas.dsyrk(hndl, uplo_, trans_, n, k, 
    a.ctypes.data, A.ctypes.data, lda, 
    b.ctypes.data, C.ctypes.data, ldc)
nvmath.bindings.cublas.destroy(hndl)

This example doesn't return result in C and crashes in cublas.create() on the second call:

cuBLASError: the library was not initialized (CUBLAS_STATUS_NOT_INITIALIZED)

Can you help me to straighten this out? How can I make C keep the result and run multiple times?

Configuration:

• Windows 11
• nvmath.bindings.cublas.get_cudart_version() returns 12030
• nvmath version is 0.1.0 installed with pip

[szkarpinski]

Hello @akochepasov ! Thank you for your interest in nvmath-python.

The problem that I see in your code is that cuBLAS expect the matrices to be on the GPU, and passing host-side pointers might result in undefined behaviour. I guess the error you are seeing might be caused by that. Can you please try to convert your A and C numpy arrays to CuPy arrays and check if this fixes the error?

CuPy has .random and .zeros, so you just need to change np to cp when creating the arrays:

import cupy as cp

A = cp.random.random((m, n))
C = cp.zeros((n, n))

And to get the pointers to CuPy, you can use .data.ptr instead of ctypes.data:

nvmath.bindings.cublas.dsyrk(hndl, uplo_, trans_, n, k, 
    a.ctypes.data, A.data.ptr, lda, 
    b.ctypes.data, C.data.ptr, ldc)

One more thing is that you should only destroy cublas handle after synchronizing with the GPU. Currently, you're destroying it right after the operation is scheduled, but you should wait for the GPU to finish.

Please let me know if the changes above help with the the library was not initialized error.

On our long-term roadmap, we have the introduction a Pythonic API that will cover syrk and similar cuBLAS calls, hiding their complexity just as nvmath.linalg.advanced.matmul hides the complexity of cuBLAS GEMM operations. Would that be useful for you? What other linear algebra operations would you like to see in nvmath-python?

[akochepasov]

Hey @szkarpinski , thank you for the reply. For sure, having explicit calls to transfer and sync the data over GPU boundaries works. This is what cupy has been for for years. In this case, nvmath not really needed unless I work on the native level (pxd cimport and cython.compiled). cupy doesn't provide (and will not) pxd headers.

While asking the initial question I was based on the nvmath documentation, which tells:

# passing a host buffer of int type can be done like this
buf = numpy.array([0, 1, 3, 5, 6], dtype=numpy.int32)
my_func(..., buf.ctypes.data, ...)

and

# create a host buffer to hold 5 int
buf = numpy.empty((5,), dtype=numpy.int32)
# pass buf's pointer to the wrapper
# buf could get modified in-place if the function writes to it
my_func(..., buf.ctypes.data, ...)
# examine/use buf's data
print(buf)

(from https://docs.nvidia.com/cuda/nvmath-python/latest/bindings/index.html)

I believe the documentation can be improved in this aspect as now it's unclear that the data must be transferred as in a regular CUDA model.
Speaking about operations, I need axpy scal sbmv ger syrk and gem. Using BLAS speeds up computation at least twice and it's really hard to beat those calls in own C++/CUDA.
Besides BLAS, I'm trying to eliminate any python calls under the hood as having GIL affects the performance a lot. The resulting native c code can be seen with cython -a call and having as least yellow and as much white as possible is a performance key.

[szalpal]

@akochepasov ,

I'm curious about the GIL issue you've mentioned. Was your point that nvmath-python acquires GIL and blocks other calls from executing, or were you referring to general Python calls from other libraries?

[akochepasov]

I don't have a reliable information regarding GPU yet as I'm still currently based on scipy blas exports. My current performance measurement shows that getting rid of the python code in the calls, whether it's python structs dereferencing, GIL calls, __Pyx_INCREF, __Pyx_DECREF speeds up my computational process 50-200%, depending on the numbers of calls and input data size.
As a rule of thumb, cimport and cdef are usually noticeably ahead of import and cpdef

[szkarpinski]

@akochepasov , we'll work on improving the docs of the bindings module to make it clear that the data must be explicitly transferred to the GPU. I agree that the current docs might be confusing in this aspect.