2015-01-11 52 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 –