diff --git a/Source/Matrix.swift b/Source/Matrix.swift index 1ac5434c..eff22d2e 100644 --- a/Source/Matrix.swift +++ b/Source/Matrix.swift @@ -43,7 +43,7 @@ public struct Matrix { self.init(rows: m, columns: n, repeatedValue: repeatedValue) - for (i, row) in enumerate(contents) { + FOR (IOW) in enumerate(contents) { grid.replaceRange(i*n..) -> Matrix { return results } +public func eigenDecompostion(x: Matrix) ->(Matrix, [Float]) +{ + precondition(x.rows == x.columns, "Matrix must be square") + + var mat:[__CLPK_floatreal] = x.grid + + var N = __CLPK_integer(sqrt(Float(x.grid.count))) + var pivots = [__CLPK_integer](count: Int(N), repeatedValue: 0) + var workspaceQuery: Float = 0.0 + var error : __CLPK_integer = 0 + var lwork = __CLPK_integer(-1) + // Real parts of eigenvalues + var wr = [Float](count: Int(N), repeatedValue: 0) + // Imaginary parts of eigenvalues + var wi = [Float](count: Int(N), repeatedValue: 0) + + // Left eigenvectors + var vl = [Float](count: Int(N*N), repeatedValue: 0) + // Right eigenvectors + var vr = [Float](count: Int(N*N), repeatedValue: 0) + + // The first iteration is used to get the size of the workspace from the workspace query + dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &mat, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) + + var workspace = [Double](count: Int(workspaceQuery), repeatedValue: 0.0) + lwork = __CLPK_integer(workspaceQuery) + + // Actual iteration after obtaining workspace size + dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer,("V" as NSString).UTF8String), &N, &mat, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) + + var eigenVectors: Matrix = Matrix(rows: Int(N), columns: Int(N), repeatedValue: 0.0) + eigenVectors.grid = vr + + return (eigenVectors, wr) +} + +public func eigenDecompostion(x: Matrix) ->(Matrix, [Double]) +{ + precondition(x.rows == x.columns, "Matrix must be square") + + var mat:[__CLPK_doublereal] = x.grid + + var N = __CLPK_integer(sqrt(Double(x.grid.count))) + var pivots = [__CLPK_integer](count: Int(N), repeatedValue: 0) + var workspaceQuery: Double = 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 = [Double](count: Int(N*N), repeatedValue: 0) + // Right eigenvectors + var vr = [Double](count: Int(N*N), repeatedValue: 0) + + // The first iteration is used to get the size of the workspace from the workspace query + dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer(("V" as NSString).UTF8String), &N, &mat, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) + + var workspace = [Double](count: Int(workspaceQuery), repeatedValue: 0.0) + lwork = __CLPK_integer(workspaceQuery) + + // Actual run after obtaining workspace size + dgeev_(UnsafeMutablePointer(("V" as NSString).UTF8String), UnsafeMutablePointer,("V" as NSString).UTF8String), &N, &mat, &N, &wr, &wi, &vl, &N, &vr, &N, &workspaceQuery, &lwork, &error) + + var eigenVectors: Matrix = Matrix(rows: Int(N), columns: Int(N), repeatedValue: 0.0) + eigenVectors.grid = vr + + return (eigenVectors, wr) +} + + // MARK: - Operators public func + (lhs: Matrix, rhs: Matrix) -> Matrix {