1

我想求解这个线性系统 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 指针时，它给出了地址的值。

4

2 回答 2

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("%d\n", info);

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

然后产生：

-3.775989
-1.281561
-1.859751
0.421568

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

于 2017-03-13T11:10:22.837 回答
0

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

于 2017-03-20T12:10:11.033 回答