阅读 114

S03E13: 用 LAPACK 求解线性方程组

说明

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

用法

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

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

[102111120][xyz]=[123]\left[ \begin{matrix} 1 & 0 & 2 \\ 1 & -1 & -1 \\ 1 & 2 & 0 \end{matrix} \right] * \left[ \begin{matrix} x \\ y \\ z \\ \end{matrix} \right] = \left[ \begin{matrix} 1 \\ 2 \\ 3 \\ \end{matrix} \right]
[102111120][xyz]=[321]\left[ \begin{matrix} 1 & 0 & 2 \\ 1 & -1 & -1 \\ 1 & 2 & 0 \end{matrix} \right] * \left[ \begin{matrix} x \\ y \\ z \\ \end{matrix} \right] = \left[ \begin{matrix} 3 \\ 2 \\ 1 \\ \end{matrix} \right]
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

文章分类
iOS
文章标签