使用 dgtsv_ 或 sgtsv_ 求解类型 A*X = B 的方程

Question

我正在尝试求解以下类型的线性方程组：A*X = B in SWIFT。

我已经能够使用消耗 O(N^2) 内存的基于 LU 分解的算法来做到这一点。

由于我的数组通常很大（10000 个样本或更多），我正在查看 LAPACK，它具有一些特定于三对角矩阵的函数，这些函数仅消耗 O(N) 内存 space 并且效率更高。

http://www.netlib.org/lapack/explore-html-3.4.2/d4/d62/group__double_g_tsolve.html#

基本上，我希望使用上面的 dgtsv_ 或 sgtsv_ 函数求解方程。但是我找不到例子。

由于我是 SWIFT 的新手，我很难传递函数要求的 8 个输入参数。哪里有例子吗？

我粘贴在我的工作代码下方（使用 LU 分解）。

import Accelerate

func solve( A:[Double], _ B:[Double] ) -> [Double] {

var inMatrix:[Double] = A

var solution:[Double] = B

// Get the dimensions of the matrix. An NxN matrix has N^2
// elements, so sqrt( N^2 ) will return N, the dimension
var N:__CLPK_integer = __CLPK_integer( sqrt( Double( A.count ) ) )

// Number of columns on the RHS
var NRHS:__CLPK_integer = 1

// Leading dimension of A and B
var LDA:__CLPK_integer = N

var LDB:__CLPK_integer = N

// Initialize some arrays for the dgetrf_(), and dgetri_() functions
var pivots:[__CLPK_integer] = [__CLPK_integer](repeating: 0, count: Int(N))

var error: __CLPK_integer   = 0

// Perform LU factorization
dgetrf_(&N, &N, &inMatrix, &N, &pivots, &error)

// Calculate solution from LU factorization
_ = "T".withCString {
    dgetrs_( UnsafeMutablePointer(mutating: [=10=]), &N, &NRHS, &inMatrix, &LDA, &pivots, &solution, &LDB, &error )
}
return solution
  }


  //Call the function
  var A: [Double] = [
      1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]

  var b: [Double] = [0, -15, -15, -3, -3, 45, -12, -6, 0, 0]

  var cj = solve(A: A, b)

  print( cj ) // --> [0.0, -2.9185349611542728, -3.3258601553829079,   1.2219755826859044, -4.5620421753607099, 14.026193118756936, -6.5427302996670358, 0.14472807991120964, -0.036182019977802411, 0.0]
  //Call the function


  //TRY LAPACK (need examples to get above solution)
  let xx = dgtsv_(<#T##__n:    UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__nrhs: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__dl: UnsafeMutablePointer<__CLPK_doublereal>!##UnsafeMutablePointer<__CLPK_doublereal>!#>, <#T##__d__: UnsafeMutablePointer<__CLPK_doublereal>!##UnsafeMutablePointer<__CLPK_doublereal>!#>, <#T##__du: UnsafeMutablePointer<__CLPK_doublereal>!##UnsafeMutablePointer<__CLPK_doublereal>!#>, <#T##__b: UnsafeMutablePointer<__CLPK_doublereal>!##UnsafeMutablePointer<__CLPK_doublereal>!#>, <#T##__ldb: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__info: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>)

  let xx2 = sgtsv_(<#T##__n: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__nrhs: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__dl: UnsafeMutablePointer<__CLPK_real>!##UnsafeMutablePointer<__CLPK_real>!#>, <#T##__d__: UnsafeMutablePointer<__CLPK_real>!##UnsafeMutablePointer<__CLPK_real>!#>, <#T##__du: UnsafeMutablePointer<__CLPK_real>!##UnsafeMutablePointer<__CLPK_real>!#>, <#T##__b: UnsafeMutablePointer<__CLPK_real>!##UnsafeMutablePointer<__CLPK_real>!#>, <#T##__ldb: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>, <#T##__info: UnsafeMutablePointer<__CLPK_integer>!##UnsafeMutablePointer<__CLPK_integer>!#>)
  //TRY LAPACK (need examples to get above solution)

Answer 1

dgtsv_() 期望 tri-diagonal 的 lower/main/upper 对角线矩阵作为单独的参数。您可以使用 &.

传递变量数组的地址

所有整数参数都是 __CLPK_integer 又名 Int32 的地址变量。

right-hand 边向量 b 被解决方案 x 覆盖为 A x = b 等式。描述 A 的三个向量被覆盖同样，因此您可能想要复制原始数据。

示例：

import Swift
import Accelerate

var mainDiagA = [ 1.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 1.0 ]
var upperDiagA = [ 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ]
var lowerDiagA = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0 ]

var b = [0.0, -15.0, -15.0, -3.0, -3.0, 45.0, -12.0, -6.0, 0.0, 0.0 ]

var n = Int32(mainDiagA.count) // Order of matrix A
var nrhs = Int32(1) // Number of right-hand sides
var info = Int32(0) // Result code

dgtsv_(&n, &nrhs, &lowerDiagA, &mainDiagA, &upperDiagA, &b, &n, &info)
if info == 0 { // success
    print(b)
    // [0.0, -2.9185349611542732, -3.3258601553829075, 1.2219755826859044, -4.5620421753607099, 14.026193118756938, -6.5427302996670367, 0.14472807991120964, -0.036182019977802411, 0.0]

}

使用 dgtsv_ 或 sgtsv_ 求解类型 A*X = B 的方程

Solve Equations Of Type A*X = B Using dgtsv_ or sgtsv_

math

lapack

swift