拉帕克+ C、奇怪的行为

发布于 2024-08-19 02:16:46 字数 797 浏览 11 评论 0原文

我正在尝试使用 LAPACK 求解简单的线性方程组。我使用针对带状矩阵优化的 dbsvg 方法。我观察到一种非常奇怪的行为。当我这样填充 AT 矩阵时：

for(i=0; i<DIM;i++) AB[0][i] = -1;
for(i=0; i<DIM;i++) AB[1][i] = 2;
for(i=0; i<DIM;i++) AB[2][i] = -1;
for(i=0; i<3; i++)
    for(j=0;j<DIM;j++) {
        AT[i*DIM+j]=AB[i][j];
    }

并调用：

dgbsv_(&N, &KL, &KU, &NRHS, AT, &LDAB, myIpiv, x, &LDB, &INFO);

它工作得很好。然而，当我这样做时：

for(i=0; i<DIM;i++) AT[i] = -1;
for(i=0; i<DIM;i++) AT[DIM+i] = 2;
for(i=0; i<DIM;i++) AT[2*DIM+i] = -1;

它会产生一个填充 NaN 的向量。以下是声明：

double AB[3][DIM], AT[3*DIM];
double x[DIM];
int myIpiv[DIM];
int N=DIM, KL=1, KU=1, NRHS=1, LDAB=DIM, LDB=DIM, INFO;

有什么想法吗？

原文

I am trying to solve a simple linear equations system using LAPACK. I use dbsvg method which is optimised for banded matrices. I've obsereved a realy strange behaviour. When I fill the AT matrix this way:

for(i=0; i<DIM;i++) AB[0][i] = -1;
for(i=0; i<DIM;i++) AB[1][i] = 2;
for(i=0; i<DIM;i++) AB[2][i] = -1;
for(i=0; i<3; i++)
    for(j=0;j<DIM;j++) {
        AT[i*DIM+j]=AB[i][j];
    }

And call:

dgbsv_(&N, &KL, &KU, &NRHS, AT, &LDAB, myIpiv, x, &LDB, &INFO);

It works perfectly. However, when I do it this way:

for(i=0; i<DIM;i++) AT[i] = -1;
for(i=0; i<DIM;i++) AT[DIM+i] = 2;
for(i=0; i<DIM;i++) AT[2*DIM+i] = -1;

It results with a vector filled with NaN. Here are the declarations:

double AB[3][DIM], AT[3*DIM];
double x[DIM];
int myIpiv[DIM];
int N=DIM, KL=1, KU=1, NRHS=1, LDAB=DIM, LDB=DIM, INFO;

Any ideas?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

东北女汉子 2024-08-26 02:16:46

您没有正确布置带存储中的条目；之前它能工作是因为一次意外的意外。 LAPACK 文档说：

<块引用>

 输入时，矩阵 A 在带存储中，在行 KL+1 到
    2*KL+KU+1；无需设置数组的第 1 行至 KL 行。
    A 的第 j 列存储在
    数组AB如下：
    AB(KL+KU+1+ij,j) = A(i,j) 对于 max(1,j-KU)<=i<=min(N,j+KL)
    退出时，分解的详细信息：U 存储为
    具有 KL+KU 上对角线的上三角带矩阵
    第 1 行到 KL+KU+1，以及期间使用的乘数
    因式分解存储在行 KL+KU+2 到 2*KL+KU+1 中。
    请参阅下文了解更多详情。

因此，如果您想要一个对角线上为 2、上下为 -1 的三对角矩阵，则布局应为：

 *  *  *  *  *  *  *  ...  *  *  *  *
 * -1 -1 -1 -1 -1 -1  ... -1 -1 -1 -1
 2  2  2  2  2  2  2  ...  2  2  2  2
-1 -1 -1 -1 -1 -1 -1  ... -1 -1 -1  *

在本例中 LDAB 应为 4。请记住，LAPACK 使用列优先布局，因此实际数组在内存中应如下所示：

{ *, *, 2.0, -1.0, *, -1.0, 2.0, -1.0, *, -1.0, 2.0, -1.0, ... }

dgbsv 为两个相同的数组提供了不同的结果，因为它正在读取数组的末尾您布置的数组。

