Calling MKL from Python : DSTEVR
<p><strong>dstevr</strong> 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?!</p>

<pre><code>  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
</code></pre>

<p>Thanks,</p>

<p>You can use instead a cython wrapper in cython_lapack, <a href="http://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html" rel="nofollow">http://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html</a></p>

<p>You'll need cython though.</p>

python scipy eigenvalue intel-mkl