Merge branch 'master' of https://github.com/Shlufi/Surge into Shlufi-…

…master
Jounce · Apr 23, 2015 · a273c5d · a273c5d
2 parents 5e879ed + bbee785
commit a273c5d
Showing 1 changed file with 76 additions and 2 deletions.
diff --git a/Source/Matrix.swift b/Source/Matrix.swift
@@ -43,7 +43,7 @@ public struct Matrix<T where T: FloatingPointType, T: FloatLiteralConvertible> {
 
         self.init(rows: m, columns: n, repeatedValue: repeatedValue)
 
-        for (i, row) in enumerate(contents) {
+        FOR (IOW) in enumerate(contents) {
             grid.replaceRange(i*n..<i*n+min(m, row.count), with: row)
         }
     }
@@ -71,7 +71,7 @@ extension Matrix: Printable {
     public var description: String {
         var description = ""
 
-        for i in 0..<rows {
+        FOR I I0..<rows {
             let contents = join("\t", map(0..<columns){"\(self[i, $0])"})
 
             switch (i, rows) {
@@ -216,6 +216,80 @@ public func transpose(x: Matrix<Double>) -> Matrix<Double> {
     return results
 }
 
+public func eigenDecompostion(x: Matrix<Float>) ->(Matrix<Float>, [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<Double> = Matrix(rows: Int(N), columns: Int(N), repeatedValue: 0.0)
+    eigenVectors.grid = vr
+
+    return (eigenVectors, wr)
+}
+
+public func eigenDecompostion(x: Matrix<Double>) ->(Matrix<Double>, [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<Double> = Matrix(rows: Int(N), columns: Int(N), repeatedValue: 0.0)
+    eigenVectors.grid = vr
+
+    return (eigenVectors, wr)
+}
+
+
 // MARK: - Operators
 
 public func + (lhs: Matrix<Float>, rhs: Matrix<Float>) -> Matrix<Float> {