PythonからMKLを呼び出す：DSTEVR

Question

dstevrは、三角対称行列の固有値を計算します。クール。それ以外は、SCIPYへのラッパーで移植されたルーチンの1つではありませんでした。だから、私はPythonからMKLを直接呼び出す方法の指示に従ってきました。添付されているのは正しい答えを与えるようです。しかし、グッシュ....これをきれいにするために、どこかありますか？あなたはcython_lapackに代わりcythonラッパーを使用することができます

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 

おかげで、

Answer 1

、http://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html

あなたはしかしcythonが必要になります。