2016-08-29 1 views
0

dstevr은 삼각형 대칭 행렬의 고유 해를 계산합니다. 시원한. 그것은 SCIPY의 래퍼 (wrapper)로 포팅 된 루틴 중 하나가 아니었다. 그래서 파이썬에서 직접 MKL을 호출하는 방법에 대한 지침을 따라 왔으며 첨부 된 내용이 정확한 대답을주는 것처럼 보입니다. 하지만 어 ... 이걸 정리할 여유가 있니?!Python에서 MKL 호출 : DSTEVR

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 

감사합니다,

답변

관련 문제