LAPACK + C, странное поведение

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; 

Это результат с вектором, заполненным 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; 

Любые идеи?

Answer 1

Вы не храните записи в памяти группы должным образом; он работал раньше счастливой случайностью. LAPACK docs говорят:

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"); 

и получил этот выход

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

По-видимому, A T и AT2 одинаковы. Который я ожидал бы.