X 的 cula "culaSgesv" 答案在哪里？答案

【问题标题】：where is cula "culaSgesv" answer for X?X 的 cula "culaSgesv" 答案在哪里？
【发布时间】：2014-03-09 11:04:30
【问题描述】：

我刚刚下载了 Cula，我想使用它实现的函数来求解线性方程组在 X 中复制 B，因为 A 是单位对角矩阵，所以答案是“B”，在这行代码中没有任何反应

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

【问题讨论】：

标签： lapack equation-solving cula

【解决方案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 的变量，并在返回结果后（在同一个变量 X 中）将其与名为 B 的变量进行比较，这可能令人困惑，但概念是一样的。

由于例子中的A矩阵是单位矩阵，所以Ax = b的正确解是x=b

对于一般情况，您应该在第 6 个参数位置传递 RHS 向量的 B 矩阵。函数完成后，结果（X）将在同一个参数中返回。

【讨论】：