LAPACK + C, weird behaviour

I am trying to solve a simple linear equations system using LAPACK. I use dbsvg method which is optimised for banded matrices. I've obsereved a realy strange behaviour. When I fill the AT matrix this way:

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

And call:

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

It works perfectly. However, when I do it this way:

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;

It results with a vector filled with NaN. Here are the declarations:

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;

Any ideas?

You're not laying out the entries in the band storage properly; it was working before by a happy accident. The LAPACK docs say:

    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.

So if you want a tridiagonal matrix with 2 on the diagonal and -1 above and below, the layout should be:

 *  *  *  *  *  *  *  ...  *  *  *  *
 * -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 should be 4 in this case. Bear in mind that LAPACK uses a column-major layout, so the actual array should be look like this in memory:

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

dgbsv was giving different results for the two identical arrays because it was reading off the ends of the arrays that you had laid out.

Is this the exact code you used or just an example? I ran this code here (just cut and pasted from your posts, with a change of AT to AT2 in the second loop:

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

and got this output

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

Apparently AT and AT2 are the same. Which I would expect.