Lapack + C, comportamiento extraño

Question

Estoy tratando de resolver un sistema de ecuaciones lineales simples usando Lapack. Utilizo el método DBSVG que está optimizado para matrices con bandas. He observado un comportamiento realmente extraño. Cuando lleno el AT Matrix de esta manera:

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

Y llama:

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

Funciona perfectamente. Sin embargo, cuando lo hago de esta manera:

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;

Resulta con un vector lleno de Nan. Aquí están las declaraciones:

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;

¿Algunas ideas?

Answer 1

No estás diseñando las entradas en el almacenamiento de la banda correctamente; Estaba funcionando antes por un feliz accidente. los Docios de Lapack decir:

    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.

Entonces, si desea una matriz tridiagonal con 2 en la diagonal y -1 arriba y abajo, el diseño debe ser:

 *  *  *  *  *  *  *  ...  *  *  *  *
 * -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 debe ser 4 en este caso. Tenga en cuenta que Lapack usa un diseño de columna mayor, por lo que la matriz real debe verse así en la memoria:

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

dgbsv estaba dando diferentes resultados para las dos matrices idénticas porque estaba leyendo los extremos de las matrices que había presentado.

Answer 2

¿Es este el código exacto que usó o simplemente un ejemplo? Ejecuté este código aquí (solo corté y pegé de tus publicaciones, con un cambio de AT a AT2 en el segundo bucle:

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

y obtuve esta salida

En: -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000000000S -1.0000 0000 -1.000000 -0000000000ES 2.000000 2.000000 2.000000 2.0000000000ES --1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1. 1.000000 -1.000000 -1.000000

At2: -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.0000000000 -1.000 000 -1.000000 -1.000000 2.00000000 2.000000 2.000000 2.0000000000 --1. -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

Aparentemente en y AT2 son los mismos. Lo que esperaría.