Wrapping magma_dpotrf for python

[psims01]

Hi,

I am trying to perform an out-of-memory Cholesky factorisation of a real positive definite matrix on a GPU from python. magma_dpotrf is a fantastic solution but I'm having trouble accessing it. This is an attempt to access it using ctypes (a derivative of an attempt to achieve something similar I found browsing for solutions online).

Code: Select all

import pycuda.driver as cuda
import pycuda.autoinit
import numpy as np
import ctypes

magma_path = '/path/to/libmagma.so'
libmagma = ctypes.cdll.LoadLibrary(magma_path)

int_pointer = lambda x: ctypes.pointer(ctypes.c_int(x))
char_pointer = lambda x: ctypes.pointer(ctypes.c_char(x))

_magma_dpotrf = libmagma.magma_spotrf_gpu     
_magma_dpotrf.restype = ctypes.c_uint
_magma_dpotrf.argtypes = [ctypes.c_char_p,
                               ctypes.POINTER(ctypes.c_int),
                               ctypes.POINTER(ctypes.c_double),
                               ctypes.POINTER(ctypes.c_int),
                               ctypes.POINTER(ctypes.c_int)]

#magma_get_spotrf_nb = libmagma.magma_get_spotrf_nb
#magma_get_spotrf_nb.restype = ctypes.c_uint
#magma_get_spotrf_nb.argtypes = [ctypes.POINTER(ctypes.c_int)]


def magma_dpotrf(uplo, N, h_R, lda, info):
   _magma_dpotrf(char_pointer(uplo),
                       int_pointer(N),
                       (ctypes.c_double * len(h_R))(*h_R),
                       int_pointer(lda),
                       int_pointer(info))
   return info    

if __name__ == '__main__':
   N = 5
   lda=N
   
   # Create matrix to be factored
   h_R = (np.eye(N) + np.ones((N,N))*0.5).astype('double').flatten()

   libmagma.magma_init()
   info = 1

   # # Do Cholesky factorization
   info = magma_dpotrf('L', N, h_R, lda, info)



magma_spotrf_gpu(char *uplo, int *n, float *a, int *lda,
                     float *work, int *info)

magma_dpotrf( opts.uplo, N, h_R, lda, &info );

Everything is fine until (the key step) the final function call

Code: Select all

_magma_dpotrf(char_pointer(uplo), int_pointer(N), (ctypes.c_double * len(h_R))(*h_R), int_pointer(lda), int_pointer(info))

which segfaults and gives,

