I am using the dgeev
algorithm from LAPACK in the Accelerate framework to calculate eigenvalues and eigenvectors of a matrix. Here's my code:
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)")
This prints only arrays containing zeroes which means they aren't being modified by the function. What am I doing wrong?
UPDATE 1
I now have this:
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)")
Still it prints zeroes.
The problem’s with your lwork
variable. This is supposed to be the size of the workspace you supply, with -1
meaning you’re performing a “workspace query”:
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.
So you probably want something like this:
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)")