2015-01-11 86 views
1

我在加速框架中使用LAPACK中的dgeev算法來計算矩陣的特徵值和特徵向量。這是我的代碼:Swift中加速框架的問題

var matrix:[__CLPK_doublereal] = [1,2,3,4,5,6,7,8,9] 

    var N = __CLPK_integer(sqrt(Double(matrix.count))) 
    var pivots = [__CLPK_integer](count: Int(N), repeatedValue: 0) 
    var workspace = [Double](count: Int(N), repeatedValue: 0.0) 
    var error : __CLPK_integer = 0 
    var lwork = __CLPK_integer(-1) 
    // Real parts of eigenvalues 
    var wr = [Double](count: Int(N), repeatedValue: 0) 
    // Imaginary parts of eigenvalues 
    var wi = [Double](count: Int(N), repeatedValue: 0) 
    // Left eigenvectors 
    var vl = [__CLPK_doublereal](count: Int(N*N), repeatedValue: 0) 
    // Right eigenvectors 
    var vr = [__CLPK_doublereal](count: Int(N*N), repeatedValue: 0) 

    dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &matrix, &N, &wr, &wi, &vl, &N, &vr, &N, &workspace, &lwork, &error) 

    println("\(wr), \(vl), \(vr)") 

這隻打印包含零的數組,這意味着它們不被函數修改。我究竟做錯了什麼?

更新1

我現在有這樣的:

var matrix:[__CLPK_doublereal] = [1,2,3,4,5,6,7,8,9] 

    var N = __CLPK_integer(sqrt(Double(matrix.count))) 
    var pivots = [__CLPK_integer](count: Int(N), repeatedValue: 0) 
    var workspaceQuery = [Double](count: 1, repeatedValue: 0.0) 
    var error : __CLPK_integer = 0 
    var lwork = __CLPK_integer(-1) 
    // Real parts of eigenvalues 
    var wr = [Double](count: Int(N), repeatedValue: 0) 
    // Imaginary parts of eigenvalues 
    var wi = [Double](count: Int(N), repeatedValue: 0) 
    // Left eigenvectors 
    var vl = [__CLPK_doublereal](count: Int(N*N), repeatedValue: 0) 
    // Right eigenvectors 
    var vr = [__CLPK_doublereal](count: Int(N*N), repeatedValue: 0) 

    dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &matrix, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) 
    var workspace = [Double](count: Int(workspaceQuery[0]), repeatedValue: 0.0) 

    dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &matrix, &N, &wr, &wi, &vl, &N, &vr, &N, &workspace, &lwork, &error) 

    println("\(wr), \(vl), \(vr)") 

不過它打印零。

回答

3

該問題與lwork變量相關。這應該是您提供的工作空間的大小,與-1這意味着你正在執行「工作空間查詢」:

LWORK (input) INTEGER 
     The dimension of the array WORK. LWORK >= max(1,3*N), and 
     if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good 
     performance, LWORK must generally be larger. 

     If LWORK = -1, then a workspace query is assumed; the routine 
     only calculates the optimal size of the WORK array, returns 
     this value as the first entry of the WORK array, and no error 
     message related to LWORK is issued by XERBLA. 

所以,你可能想是這樣的:

var workspaceQuery: Double = 0.0 
dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &matrix, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) 

// prints "102.0" 
println("\(workspaceQuery)") 

// size workspace per the results of the query: 
var workspace = [Double](count: Int(workspaceQuery), repeatedValue: 0.0) 
lwork = __CLPK_integer(workspaceQuery) 

dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &matrix, &N, &wr, &wi, &vl, &N, &vr, &N, &workspace, &lwork, &error) 

// this now prints non-zero values 
println("\(wr), \(vl), \(vr)") 
+0

你值得你名稱:)我剛剛得出同樣的結論... –

+0

那麼,應該在第二次電話會議中「lwork」的價值是什麼? –

+0

@YoussefSami:我認爲它應該是工作區的實際大小:'lwork = __CLPK_integer(workspace.count)'。我發現這個C代碼的例子,它也演示了這個用法:https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev_ex.c.htm –