python: constants.cpp:357: const char* lapack_uplo_const(magma_uplo_t): Assertion `magma_const <= MagmaFull' failed.
Aborted

I'm using MAGMA 1.6.2, cuda 6.5, python 2.7 and pycuda 2014.1.

Could anyone familiar with Magma and python help me with how to correct the issue(s) causing the segmentation fault?
Alternatively, I'm not wedded to ctypes for wrapping Magma (I initially tried SWIG but without success), if anyone could point me to a working example of wrapping Magma functions for Python that would also be much appreciated!

Thank you in advance,

Peter

[mgates3]

Two points:

1) In MAGMA, uplo is an enum, not a char. You can either use the constants MagmaLower and MagmaUpper, or convert characters 'L' and 'U' using magma_uplo_const( ).

Code: Select all

magma_dpotrf( MagmaLower, N, A, lda, &info );
magma_dpotrf( magma_uplo_const('L'), N, A, lda, &info );

2) In case you didn't see, there are two functions:
magma_dpotrf takes A in CPU host memory. It is a drop-in replacement for LAPACK's dpotrf.
magma_dpotrf_gpu takes dA in GPU device memory. The usual convention in MAGMA is to prefix arrays with "d" when they are in GPU device memory.

-mark

[psims01]

Thanks Mark.

In reply to 2.) I'm using an NVIDIA K20 GPU so I've got 5 GB of memory. The matrix I want to factorise is ~10 GB in size so I don't think I can use magma_dpotrf_gpu because I can't allocate that memory. I thought magma_dpotrf moved parts of the matrix on and off the GPU in blocks and was for the out-of-GPU-memory case?

I'll have a look at 1.) next week and post the result.

[mgates3]

Yes, magma_dpotrf should work for the out-of-GPU-memory case. The A matrix is in the CPU host memory. Internally, parts of it are put into GPU memory. I just wanted to clarify since it wasn't clear where your matrix resided, and calling the wrong routine is a common source of segfaults.

-mark

[psims01]

Hi Mark,

Looking at your point 1.) I've switched, char_pointer(uplo) with uplo='L' to libmagma.magma_uplo_const(ctypes.c_char(uplo)) which is equal to 122. That improves the situation and there is no longer a segmentation fault. However I'm still not getting back the factorised matrix and magma_dpotrf returns with info =-113. The documentation has '< 0: if INFO = -i, the i-th argument had an illegal value or another error occured, such as memory allocation failed. ' so this falls into the 'another error' category. I'm not certain where to go from here.

Peter

[mgates3]

Errors are discussed in the documentation:
http://icl.cs.utk.edu/projectsfiles/mag ... rrors.html

magma_strerror( err ) will provide a description of the error code. It's best to use the symbolic names, but if needed you can look up the numeric values of constants in the header file.

Code: Select all

magma>  grep ERR include/magma_types.h 
// LAPACK argument errors are < 0 but > MAGMA_ERR.
// MAGMA errors are < MAGMA_ERR.
#define MAGMA_ERR                  -100
#define MAGMA_ERR_NOT_INITIALIZED  -101
#define MAGMA_ERR_REINITIALIZED    -102
#define MAGMA_ERR_NOT_SUPPORTED    -103
#define MAGMA_ERR_ILLEGAL_VALUE    -104
#define MAGMA_ERR_NOT_FOUND        -105
#define MAGMA_ERR_ALLOCATION       -106
#define MAGMA_ERR_INTERNAL_LIMIT   -107
#define MAGMA_ERR_UNALLOCATED      -108
#define MAGMA_ERR_FILESYSTEM       -109
#define MAGMA_ERR_UNEXPECTED       -110
#define MAGMA_ERR_SEQUENCE_FLUSHED -111
#define MAGMA_ERR_HOST_ALLOC       -112
#define MAGMA_ERR_DEVICE_ALLOC     -113
#define MAGMA_ERR_CUDASTREAM       -114
#define MAGMA_ERR_INVALID_PTR      -115
#define MAGMA_ERR_UNKNOWN          -116
#define MAGMA_ERR_NOT_IMPLEMENTED  -117

It looks like you have a device memory allocation failure (-113). This is a little odd, because the out-of-GPU-memory algorithm should kick in. See results below on my machine, where N=15000 fails for GPU interface, but succeeds for CPU interface with out-of-GPU-memory algorithm. What size problem are you attempting? What GPU are you using? Can you run the problem using the MAGMA tester (testing/testing_dpotrf and testing/testing_dpotrf_gpu) to see its results, and post the complete input & output.

Code: Select all

magma-trunk/testing> ./testing_dpotrf_gpu --range 14000:17000:1000 
% MAGMA 1.6.2 svn compiled for CUDA capability >= 3.0
% CUDA runtime 5050, driver 6050. MAGMA not compiled with OpenMP. MKL 11.2.3, MKL threads 4. 
% device 0: GeForce GT 750M, 925.5 MHz clock, 2047.6 MB memory, capability 3.0
% Mon Jun 15 12:24:35 2015
% Usage: ./testing_dpotrf_gpu [options] [-h|--help]

% uplo = Lower
% N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R_magma - R_lapack||_F / ||R_lapack||_F
%=======================================================
14000     ---   (  ---  )     26.46 (  34.57)     ---  
!!!! magma_malloc failed for: d_A

magma-trunk/testing> ./testing_dpotrf --range 14000:17000:1000
% MAGMA 1.6.2 svn compiled for CUDA capability >= 3.0
% CUDA runtime 5050, driver 6050. MAGMA not compiled with OpenMP. MKL 11.2.3, MKL threads 4. 
% device 0: GeForce GT 750M, 925.5 MHz clock, 2047.6 MB memory, capability 3.0
% Mon Jun 15 12:25:55 2015
% Usage: ./testing_dpotrf [options] [-h|--help]

% ngpu = 1, uplo = Lower
%   N   CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R_magma - R_lapack||_F / ||R_lapack||_F
%=======================================================
14000     ---   (  ---  )     26.17 (  34.95)     ---  
15000     ---   (  ---  )     26.21 (  42.93)     ---  
16000     ---   (  ---  )     26.27 (  51.97)     ---  
17000     ---   (  ---  )     26.23 (  62.45)     ---  

[psims01]

Hi Mark,

Thanks, that's good to know. The GPU in an NVIDIA K20 and the matrix is N=5 so size is definitely not an issue here, I'm simply attempting to test the wrapper. I just left the lab but I'll run and post the testing result once I'm back.

Peter

[psims01]

Here are the testing_dpotrf and testing_dpotrf_gpu outputs I get (currently running on a cluster login node which is where I am testing the python wrapper but not the GPUs I'll run the final product):

Code: Select all

 ../testing/testing_dpotrf
MAGMA 1.6.2  compiled for CUDA capability >= 2.0 
CUDA runtime 6000, driver 6050. OpenMP threads 1. MKL 10.3.10, MKL threads 1. 
ndevices 8                                                                    
device 0: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 1: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 2: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 3: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 4: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 5: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 6: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
device 7: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0          
Usage: ../testing/testing_dpotrf [options] [-h|--help]                        

ngpu = 1, uplo = Lower
    N   CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R_magma - R_lapack||_F / ||R_lapack||_F
========================================================                                 
 1088     ---   (  ---  )      4.43 (   0.10)     ---                                    
 2112     ---   (  ---  )     10.03 (   0.31)     ---                                    
 3136     ---   (  ---  )     10.93 (   0.94)     ---                                    
 4160     ---   (  ---  )     11.30 (   2.12)     ---                                    
 5184     ---   (  ---  )     11.53 (   4.03)     ---                                    
 6208     ---   (  ---  )     11.71 (   6.81)     ---                                    
 7232     ---   (  ---  )     11.84 (  10.65)     ---                                    
 8256     ---   (  ---  )     11.94 (  15.72)     ---                                    
 9280     ---   (  ---  )     12.03 (  22.16)     ---
10304     ---   (  ---  )     12.09 (  30.16)     ---

and,

Code: Select all

../testing/testing_dpotrf_gpu
MAGMA 1.6.2  compiled for CUDA capability >= 2.0
CUDA runtime 6000, driver 6050. OpenMP threads 1. MKL 10.3.10, MKL threads 1.
ndevices 8
device 0: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 1: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 2: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 3: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 4: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 5: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 6: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
device 7: GRID K1, 849.5 MHz clock, 4095.8 MB memory, capability 3.0
Usage: ../testing/testing_dpotrf_gpu [options] [-h|--help]

uplo = Lower
  N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R_magma - R_lapack||_F / ||R_lapack||_F
========================================================
 1088     ---   (  ---  )      6.13 (   0.07)     ---
 2112     ---   (  ---  )     10.10 (   0.31)     ---
 3136     ---   (  ---  )     11.02 (   0.93)     ---
 4160     ---   (  ---  )     11.36 (   2.11)     ---
 5184     ---   (  ---  )     11.58 (   4.01)     ---
 6208     ---   (  ---  )     11.75 (   6.79)     ---
 7232     ---   (  ---  )     11.87 (  10.62)     ---
 8256     ---   (  ---  )     11.97 (  15.67)     ---
 9280     ---   (  ---  )     12.05 (  22.11)     ---
10304     ---   (  ---  )     12.12 (  30.08)     ---

[mgates3]

It looks like MAGMA is working, though the performance is quite low.

Two possible reasons for the low performance are:
1) Using 1 thread for MKL. Try changing $OMP_NUM_THREADS or $MKL_NUM_THREADS to something higher.
2) It appears that the K1 GPU is primarily a graphics GPU. I can't find the double precision specs for that card, but it may be significantly slower than single precision (like 10x slower). Try running testing_sgemm and testing_dgemm to see the difference.

Did you resolve the info = -113 (device malloc failed) issue? If not, can you post your make.inc file and how you compile your own application? Since you said it was for a small N, it's possible there's something weird going on with how things are compiled and linked together.

-mark

[psims01]

Yes, I just remade the latest version of MAGMA and now magma_dpotrf is giving what I expect! Thanks for your help! I now have a related question: the problem I would actually like the solution to is, x, for Ax=b with N~30000. Now that I have the Cholesky decomposition of A I can use lapack pdotrs on the CPU for x but are there any (out-of-gpu-memory) gpu alternatives for this step?

Peter