2014-02-11 1 views
0

방금 ​​Cula를 다운로드했고 선형 방정식의 시스템을 해결하기 위해 구현 된 함수를 사용하고자합니다. 예제 디렉토리를 살펴 보았지만 코드 아래에서 보았지만 A * X = B의 X 솔루션을 얻고 싶을 때 매우 혼란 스럽습니다. 대답은, "B"있는 코드 아무것도이 줄을어디에 X에 대한 cula "culaSgesv"대답은 무엇입니까?

status = culaSgesv(N, NRHS, A, N, IPIV, X, N); 

발생하므로 X에서 그리고 이후 B 정체성 대각 행렬이다 (도움이되지 않았다 B로 변경 X를!)

는 것 무슨 일인지 말해 줄래? 이 질문에 "X"라는 대답을 어떻게 얻을 수 있습니까?

누군가 더 자세한 정보가 필요하면 알려주세요.

#ifdef CULA_PREMIUM 
void culaDoubleExample() 
{ 
#ifdef NDEBUG 
int N = 4096; 
#else 
int N = 780; 
#endif 
int NRHS = 1; 
int i; 

culaStatus status; 

culaDouble* A = NULL; 
culaDouble* B = NULL; 
culaDouble* X = NULL; 
culaInt* IPIV = NULL; 

culaDouble one = 1.0; 
culaDouble thresh = 1e-6; 
culaDouble diff; 

printf("-------------------\n"); 
printf("  DGESV\n"); 
printf("-------------------\n"); 

printf("Allocating Matrices\n"); 
A = (culaDouble*)malloc(N*N*sizeof(culaDouble)); 
B = (culaDouble*)malloc(N*sizeof(culaDouble)); 
X = (culaDouble*)malloc(N*sizeof(culaDouble)); 
IPIV = (culaInt*)malloc(N*sizeof(culaInt)); 
if(!A || !B || !IPIV) 
    exit(EXIT_FAILURE); 

printf("Initializing CULA\n"); 
status = culaInitialize(); 
checkStatus(status); 

// Set A to the identity matrix 
memset(A, 0, N*N*sizeof(culaDouble)); 
for(i = 0; i < N; ++i) 
    A[i*N+i] = one; 

// Set B to a random matrix (see note at top) 
for(i = 0; i < N; ++i) 
    B[i] = (culaDouble)rand(); 
memcpy(X, B, N*sizeof(culaDouble)); 

memset(IPIV, 0, N*sizeof(culaInt)); 

printf("Calling culaDgesv\n"); 
DWORD dw1 = GetTickCount(); 
status = culaDgesv(N, NRHS, A, N, IPIV, X, N); 

DWORD dw2 = GetTickCount(); 
cout<<"Time difference is "<<(dw2-dw1)<<" milliSeconds"<<endl; 
if(status == culaInsufficientComputeCapability) 
{ 
    printf("No Double precision support available, skipping example\n"); 
    free(A); 
    free(B); 
    free(IPIV); 
    culaShutdown(); 
    return; 
} 
checkStatus(status); 

printf("Verifying Result\n"); 
for(i = 0; i < N; ++i) 
{ 
    diff = X[i] - B[i]; 
    if(diff < 0.0) 
     diff = -diff; 
    if(diff > thresh) 
     printf("Result check failed: i=%d X[i]=%f B[i]=%f", i, X[i], B[i]); 
} 

printf("Shutting down CULA\n\n"); 
culaShutdown(); 

free(A); 
free(B); 
free(IPIV); 
} 

답변

1

당신은 Sgesv 언급하지만 당신이 보여준 샘플 코드는 Dgesv입니다. 그럼에도 불구하고 대답은 같습니다.

Netlib LAPACK reference 따르면, RHS 벡터의 B 행렬 6 파라미터로 함수에 전달된다

B

 B is DOUBLE PRECISION array, dimension (LDB,NRHS) 
     On entry, the N-by-NRHS matrix of right hand side matrix B. 
     On exit, if INFO = 0, the N-by-NRHS solution matrix X. 

[아웃에] 그리고 X 행렬 동일한 매개 변수 위치에 반환됩니다. 따라서 B에 함수에 전달할 때 NxNRHS 오른쪽 벡터가 포함되고 동일한 매개 변수는 X 결과를 반환합니다. 당신이 실제로 변수가 X라고하고, 결과는 개념은 그들이 아마도 혼란 B라는 이름의 변수에 대해 그것을 비교되지만이다 ((가) X 변수와 동일한에서) 반환 후 전달하는 표시 한 코드에서

똑같다. 예에서 A 행렬은 단위 행렬이기 때문에

Ax = b 정확한 솔루션이 일반적인 경우를 들어 6 번째 파라미터 위치 RHS 벡터중인 B 매트릭스를 통과해야 x=b

이다. 함수가 완료되면 결과 (X)가 동일한 매개 변수로 반환됩니다.

관련 문제