Add dgetri_

Author: ptoffy

Sources/AccelerateLinux/MatrixOps/LAPACK.swift

@@ -23,7 +23,7 @@ public func dgesv_(
     _ __info: UnsafeMutablePointer<__CLPK_integer>!
 )
 
-/// DGESVD computes the singular value decomposition (SVD) of a real
+/// Computes the singular value decomposition (SVD) of a real
 /// M-by-N matrix A, optionally computing the left and/or right singular
 /// vectors. The SVD is written
 ///
@@ -55,7 +55,7 @@ public func dgesvd_(
     _ __info: UnsafeMutablePointer<__CLPK_integer>!
 )
 
-/// DGETRF computes an LU factorization of a general M-by-N matrix A
+/// Computes an LU factorization of a general M-by-N matrix A
 /// using partial pivoting with row interchanges.
 ///
 /// The factorization has the form
@@ -74,4 +74,20 @@ public func dgetrf_(
     _ __ipiv: UnsafeMutablePointer<__CLPK_integer>!,
     _ __info: UnsafeMutablePointer<__CLPK_integer>!
 ) -> Int32
+
+/// Computes the inverse of a matrix using the LU factorization
+/// computed by ``dgterf_``.
+///
+/// This method inverts U and then computes `inv(A)` by solving the system
+/// `inv(A)*L = inv(U) for inv(A)`.
+@_silgen_name("dgetri_")
+public func dgetri_(
+    _ __n: UnsafeMutablePointer<__CLPK_integer>!,
+    _ __a: UnsafeMutablePointer<__CLPK_doublereal>!,
+    _ __lda: UnsafeMutablePointer<__CLPK_integer>!,
+    _ __ipiv: UnsafeMutablePointer<__CLPK_integer>!,
+    _ __work: UnsafeMutablePointer<__CLPK_doublereal>!,
+    _ __lwork: UnsafeMutablePointer<__CLPK_integer>!,
+    _ __info: UnsafeMutablePointer<__CLPK_integer>!
+) -> Int32
 #endif  // canImport(Accelerate)

Tests/AccelerateLinuxTests/MatrixTests/LAPACKTests.swift

@@ -148,4 +148,45 @@ struct LAPACKTests {
             ]
         )
     }
+
+    // https://numericalalgorithmsgroup.github.io/LAPACK_Examples/examples/doc/dgetri_example.html
+    @Test("dgetri_")
+    func test_dgetri_() {
+        let M = 4
+        let N = 4
+        let LDA = M
+
+        var m = __CLPK_integer(M)
+        var n = __CLPK_integer(N)
+        var lda = __CLPK_integer(LDA)
+        var info = __CLPK_integer(0)
+
+        var ipiv = [__CLPK_integer](repeating: 0, count: Int(min(m, n)))
+
+        var a: [__CLPK_doublereal] = [
+            1.80, 5.25, 1.58, -1.11,
+            2.88, -2.95, -2.69, -0.66,
+            2.05, -0.95, -2.90, -0.59,
+            -0.89, -3.80, -1.04, 0.80,
+        ]
+
+        // A is nonsymmetric and must be factorized first
+        _ = dgetrf_(&m, &n, &a, &lda, &ipiv, &info)
+
+        var work = [__CLPK_doublereal](repeating: 0.0, count: Int(m))
+        var lwork = __CLPK_integer(64 * n)
+
+        _ = dgetri_(&n, &a, &lda, &ipiv, &work, &lwork, &info)
+
+        #expect(
+            a.map {
+                ($0 * pow(10, 2)).rounded() / pow(10, 2)
+            } == [
+                1.77, -0.12, 0.18, 2.49,
+                0.58, -0.45, 0.45, 0.76,
+                0.08, 0.41, -0.67, -0.04,
+                4.82, -1.71, 1.48, 7.61,
+            ]
+        )
+    }
 }