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

그리고 전화 : 나는 AT 행렬이 방법을 채울 때

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. 

그래서 당신이 대각선 -1 위 아래 2와 삼중 대각 행렬을 원하는 경우, 레이아웃은 다음과 같아야합니다

* * * * * * * ... * * * * 
* -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하는 AT의 변화와 함께, 귀하의 게시물에서 잘라 붙여 실행 :

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

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 0.000000 0.000000

분명히 A T와 AT2는 동일합니다. 나는 어느 것을 기대할 것인가.