在 lapack 求解器 dgbsv 的行主要布局中编写带状矩阵

Question

我想求解这个线性方程组 Ax=b，其中：

A =  
    [-0.23  2.54  -3.66  0;  
     -6.98  2.46  -2.73  -2.73;  
      0     2.56   2.46   4.07;  
      0     0     -4.78   3.82]

b = [4.42 27.13 -6.14  10.5]

解决方案应该是

x = [-2  3  1  -4]

这是一个带状矩阵，下带等于 1，上带等于 2

使用DGBSV求解器求解如下

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#define N 4
int main() { 
    int i;
    MKL_INT ipiv[N];
    double a[16] = { 0,    -0.23,  2.54, -3.66,
                    -6.98,  2.46, -2.73, -2.13,
                     2.56,  2.46,  4.07,  0,
                    -4.78,  3.82,  0,     0};

    double b[4] =  {4.42, 27.13, -6.14, 10.5};

    LAPACKE_dgbsv( LAPACK_ROW_MAJOR, N, 1, 2, 1, a, N, ipiv, b, 1);

    for(i=0;i<N;i++)
        printf("%f\n", b[i]);
}

代码正在 dgbsv 解算器处中止。当我写矩阵 a 和 b 指针时，它给出了地址的值。

Answer 1

对于您问题中所述的输入，即

A =  
[-0.23  2.54  -3.66  0;  
 -6.98  2.46  -2.73  -2.73;  
  0     2.56   2.46   4.07;  
  0     0     -4.78   3.82]

b = [4.42 27.13 -6.14  10.5]

我得到的解决方案（将系统求解为密集系统）是：

[-3.77599, -1.28156, -1.85975, 0.421568]

不过，至于代码，有几点值得一提。什么函数LAPACKE_dgbsv essentially does is that it checks the validity of the input and in turn calls the function LAPACKE_dgbsv_work. If this function detects that the supplied layout is LAPACK_ROW_MAJOR, it merely transposes its input (matrix, right-hand sides) and passes all this to LAPACKE_dgbsv which expects column-major version of the packed matrix.

因此，如果您的代码指定 LAPACK_ROW_MAJOR，数组 a 应包含上面 link 中指定的压缩矩阵的行主版本。此外，也许更重要的是，LAPACKE_dgbsv 需要在数组 a 中添加额外的 space，以便它可以存储 LU 分解的结果。具体来说，必须有额外的 kl 行。

因此

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#define N 4
int main() {

    int i;
    MKL_INT ipiv[N];

    double b[4] =  {4.42, 27.13, -6.14, 10.5};

    double a[] = {
        0, 0, 0, 0,
        0, 0, -3.66, -2.73,
        0, 2.54, -2.73, 4.07,
        -0.23, 2.46, 2.46, 3.82,
        -6.98, 2.56, -4.78, 0};

    MKL_INT info = LAPACKE_dgbsv(LAPACK_ROW_MAJOR, N, 1, 2, 1, a, N, ipiv, b, 1);
    //printf("%d\n", info);

    for(i=0;i<N;i++) printf("%f\n", b[i]);
}

然后生成：

这与使用密集求解器获得的解一致。

Answer 2

可能的问题是gbsv functions expect matrix A to have space for 2KL + KU + 1 rows, whereas your code allocates only KL + KU + 1 rows. If this is wrong guess, you can look at the implementation of LAPACKE_dgbsv wrapper here。

在 lapack 求解器 dgbsv 的行主要布局中编写带状矩阵

writing a banded matrix in a row major layout for lapack solver dgbsv

linear-algebra

intel-mkl

lapacke