Lapack + C、奇妙な動作

StackOverflow https://stackoverflow.com/questions/2121199

  •  22-09-2019
  •  | 
  •  

質問

Lapackを使用して、単純な線形方程式システムを解決しようとしています。バンドマトリックスに最適化されたDBSVGメソッドを使用します。私は本当に奇妙な行動を観察しました。この方法でマトリックスを埋めるとき：

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;

ナンで満たされたベクトルで生じます。宣言は次のとおりです。

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;

何か案は？

役に立ちましたか？

解決

バンドストレージのエントリを適切にレイアウトしているわけではありません。それは以前に幸せな事故で働いていました。 ラパックドキュメント いう：

    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.

したがって、対角線に2が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 レイアウトしたアレイの端から読み取られていたため、2つの同一の配列に異なる結果をもたらしていました。

他のヒント

これはあなたが使用した正確なコードですか、それとも単なる例ですか？ここでこのコードを実行しました（あなたの投稿から切り取って貼り付けただけで、2番目のループでAT2に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.0000000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000 000 -1.000000 -1.000000 2.0000.000000.000000000000 -10000 -10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000- -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

どうやらATとAT2は同じです。私が期待するでしょう。

ライセンス： CC-BY-SA と 帰属
所属していません StackOverflow
scroll top