float support for lapack functions. (#3)

* float support for lapack functions.

* adjusted lint rules.

Author: dastrobu

.swiftlint.yml

@@ -3,6 +3,7 @@ disabled_rules: # rule identifiers to exclude from running
   - identifier_name
   - force_try
   - line_length
+  - file_length
 function_body_length:
   - 100
 excluded:

AccelerateArray.podspec

@@ -16,7 +16,7 @@ Pod::Spec.new do |s|
   #
 
   s.name         = "AccelerateArray"
-  s.version      = "0.2.0"
+  s.version      = "0.3.0"
   s.summary      = "Swift Array Extensions for the Apple Accelerate Framework"
 
   # This description is used to generate tags and improve search results.

README.md

@@ -23,7 +23,7 @@ additional types, which can be easily built on top of this package.
 
 ### Swift Package Manager
     dependencies: [
-            .package(url: "https://github.com/dastrobu/AccelerateArray.git", from: "0.2.0"),
+            .package(url: "https://github.com/dastrobu/AccelerateArray.git", from: "0.3.0"),
         ],
         
 ### Cocoa Pods

Sources/AccelerateArray/lapack.swift

@@ -15,7 +15,207 @@ public enum LapackError: Error {
 ///
 /// Float array extension
 public extension Array where Element == Float {
+    /// SGETRF computes an LU factorization of a general M-by-N matrix A
+    /// using partial pivoting with row interchanges.
+    ///
+    /// The factorization has the form
+    ///    A = P * L * U
+    /// where P is a permutation matrix, L is lower triangular with unit
+    /// diagonal elements (lower trapezoidal if m > n), and U is upper
+    /// triangular (upper trapezoidal if m < n).
+    ///
+    /// This is the right-looking Level 3 BLAS version of the algorithm.
+    ///
+    /// This array must be in column major storage.
+    ///
+    /// http://www.netlib.org/lapack/explore-html/d8/ddc/group__real_g_ecomputational_ga8d99c11b94db3d5eac75cac46a0f2e17.html#ga8d99c11b94db3d5eac75cac46a0f2e17
+    ///
+    /// - Parameters:
+    ///     - m: number of rows
+    ///     - n: number of columns
+    ///
+    /// - Returns: The pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row IPIV(i).
+    mutating func getrf(m: Int, n: Int) throws -> [Int32] {
+        var ipiv = [Int32](repeating: 0, count: Swift.min(m, n))
+        try getrf(m: m, n: n, ipiv: &ipiv)
+        return ipiv
+    }
+
+    /// SGETRF computes an LU factorization of a general M-by-N matrix A
+    /// using partial pivoting with row interchanges.
+    ///
+    /// The factorization has the form
+    ///    A = P * L * U
+    /// where P is a permutation matrix, L is lower triangular with unit
+    /// diagonal elements (lower trapezoidal if m > n), and U is upper
+    /// triangular (upper trapezoidal if m < n).
+    ///
+    /// This is the right-looking Level 3 BLAS version of the algorithm.
+    ///
+    /// This array must be in column major storage.
+    ///
+    /// http://www.netlib.org/lapack/explore-html/d8/ddc/group__real_g_ecomputational_ga8d99c11b94db3d5eac75cac46a0f2e17.html#ga8d99c11b94db3d5eac75cac46a0f2e17
+    ///
+    /// - Parameters:
+    ///     - m: number of rows
+    ///     - n: number of columns
+    ///     - ipiv: the pivot indices; for 1 <= i <= min(m,n), row i of the matrix was interchanged with row ipiv[i].
+    mutating func getrf(m m_: Int, n n_: Int, ipiv: inout [Int32]) throws {
+        var m = Int32(m_)
+        assert(m >= 0, "\(m) >= 0")
+        var n = Int32(n_)
+        assert(n >= 0, "\(n) >= 0")
+        // leading dimension is the number of rows in column major order
+        var lda = Int32(m)
+        assert(lda >= Swift.max(1, m), "\(lda) >= max(1, \(m)")
+        assert(count == lda * n, "\(count) == (\(lda),\(n))")
+
+        assert(ipiv.count >= Swift.min(m, n), "\(ipiv.count) > min(\(m), \(n)")
+
+        var info: Int32 = 0
+        sgetrf_(&m, &n, &self, &lda, &ipiv, &info)
+        if info != 0 {
+            throw LapackError.getrf(info)
+        }
+    }
+
+    /// SGETRI computes the inverse of a matrix using the LU factorization
+    /// computed by DGETRF.
+    ///
+    /// This method inverts U and then computes inv(A) by solving the system
+    /// inv(A)*L = inv(U) for inv(A).
+    ///
+    ///
+    mutating func getri() throws {
+        var n = Int32(Double(count).squareRoot())
+        assert(count == n * n, "\(count) == \(n) * \(n)")
+        var ipiv = try getrf(m: Int(n), n: Int(n))
+
+        var lda = Int32(count / Int(n))
+        assert(lda >= Swift.max(1, n), "\(lda) >= max(1, \(n)")
+
+        var info: Int32 = 0
+
+        // do optimal workspace query
+        var lwork: Int32 = -1
+        var work = [__CLPK_real](repeating: 0.0, count: 1)
+        sgetri_(&n, &self, &lda, &ipiv, &work, &lwork, &info)
+        if info != 0 {
+            throw LapackError.getri(info)
+        }
+
+        // retrieve optimal workspace
+        lwork = Int32(work[0])
+        work = [__CLPK_real](repeating: 0.0, count: Int(lwork))
+
+        // do the inversion
+        sgetri_(&n, &self, &lda, &ipiv, &work, &lwork, &info)
+        if info != 0 {
+            throw LapackError.getri(info)
+        }
+    }
+
+    /// DGESV computes the solution to a real system of linear equations
+    ///    A * X = B,
+    /// where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
+    ///
+    /// The LU decomposition with partial pivoting and row interchanges is
+    /// used to factor A as
+    ///    A = P * L * U,
+    /// where P is a permutation matrix, L is unit lower triangular, and U is
+    /// upper triangular.  The factored form of A is then used to solve the
+    /// system of equations A * X = B.
+    ///
+    /// This array must be in column major storage.
+    ///
+    /// http://www.netlib.org/lapack/explore-html/d8/ddc/group__real_g_ecomputational_ga1af62182327d0be67b1717db399d7d83.html#ga1af62182327d0be67b1717db399d7d83
+    mutating func gesv(B: inout [Element]) throws {
+        var ipiv: [Int32] = [Int32].init(repeating: 0, count: n)
+        try gesv(ipiv: &ipiv, B: &B)
+    }
+
+    /// SGESV computes the solution to a real system of linear equations
+    ///    A * X = B,
+    /// where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
+    ///
+    /// The LU decomposition with partial pivoting and row interchanges is
+    /// used to factor A as
+    ///    A = P * L * U,
+    /// where P is a permutation matrix, L is unit lower triangular, and U is
+    /// upper triangular.  The factored form of A is then used to solve the
+    /// system of equations A * X = B.
+    ///
+    /// This array and B must be in column major storage.
+    ///
+    /// http://www.netlib.org/lapack/explore-html/d8/ddc/group__real_g_ecomputational_ga461f4ac32685a5ca30e293ee73d32920.html#ga461f4ac32685a5ca30e293ee73d32920
+    mutating func gesv(ipiv: inout [Int32], B: inout [Element]) throws {
+        var n = Int32(self.n)
+        assert(count == n * n, "\(count) == \(n) * \(n)")
+        assert(ipiv.count == n, "\(ipiv.count) == \(n)")
 
+        var nrhs = Int32(B.count / Int(n))
+        assert(nrhs >= 1, "\(nrhs) >= 1")
+        assert(B.count == nrhs * n, "\(B.count) == \(nrhs) * \(n)")
+
+        var lda = n
+        assert(lda >= Swift.max(1, n), "\(lda) >= max(1, \(n)")
+
+        var ldb = Int32(B.count / Int(nrhs))
+        assert(ldb * nrhs == B.count, "\(ldb) * \(nrhs) == \(B.count)")
+        assert(ldb >= Swift.max(1, n), "\(ldb) >= max(1, \(n))")
+
+        var info: Int32 = 0
+        sgesv_(&n, &nrhs, &self, &lda, &ipiv, &B, &ldb, &info)
+        if info != 0 {
+            throw LapackError.dgesv(info)
+        }
+    }
+
+    /// SGTSV  solves the equation
+    ///
+    ///    A*X = B,
+    ///
+    /// where A is an n by n tridiagonal matrix, by Gaussian elimination with
+    /// partial pivoting.
+    ///
+    /// Note that the equation  A**T*X = B  may be solved by interchanging the
+    /// order of the arguments DU and DL.
+    ///
+    /// This array represents the diagonal of A.
+    ///
+    /// http://www.netlib.org/lapack/explore-html/d1/d88/group__real_g_tsolve_gae1cbb7cd9c376c9cc72575d472eba346.html#gae1cbb7cd9c376c9cc72575d472eba346
+    ///
+    /// - Parameters:
+    ///     - nrhs: The number of right hand sides, i.e., the number of columns of the matrix B.  NRHS >= 0.
+    ///     - dl: On entry, DL must contain the (n-1) sub-diagonal elements of A.
+    ///           On exit, DL is overwritten by the (n-2) elements of the
+    ///           second super-diagonal of the upper triangular matrix U from
+    ///           the LU factorization of A, in DL(1), ..., DL(n-2).
+    ///     - du: On entry, DU must contain the (n-1) super-diagonal elements of A.
+    ///           On exit, DU is overwritten by the (n-1) elements of the first
+    ///           super-diagonal of U.
+    ///     - B:   On entry, the N by NRHS matrix of right hand side matrix B.
+    ///            On exit, if no error was thrown, the N by NRHS solution matrix X.
+    ///
+    mutating func gtsv(nrhs: Int, dl: inout [Element], du: inout [Element], B: inout [Element]) throws {
+        assert(count - 1 == dl.count, "\(count) - 1 == \(dl.count)")
+        assert(count - 1 == du.count, "\(count) - 1 == \(du.count)")
+        var n = Int32(count)
+
+        var nrhs = Int32(B.count / Int(n))
+        assert(nrhs >= 1, "\(nrhs) >= 1")
+        assert(B.count == Int(nrhs) * count, "\(B.count) == \(nrhs) * \(n)")
+
+        var ldb = Int32(B.count / Int(nrhs))
+        assert(ldb * nrhs == B.count, "\(ldb) * \(nrhs) == \(B.count)")
+        assert(ldb >= Swift.max(1, n), "\(ldb) >= max(1, \(n))")
+
+        var info: Int32 = 0
+        sgtsv_(&n, &nrhs, &dl, &self, &du, &B, &ldb, &info)
+        if info != 0 {
+            throw LapackError.dgesv(info)
+        }
+    }
 }
 
 /// Array extension employing the LAPACK framework.
@@ -206,7 +406,7 @@ public extension Array where Element == Double {
     ///     - B:   On entry, the N by NRHS matrix of right hand side matrix B.
     ///            On exit, if no error was thrown, the N by NRHS solution matrix X.
     ///
-    mutating func gtsv(nrhs: Int, dl: inout [Double], du: inout [Double], B: inout [Double]) throws {
+    mutating func gtsv(nrhs: Int, dl: inout [Element], du: inout [Element], B: inout [Element]) throws {
         assert(count - 1 == dl.count, "\(count) - 1 == \(dl.count)")
         assert(count - 1 == du.count, "\(count) - 1 == \(du.count)")
         var n = Int32(count)

Tests/AccelerateArrayTests/lapack.swift

@@ -2,6 +2,144 @@ import XCTest
 @testable import AccelerateArray
 
 class LapackTests: XCTestCase {
+    func testGetrfFloat() throws {
+        // A in row major
+        let A: [Float] = [
+            1.0, 2.0,
+            3.0, 4.0,
+            5.0, 6.0,
+            7.0, 8.0
+        ]
+        // convert A to col major
+        var At = A.mtrans(m: 2, n: 4)
+        let ipiv = try At.getrf(m: 4, n: 2)
+        // convert solution to row major
+        let X = At.mtrans(m: 4, n: 2)
+        // L in row major
+        let L: [Float] = [
+            1.0, 0.0,
+            X[2], 1.0,
+            X[4], X[5],
+            X[6], X[7],
+        ]
+        // U in row major
+        let U: [Float] = [
+            X[0], X[1],
+            0.0, X[3],
+        ]
+
+        // note, the indices in ipiv are one base (fortran)
+        // construct the permutation vector
+        // see: https://math.stackexchange.com/a/3112224/91477
+        var p = [0, 1, 2, 3]
+        for i in 0..<ipiv.count {
+            p.swapAt(i, Int(ipiv[i] - 1))
+        }
+
+        let n = 4
+        var P: [Float] = Array(repeating: 0, count: n * n)
+        for i in 0..<p.count {
+            // i iterates columns of P (in row major)
+            // p[i] indicates which element in the column must be set to one, to create the permutation matrix
+            P[i + p[i] * n] = 1.0
+        }
+
+        let PLU = P.mmul(B: L.mmul(B: U, m: 4, n: 2, p: 2), m: 4, n: 2, p: 4)
+        XCTAssertEqual(A, PLU, accuracy: 1e-6)
+    }
+
+    func testGetriFloat() throws {
+        // inversion is independent of row/col major storage
+        var A: [Float] = [
+            1.0, 2.0,
+            3.0, 4.0,
+        ]
+        try A.getri()
+
+        let Ainv: [Float] = [
+            -2.0, 1.0,
+            1.5, -0.5,
+        ]
+
+        XCTAssertEqual(A, Ainv)
+    }
+
+    func testGesvFloat() throws {
+        // A in row major
+        let A: [Float] = [
+            1.0, 2.0,
+            3.0, 4.0,
+        ]
+        // convert A to col major
+        var At = A.mtrans(m: 2, n: 2)
+        let b: [Float] = [
+            1.0, 1.0
+        ]
+        // B in row major
+        let B: [Float] = [
+            1.0, 2.0, 3.0,
+            1.0, 2.0, 3.0,
+        ]
+        // B in col major
+        var Bt = B.mtrans(m: 3, n: 2)
+        // Ainv in row major
+        let Ainv: [Float] = [
+            -2.0, 1.0,
+            1.5, -0.5,
+        ]
+        let x1 = Ainv.mmul(B: b, m: 2, n: 1, p: 2)
+        // X1 is in row major
+        let X1 = Ainv.mmul(B: B, m: 2, n: 3, p: 2)
+
+        // solution is stored col major in Bt
+        try At.gesv(B: &Bt)
+        let X2 = Bt.mtrans(m: 2, n: 3)
+
+        XCTAssertEqual(x1[0], X1[0], accuracy: 1e-15)
+        XCTAssertEqual(x1[1], X1[3], accuracy: 1e-15)
+        XCTAssertEqual(X1, X2, accuracy: 1e-6)
+    }
+
+    func testGtsvFloat() throws {
+        // A in row major
+        let A: [Float] = [
+            1.0, 1.0, 0.0, 0.0,
+            -1.0, 2.0, 2.0, 0.0,
+            0.0, -2.0, 3.0, 3.0,
+            0.0, 0.0, -3.0, 4.0,
+        ]
+        // convert A to col major
+        var At = A.mtrans(m: 4, n: 4)
+
+        // diagonals of A
+        var d: [Float] = [1.0, 2.0, 3.0, 4.0, ]
+        var du: [Float] = [1.0, 2.0, 3.0, ]
+        var dl: [Float] = [-1.0, -2.0, -3.0, ]
+
+        // B in row major
+        let B: [Float] = [
+            1.0, 2.0,
+            1.0, 2.0,
+            1.0, 2.0,
+            1.0, 2.0,
+        ]
+        // B in col major
+        var Bt = B.mtrans(m: 2, n: 4)
+        // make a copy of Bt
+        var Ct = Bt
+
+        // solve with general solver
+        // solution is stored col major in Ct
+        try At.gesv(B: &Ct)
+        let X1 = Ct.mtrans(m: 4, n: 2)
+
+        // solution is stored col major in Bt
+        try d.gtsv(nrhs: 2, dl: &dl, du: &du, B: &Bt)
+        let X2 = Bt.mtrans(m: 4, n: 2)
+
+        XCTAssertEqual(X1, X2, accuracy: 1e-15)
+    }
+
     func testGetrfDouble() throws {
         // A in row major
         let A: [Double] = [
@@ -142,6 +280,10 @@ class LapackTests: XCTestCase {
 
     static var allTests: [(String, (LapackTests) -> () throws -> Void)] {
         return [
+            ("testGetrfFloat", testGetrfFloat),
+            ("testGetriFloat", testGetriFloat),
+            ("testGesvFloat", testGesvFloat),
+            ("testGtsvFloat", testGtsvFloat),
             ("testGetrfDouble", testGetrfDouble),
             ("testGetriDouble", testGetriDouble),
             ("testGesvDouble", testGesvDouble),