2015-01-11 11 views
1

Я использую алгоритм dgeev из LAPACK в структуре Accelerate для вычисления собственных значений и собственных векторов матрицы. Вот мой код:Проблема с ускорением рамки в 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)") 

Это печатает только массивы, содержащие нули, что означает, что они не изменяются функцией. Что я делаю не так?

UPDATE 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 –