LAPACK + C, comportement bizarre

Question

J'essaie de résoudre un système d'équations linéaires simple en utilisant LAPACK. J'utilise la méthode dbsvg qui est optimisée pour les matrices en bandes. J'ai observé un comportement vraiment étrange. Quand je remplis la matrice AT ainsi:

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]; 
    } 

Et appelez:

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

Il fonctionne parfaitement. Cependant, quand je le fais de cette façon:

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; 

Il en résulte un vecteur rempli de NaN. Voici les déclarations:

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; 

Des idées?

Answer 1

Vous ne disposez pas correctement les entrées dans le stockage de bande; ça fonctionnait avant par un heureux accident. Le LAPACK docs disent:

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. 

Donc, si vous voulez une matrice tridiagonal avec 2 sur la diagonale et -1 ci-dessus et ci-dessous, la mise en page doit être:

* * * * * * * ... * * * * 
* -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 devrait être 4 dans ce cas.Gardez à l'esprit que LAPACK utilise une mise en page de colonne principale, de sorte que le tableau réel devrait être ressembler à ceci en mémoire:

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

dgbsv donnait des résultats différents pour les deux tableaux identiques parce qu'il était en train de lire les extrémités du tableaux que vous aviez disposés.

Answer 2

Est-ce le code exact que vous avez utilisé ou juste un exemple? J'ai couru ce code ici (juste coupé et collé de vos messages, avec un changement d'AT à AT2 dans la deuxième boucle:

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

et a obtenu cette sortie

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

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

Apparemment A T et AT2 sont les mêmes. Ce que je m'attendrais.