You're not laying out the entries in the band storage properly; it was working before by a happy accident. The LAPACK docs say:

    On entry, the matrix A in band storage, in rows KL+1 to
    2*KL+KU+1; rows 1 to KL of the array need not be set.
    The j-th column of A is stored in the j-th column of the
    array AB as follows:
    AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL)
    On exit, details of the factorization: U is stored as an
    upper triangular band matrix with KL+KU superdiagonals in
    rows 1 to KL+KU+1, and the multipliers used during the
    factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
    See below for further details.

So if you want a tridiagonal matrix with 2 on the diagonal and -1 above and below, the layout should be:

 *  *  *  *  *  *  *  ...  *  *  *  *
 * -1 -1 -1 -1 -1 -1  ... -1 -1 -1 -1
 2  2  2  2  2  2  2  ...  2  2  2  2
-1 -1 -1 -1 -1 -1 -1  ... -1 -1 -1  *

LDAB should be 4 in this case. Bear in mind that LAPACK uses a column-major layout, so the actual array should be look like this in memory:

{ *, *, 2.0, -1.0, *, -1.0, 2.0, -1.0, *, -1.0, 2.0, -1.0, ... }

dgbsv was giving different results for the two identical arrays because it was reading off the ends of the arrays that you had laid out.

回复收藏 0 原文

淡笑忘祈一世凡恋 2024-08-26 02:16:46

这是您使用的确切代码还是只是一个示例？我在这里运行了这段代码（只是从您的帖子中剪切并粘贴，在第二个循环中将 AT 更改为 AT2：

const int DIM=10;
double AB[DIM][DIM], AT[3*DIM], AT2[3*DIM];
int i,j;

for(i=0; i<DIM;i++) AB[0][i] = -1;
for(i=0; i<DIM;i++) AB[1][i] = 2;
for(i=0; i<DIM;i++) AB[2][i] = -1;
for(i=0; i<3; i++)
        for(j=0;j<DIM;j++) {
                AT[i*DIM+j]=AB[i][j];
        }
printf("AT:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT[i]);
printf("\n\n");
for(i=0; i<DIM;i++) AT2[i] = -1;
for(i=0; i<DIM;i++) AT2[DIM+i] = 2;
for(i=0; i<DIM;i++) AT2[2*DIM+i] = -1;
printf("AT2:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT2[i]);
printf("\n\n");
printf("Diff:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT[i]-AT2[i]);
printf("\n\n");

并得到此输出

AT：-1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000
00 -1.000000 -1.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.0
00000 2.000000 2.000000 2.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000
00 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000

AT2:-1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000
000 -1.000000 -1.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.
000000 2.000000 2.000000 2.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000
000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000

差异：0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.0
00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.
000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0
.000000 0.000000 0.000000 0.000000

显然 AT 和 AT2 是相同的。这是我所期望的。

Is this the exact code you used or just an example? I ran this code here (just cut and pasted from your posts, with a change of AT to AT2 in the second loop:

const int DIM=10;
double AB[DIM][DIM], AT[3*DIM], AT2[3*DIM];
int i,j;

for(i=0; i<DIM;i++) AB[0][i] = -1;
for(i=0; i<DIM;i++) AB[1][i] = 2;
for(i=0; i<DIM;i++) AB[2][i] = -1;
for(i=0; i<3; i++)
        for(j=0;j<DIM;j++) {
                AT[i*DIM+j]=AB[i][j];
        }
printf("AT:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT[i]);
printf("\n\n");
for(i=0; i<DIM;i++) AT2[i] = -1;
for(i=0; i<DIM;i++) AT2[DIM+i] = 2;
for(i=0; i<DIM;i++) AT2[2*DIM+i] = -1;
printf("AT2:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT2[i]);
printf("\n\n");
printf("Diff:");
for (i=0;i<3*DIM;++i) printf("%lf ",AT[i]-AT2[i]);
printf("\n\n");

and got this output

AT:-1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000
00 -1.000000 -1.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.0
00000 2.000000 2.000000 2.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000
00 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000

Diff:0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.0
00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.
000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0
.000000 0.000000 0.000000 0.000000

Apparently AT and AT2 are the same. Which I would expect.

回复收藏 0 原文

~没有更多了~