LAPACK + C, le comportement bizarre

Question

Je suis en train de résoudre un système simple d'équations linéaires en utilisant LAPACK. J'utilise la méthode dbsvg qui est optimisée pour les matrices en bandes. J'ai obsereved un comportement étrange vraiment. 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 avec 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;

Toutes les idées?

Answer 1

Vous n'êtes pas exposant les entrées dans le stockage de bande correctement; il travaillait avant par un accident heureux. 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 des tableaux que vous aviez disposées.

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 seconde 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 -1,000000 -1,000000 2,000000 000 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. 0,000000 0,000000 0,000000 000000 0,000000 0,000000 0,000000 0,000000 0,000000 0 0,000000 0,000000 0,000000 0,000000

Apparemment, AT et AT2 sont les mêmes. Ce que je me attends.