为 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("%fn", 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本质上的作用是检查输入的有效性，然后调用函数LAPACKE_dgbsv_work。如果此函数检测到提供的布局LAPACK_ROW_MAJOR，它只是转置其输入(矩阵，右侧)并将所有这些传递给LAPACKE_dgbsv，后者需要打包矩阵的列主版本。

因此，如果您的代码指定LAPACK_ROW_MAJOR，则数组a应包含上面链接中指定的打包矩阵的行主版本。此外，也许更重要的是，LAPACKE_dgbsv在数组a中需要额外的空间，以便它可以存储 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("%dn", info);
for(i=0;i<N;i++) printf("%fn", b[i]);
}

然后产生：

-3.775989
-1.281561
-1.859751
0.421568

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

Answer 2

问题可能是 gbsv 函数期望矩阵 A 有 2KL + KU + 1 行的空间，而您的代码只分配 KL + KU + 1 行。如果这是错误的猜测，你可以在这里查看LAPACKE_dgbsv包装器的实现。