LAPACK dpbsv求解正定矩阵返回info=3的问题排查

我尝试使用LAPACK的带状对称矩阵求解器dpbsv，测试如下正定矩阵（经Mathematica验证，行列式为3684）：

4,  2,  0,   0, 0
2,  4,  3,   0, 0
0,  3, 11,  7,  0
0,  0,  7, 11,  5
0,  0,  0,  5, 13

用Swift构建数组：

var a: [Double] = [ 0, 2,  3,  7,  5,
                    4, 4, 11, 11, 13]
var b: [Double] = [1, 2, 3, 4, 5]

调用dpbsv_的代码如下：

var uplo = Int8("U".utf8.first!) // 设置为'U'

var n = __CLPK_integer(5)
var kd = __CLPK_integer(1)
var ldab = kd + 1
var nrhs = __CLPK_integer(1)
var ldb = __CLPK_integer(5)
var info: __CLPK_integer = 0

dpbsv_(&uplo,
       &n,
       &kd,
       &nrhs,
       &a,
       &ldab,
       &b,
       &ldb,
       &info)

if info != 0 {
    // info返回3，表示矩阵非正定
    NSLog("error \(info)")
}

但info返回3，提示矩阵非正定。更换其他Mathematica验证过的正定矩阵，结果一致。请问参数解读是否有误？问题出在哪里？

核心问题是带状矩阵的存储格式不符合LAPACK要求：

1. LAPACK的矩阵存储采用列优先顺序，而你用了行优先的逻辑构建数组a。
2. 当uplo='U'时，dpbsv要求的带状矩阵存储规则：
   • 数组a的行数为kd+1（这里kd=1，所以2行）
   • 第1行对应主对角线上方第kd条次对角线（这里是上方第1条），空缺位置补0
   • 第2行对应主对角线
   • 元素按列优先顺序排列（先存每一列的所有行，再依次存下一列）

你的原数组是按行优先写的（先写第一行的[0,2,3,7,5]，再写第二行的[4,4,11,11,13]），这完全不符合LAPACK的列优先存储要求，导致传入的矩阵并非你预期的正定矩阵，因此info=3报错。

修正后的数组 #

按照列优先规则，正确的a数组应该是：

var a: [Double] = [0, 4, 2, 4, 3, 11, 7, 11, 5, 13]

对应列优先的存储逻辑：

• 第1列：上方第1条对角线补0，主对角线元素4 → 0,4
• 第2列：上方第1条对角线元素2，主对角线元素4 → 2,4
• 第3列：上方第1条对角线元素3，主对角线元素11 → 3,11
• 第4列：上方第1条对角线元素7，主对角线元素11 → 7,11
• 第5列：上方第1条对角线元素5，主对角线元素13 → 5,13

验证 #

使用修正后的数组调用dpbsv_，info会返回0，求解过程正常。

内容的提问来源于stack exchange，提问作者dcsalmon