cula「culaSgesv」在哪裏答案？

我剛剛下載了Cula，我想用它實現的函數來解決線性方程組的系統問題。我查看了實例目錄，並且看到了下面的代碼，但是當他們想要獲得A * X = B的X解決方案時，的B X和因爲A是身份對角矩陣，因此，答案是，「B」，並在這行代碼沒有任何反應cula「culaSgesv」在哪裏答案？

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

（改變X到B並沒有幫助！）

你會請告訴我發生了什麼事？請告訴我如何從這裏得到答案「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); 
}

來源

2014-02-11 Alexander1991

你提到Sgesv但你已經顯示的示例代碼是Dgesv。儘管如此，答案是一樣的。

按照Netlib LAPACK reference，RHS矢量的B矩陣傳遞給函數作爲第六參數：

[IN，OUT]乙

 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並將結果返回後（在相同的變量X）他們正在對一個名爲B這也許是混淆變量與它，但概念是代碼一樣。

由於A矩陣的例子是單位矩陣，Ax = b正確的解決辦法是x=b

對於一般情況下，你應該通過你的B矩陣RHS向量的第6個參數的位置。完成該功能後，結果（X）將以相同的參數返回。

來源

2014-02-12 06:22:56

cula「culaSgesv」在哪裏答案？

回答

相關問題