I am trying to use scalapack functions (pdgetri, pdgetrf) in my cython code. This is my code:
scalapack_funcs.pyx
- Code: Select all
import numpy as np
cimport numpy as np
cimport cython
cimport openmp
from cpython.pythread cimport *
cdef extern from "scalapack_func.h" nogil:
void C_PDGETRI( long* N, double* A, int* IA, int* JA, int* DESCA, int* IPIV, double* WORK, long* LWORK, int* IWORK, long* LIWORK, int* INFO )
void C_PDGETRF( long* M, long* N, double* A, int* IA, int* JA, int* DESCA, int* IPIV, int* INFO)
cdef double[:,:] pdgetri_(double[:,:] A):
cdef long N = A.shape[0]
cdef long LWORK = A.shape[1]
cdef long LIWORK = A.shape[1]
cdef int INFO = 0
cdef int IA = 0 #the row index in the global array A indicating the first row of sub( A )
cdef int JA = 0 #The column index in the global array A indicating the first column of sub( A ).
cdef double[::1] WORK = np.empty(LWORK, dtype=np.float64)
cdef int[::1] IWORK = np.empty(LIWORK, dtype=np.int32)
cdef int[::1] IPIV = np.empty(N, dtype=np.int32)
cdef int[::1] DESCA = np.empty(N, dtype=np.int32)
C_PDGETRI(&N, &A[0,0], &IA, &JA, &DESCA[0], &IPIV[0], &WORK[0], &LWORK, &IWORK[0], &LIWORK, &INFO)
if INFO < 0:
raise ValueError('illegal value in %d-th argument of internal pdgetri' % -INFO)
elif INFO > 0:
raise ValueError("Diagonal number %d is exactly zero. Singular matrix." % -INFO)
return A
cdef double[:,:] pdgetrf_(double[:,:] A):
cdef long M = A.shape[0] #The number of rows to be operated on,
cdef long N = A.shape[1] #The number of columns to be operated on
cdef int INFO = 0
cdef int IA = 0 #the row index in the global array A indicating the first row of sub( A )
cdef int JA = 0 #The column index in the global array A indicating the first column of sub( A ).
cdef int[::1] IPIV = np.empty(M, dtype=np.int32)
cdef int[::1] DESCA = np.empty(M, dtype=np.int32)
C_PDGETRF(&M, &N, &A[0,0], &IA, &JA, &DESCA[0], &IPIV[0], &INFO)
if INFO < 0:
raise ValueError('illegal value in %d-th argument of internal pdgetrf' % -INFO)
elif INFO > 0:
raise ValueError("Diagonal number %d is exactly zero. Singular matrix." % -INFO)
return A
cpdef double[:,:] inverse(double[:,:] A):
if (A.shape[0] !=A.shape[1]):
raise TypeError("The input array must be a square matrix")
if not np.ascontiguousarray(A, dtype=np.float64).flags['F_CONTIGUOUS']:
A = A.copy_fortran()
pdgetrf_(A)
pdgetri_(A)
return A
scalapack_wrap.f90
- Code: Select all
subroutine C_PDGETRI( N, A, IA, JA, DESCA, IPIV, WORK, LWORK, IWORK, LIWORK, INFO ) bind(c, NAME='C_PDGETRI')
use iso_c_binding, only: c_double, c_int
implicit none
real(c_double), dimension(N, N), intent(in):: A
real(c_double), dimension(LWORK), intent(in):: WORK
integer(c_int), dimension(N), intent(in):: DESCA
integer(c_int), dimension(N), intent(in)::IPIV
integer(c_int), dimension(LWORK), intent(in):: IWORK
integer(c_int), intent(in) :: N, IA, JA, LIWORK, LWORK, INFO
call PDGETRI( N, A, IA, JA, DESCA, IPIV, WORK, LWORK, IWORK, LIWORK, INFO )
end subroutine C_PDGETRI
subroutine C_PDGETRF(M, N, A, IA, JA, DESCA, IPIV, INFO) bind(c, NAME='C_PDGETRF')
use iso_c_binding, only: c_double, c_int
implicit none
real(c_double), dimension(M, N), intent(in):: A
integer(c_int), dimension(N), intent(in):: DESCA
integer(c_int), dimension(N), intent(in)::IPIV
integer(c_int), intent(in) :: M, N, IA, JA, INFO
call PDGETRF( M, N, A, IA, JA, DESCA, IPIV, INFO )
end subroutine C_PDGETRF
scalapack_func.h
- Code: Select all
extern void C_PDGETRI( long* N, double* A, int* IA, int* JA, int* DESCA, int* IPIV, double* WORK, long* LWORK, int* IWORK, long* LIWORK, int* INFO );
extern void C_PDGETRF( long* M, long* N, double* A, int* IA, int* JA, int* DESCA, int* IPIV, int* INFO);
And finally the "setup.py" code
- Code: Select all
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
from numpy import get_include
import subprocess
import os
from os import system
import warnings
import mpi4py
os.environ["CC"] = "mpicc"
fortran_mod_comp = 'mpif90 /usr/local/scalapack-2.0.2/SRC/pdgetri.f -c -o pdgetri_.o -O3 -fPIC'
print fortran_mod_comp
system(fortran_mod_comp)
fortran_mod_comp = 'mpif90 /usr/local/scalapack-2.0.2/SRC/pdgetrf.f -c -o pdgetrf_.o -O3 -fPIC'
print fortran_mod_comp
system(fortran_mod_comp)
shared_obj_comp = 'mpif90 scalapack_wrap.f90 -c -o scalapack_wrap.o -O3 -fPIC'
print shared_obj_comp
system(shared_obj_comp)
def runcommand(cmd):
process = subprocess.Popen(cmd.split(), shell=False, stdout=subprocess.PIPE,
stderr=subprocess.STDOUT, universal_newlines=True)
c = process.communicate()
if process.returncode != 0:
raise Exception("Something went wrong whilst running the command: %s" % cmd)
return c[0]
def whichscalapack():
if 'MKLROOT' in os.environ:
return 'intelmkl'
else:
return 'netlib'
def whichmpi():
import re
try:
mpiv = runcommand('mpirun -V')
if re.search('Intel', mpiv):
return 'intelmpi'
elif re.search('Open MPI', mpiv):
return 'openmpi'
except:
return 'mpich'
warnings.warn('Unknown MPI environment.')
return None
scalapackversion = whichscalapack()
mpiversion = whichmpi()
if mpiversion == 'openmpi':
mpilinkargs = runcommand('mpicc -showme:link').split()
mpicompileargs = runcommand('mpicc -showme:compile').split()
if scalapackversion == 'intelmkl':
scl_lib = ['mkl_scalapack_lp64', 'mkl_intel_lp64','mkl_core', 'mkl_sequential', 'mkl_blacs_'+mpiversion+'_lp64', 'iomp5', 'pthread','mkl_intel_thread',"openblas", "blas", "gfortran"]
scl_libdir = [os.environ['MKLROOT']+'/lib/intel64' if 'MKLROOT' in os.environ else '']
elif scalapackversion == 'netlib':
scl_lib = ['scalapack', 'gfortran']
scl_libdir = [ os.path.dirname(runcommand('gfortran -print-file-name=libgfortran.a')) ]
else:
raise Exception("Scalapack distribution unsupported. Please modify setup.py manually.")
ext_modules = [Extension(# module name:
'scalapack_funcs',
sources=['scalapack_funcs.pyx'],
include_dirs = [get_include(), mpi4py.get_include(),os.environ['MKLROOT']+"/include"],
library_dirs=scl_libdir, libraries=scl_lib,
extra_compile_args=[ "-fopenmp" ] + mpicompileargs,
extra_link_args=['pdgetrf_.o', 'pdgetri_.o', 'scalapack_wrap.o',"-fopenmp"] + mpilinkargs)]
setup(
cmdclass = {'build_ext': build_ext},
ext_modules = ext_modules
)
The code gets compiled but when I import it in python I get the following error:
- Code: Select all
>>> import scalapack_funcs
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ImportError: /opt/intel/mkl/lib/intel64/libmkl_scalapack_lp64.so: undefined symbol: dapply_2hv
I could not find anythin related to dapply_2hv in the web. Why do I get this error? Could anybody explain it? Thanks in advance.
