Lapack + c, странное поведение
-
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;
Это приводит к вектору, заполненному 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;
Любые идеи?
Решение
Вы не выкладываете записи в хранилище полосы должным образом; Это работало до счастливой аварии. А Лапак документов сказать:
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 в этом случае. Имейте в виду, что Лапак использует макет колонн-мажор, поэтому фактический массив должен выглядеть в памяти:
{ *, *, 2.0, -1.0, *, -1.0, 2.0, -1.0, *, -1.0, 2.0, -1.0, ... }
dgbsv
давал разные результаты для двух идентичных массивов, потому что он читал концы массивов, которые вы изложили.
Другие советы
Это точный код, который вы использовали, или просто пример? Я запустил этот код здесь (просто вырезал и вставил из ваших постов, с изменением на 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,0000000000000000000000002,00000000 2,000000 -100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000:00000000000000000000. -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
Видимо в и AT2 одинаковы. Что я ожидал бы.