لاباك + ج ، سلوك غريب

Question

أحاول حل نظام المعادلات الخطية البسيطة باستخدام 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;

أيه أفكار؟

Answer 1

أنت لا تضع الإدخالات في تخزين النطاق بشكل صحيح ؛ كان يعمل من قبل بحادث سعيد. ال مستندات لاباك قل:

    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 على القطري و -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 كان يعطي نتائج مختلفة للصفعتين المتطابقتين لأنه كان يقرأ من نهايات المصفوفات التي وضعتها.

Answer 2

هل هذا هو الرمز الدقيق الذي استخدمته أم مجرد مثال؟ قمت بتشغيل هذا الرمز هنا (فقط قطع ولصق من مشاركاتك ، مع تغيير في 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");

وحصلت على هذا المخرج

في: -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000 00 -1.000000 -1.000000 2.000000000000000000000000000000000000 -000000000000000000000000000000 2.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.00000000000000000000 -000000000000000000000000000000000000000000000000000000000000 2.000000. -1.000000 -1.000000 -1.000000

DIF: 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.0000000000000000000000000000000000000000000000000000000000 0.000000 0 0.0000

على ما يبدو في و AT2 هي نفسها. الذي أتوقعه.