S03E13： 用 LAPACK 求解线性方程组

说明

我们已经知道可以用 simd 求解线性方程组，但是 simd 最多只支持 4x4，当需要进行更大范围的方程组求解时（比如 AutoLayout 的布局计算）就会需要求解更大的方程组，当然还是用 LAPACK 来求解了，毕竟苹果自家的Accelerate框架中自带了相关函数。

用法

用 LAPACK 解方程组，可以使用sgesv_() -> Int32函数，开头的s代表着单精度，即Float类型，如果需要使用双精度Double类型，也可以使用dgesv_() -> Int32函数。

和 SVD 类似，我们还需要使用苹果封装好的类Matrix。sgesv_()函数的另一个强大之处，在于可以同时求多组结果。用几何的方式描述就是：已知变换矩阵的情况下，可以同时求出多个点，在变换前的坐标。用代数方式描述：如下，当系数矩阵相同时，可以同时求出多组(x, y, z)对应的解。

extension Matrix {
    /// 方程组求解，a为系数矩阵，b 为结果（可多组），返回值实际存储在 b 中
    public static func sv(a:Matrix,b:Matrix) -> Matrix {
        //info = 0, 程序正常运行结束。
        var info = Int32(0)
        //矩阵的size，系数矩阵应该是方阵，一个5x5的矩阵，则n=5；
        var n = a.n
        //rhs的列数，LAPACK可以同时对一个系数矩阵，多个rhs进行求解，因为对系数矩阵只需要进行一次LU分解，多个rhs一起求解更方便；
        var nrhs = Int32(b.n)
        //leading dimension of a，lda>=max(1,n)；
        var lda = n
        //leading dimension of b, ldb>=max(1,n)
        var ldb = n
        //pivot的行交换记录，第i行被交换到第ipiv[i]行；
        var ipiv = Array(repeating: Int32(0), count: Int(n))
        //求解后的结果存储在 b 中
        sgesv_(&n, &nrhs, a.data.baseAddress, &lda, &ipiv, b.data.baseAddress, &ldb, &info)
        
        if info < 0 {
            fatalError("`sgesv_` info = \(info)")
        } else if info > 0 {
            fatalError("algorithmDidNotConverge")
        }
        return b
    }
}

代码

使用方法如下：

let nums:[Float] = [1,1,1,0,-1,2,2,-1,0]
let m = Matrix(source: nums, rowCount: 3, columnCount: 3)
let r = Matrix(source: [1,2,3,3,2,1], rowCount: 3, columnCount: 2)
let x = Matrix.sv(a: m, b: r)
print(x)
/*
 2.00 2.00
0.50 -0.50
-0.50 0.50
 */

项目代码

本系列文章代码已发布在 github：ComputationalGeometrySwift