CODEWITHSUNDEEP
pythonscipyeigenvalueintel-mkl
Tunneller
Tunneller

Reputation: 589

Calling MKL from Python : DSTEVR

dstevr computes eigensolutions of triangular symmetric matrix. Cool. Except it was not one of the routines ported with a wrapper to SCIPY. So I've followed the instructions on how to call MKL directly from Python and the attached seems to give the correct answer. But gosh.... Is there someway to clean this up?!

  import numpy as np 
  from scipy import linalg
  from ctypes import * 

  c_double_p = POINTER(c_double)
  c_int_p = POINTER(c_int)
  c_char_p = POINTER(c_char)

  mkl = CDLL('mkl_rt.dll')
  dstevr = mkl.dstevr

  #SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, M, 
  # *      W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO)
  #  CHARACTER * 1 JOBZ, RANGE
  #  INTEGER N, IL, IU, M, LDZ, LWORK, LIWORK, INFO
  #  INTEGER ISUPPZ(*), IWORK(*)
  #  DOUBLE PRECISION VL, VU, ABSTOL
  #  DOUBLE PRECISION D(*), E(*), W(*), Z(LDZ,*), WORK(*)

  dstevr.argtypes = [ c_char_p, c_char_p, c_int_p, c_double_p, c_double_p,  c_double_p, c_double_p, c_int_p, c_int_p, c_double_p,  \
 c_int_p, c_double_p, c_double_p, c_int_p, c_int_p, c_double_p, c_int_p, c_int_p, c_int_p, c_int_p]


  sv = "v"
  eig_j = c_char_p(c_char(sv[0]))
  sr = "a"
  eig_r = c_char_p(c_char(sr[0]))

  vl = c_double(0.0)
  vu = c_double(0.0)
  il = c_int(0)
  iu = c_int(0)
  abstol = c_double(0.0)
  cm = c_int(0)
  eig_info = c_int(0)

  N = 6
  cn = c_int(N)
  ldz = cn 

  lwork = c_int(20*N)
  liwork = c_int(10*N)

  diag = np.ascontiguousarray(np.ones(N)*2)
  diag_p = diag.ctypes.data_as(c_double_p)

  offdiag = np.ascontiguousarray(np.ones(N-1)*(-1))
  offdiag_p = offdiag.ctypes.data_as(c_double_p)

  isuppz = np.ascontiguousarray(np.ones(N*2),dtype=int)
  isuppz_p = isuppz.ctypes.data_as(c_int_p)

  eigw = np.ascontiguousarray(np.zeros(N))
  eigw_p = eigw.ctypes.data_as(c_double_p)

  workz = np.ascontiguousarray(np.ones(20*N))
  workz_p = workz.ctypes.data_as(c_double_p)

  iworkz = np.ascontiguousarray(np.ones(10*N),dtype=int)
  iworkz_p = iworkz.ctypes.data_as(c_int_p)

  eigz = np.ascontiguousarray(np.ones(N*N))
  eigz_p = eigz.ctypes.data_as(c_double_p)


  dstevr( eig_j, eig_r, byref(cn), diag_p, offdiag_p, byref(vl),byref(vu),byref(il),byref(iu),byref(abstol),byref(cm),\
  eigw_p, eigz_p, byref(ldz), isuppz_p, workz_p, byref(lwork), iworkz_p, byref(liwork), byref(eig_info)) 

  print "Eig_Info", eig_info

  print eigz

  A = np.eye(N,N,k=-1)*(-1) + np.eye(N,N)*2 + np.eye(N,N,k=1)*(-1)

  w,v= linalg.eigh(A)

  print v.T

Thanks,

Upvotes: 2

Views: 320

Answers (1)

ev-br
ev-br

Reputation: 26050

You can use instead a cython wrapper in cython_lapack, http://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html

You'll need cython though.

Upvotes: 0

Related Questions