小贝子编程

CUDA cuSolver gesvdj with large matrix



我正在运行在"G.2.带有奇异矢量的 SVD(通过雅可比方法("，在 NVIDIA P6000 上。稍作修改是为了在堆中为 A、U 和 V 向量动态分配内存，并使用依赖于索引的值填充指定大小的 A 向量到 A 中。我还转换了从双打到浮动的所有内容。最后一个修改是对 gesvdj 调用本身的循环和收敛检查，最多进行一些 # 次迭代(在我的例子中为 10 次(。

通过这些细微的修改，我能够克服在大于 ~1000x1000 的对称阵列上执行 SVD 的第一个障碍。我最终需要在大小为 1048576x20 的数组上运行 SVD。

目前，该算法适用于大小为 10000x20 的数组，但当我转到 50000x20 时失败。

这个问题似乎源于gesvdj调用本身。调用 gesvdj 后的同步调用失败并返回泛型访问错误。

如果我使用 cuda-memcheck 运行程序，我会在同一块中的不同线程收到一系列以下错误：

Invalid __global__ write of size 4
=========     at 0x00000108 in void eye_kernel<float, int=5, int=3>(int, int, float*, int)
=========     by thread (0,7,0) in block (16,1342,0)
=========     Address 0x7f4ed40414c0 is out of bounds
=========     Saved host backtrace up to driver entry point at kernel launch time
=========     Host Frame:/lib64/libcuda.so.1 (cuLaunchKernel + 0x2c5) [0x269e85]
=========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x100c822]
=========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x100ca17]
=========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x1040dd5]
=========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x235c5d]
=========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 (cusolverDnSgesvdj + 0x508) [0x21f9a8]
=========     Host Frame:./gesvdj_example [0x4518]
=========     Host Frame:/lib64/libc.so.6 (__libc_start_main + 0xf5) [0x223d5]
=========     Host Frame:./gesvdj_example [0x3ab9]

我想知道我是否遇到了某种客户内部限制？有人有什么想法吗？如有必要，我可以提供我的确切代码，但它与示例非常相似，以至于我想我只会将人们指向那里。

谢谢！

编辑以从我链接到的示例中添加有问题的代码，算法在 assert(CUSOLVER_STATUS_SUCCESS == 状态(行处失败。我真的很陌生地使用 C 和 CUDA 编码，所以如果我遗漏了一些明显的调试信息，我很抱歉。

/* step 5: compute SVD */
status = cusolverDnDgesvdj(
cusolverH,
jobz,  /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
/* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
econ,  /* econ = 1 for economy size */
m,     /* nubmer of rows of A, 0 <= m */
n,     /* number of columns of A, 0 <= n  */
d_A,   /* m-by-n */
lda,   /* leading dimension of A */
d_S,   /* min(m,n)  */
/* the singular values in descending order */
d_U,   /* m-by-m if econ = 0 */
/* m-by-min(m,n) if econ = 1 */
lda,   /* leading dimension of U, ldu >= max(1,m) */
d_V,   /* n-by-n if econ = 0  */
/* n-by-min(m,n) if econ = 1  */
lda,   /* leading dimension of V, ldv >= max(1,n) */
d_work,
lwork,
d_info,
gesvdj_params);
cudaStat1 = cudaDeviceSynchronize();
assert(CUSOLVER_STATUS_SUCCESS == status);
assert(cudaSuccess == cudaStat1);

编辑 2添加我的代码...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <cusolverDn.h>
#include<iostream>
#include<iomanip>
#include<assert.h>
#include <cuda_runtime_api.h>
#include "cuda_runtime.h"
int main(int argc, char*argv[])
{
//Setting which device to run on
const int device = 0;
cudaSetDevice(device);
cusolverDnHandle_t cusolverH = NULL;
cudaStream_t stream = NULL;
gesvdjInfo_t gesvdj_params = NULL;
cusolverStatus_t status = CUSOLVER_STATUS_SUCCESS;
cudaError_t cudaStat1 = cudaSuccess;
cudaError_t cudaStat2 = cudaSuccess;
cudaError_t cudaStat3 = cudaSuccess;
cudaError_t cudaStat4 = cudaSuccess;
cudaError_t cudaStat5 = cudaSuccess;
const long int m = 50000;
const long int n = 20;
const int lda = m;
// --- Setting the host, Nrows x Ncols matrix
float *A = (float *)malloc(m * n * sizeof(float));
for(long int j = 0; j < m; j++)
for(long int i = 0; i < n; i++)
A[j + i*m] = sqrt((float)(i + j));
float *U = (float *)malloc(m * m * sizeof(float)); /* m-by-m unitary matrix, left singular vectors  */
float *V = (float *)malloc(m * n * sizeof(float)); /* n-by-n unitary matrix, right singular vectors */
float S[n];     /* numerical singular value */ 
float *d_A = NULL;  /* device copy of A */
float *d_S = NULL;  /* singular values */
float *d_U = NULL;  /* left singular vectors */
float *d_V = NULL;  /* right singular vectors */
int *d_info = NULL;  /* error info */
int lwork = 0;       /* size of workspace */
float *d_work = NULL; /* devie workspace for gesvdj */
int info = 0;        /* host copy of error info */
/* configuration of gesvdj  */
const double tol = 1.e-7;
const int max_sweeps = 100;
const cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvectors.
const int econ = 0 ; /* econ = 1 for economy size */
/* numerical results of gesvdj  */
double residual = 0;
int executed_sweeps = 0;
printf("tol = %E, default value is machine zero n", tol);
printf("max. sweeps = %d, default value is 100n", max_sweeps);
printf("econ = %d n", econ);
printf("=====n");
/* step 1: create cusolver handle, bind a stream */
status = cusolverDnCreate(&cusolverH);
assert(CUSOLVER_STATUS_SUCCESS == status);
cudaStat1 = cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);
assert(cudaSuccess == cudaStat1);
status = cusolverDnSetStream(cusolverH, stream);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* step 2: configuration of gesvdj */
status = cusolverDnCreateGesvdjInfo(&gesvdj_params);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* default value of tolerance is machine zero */
status = cusolverDnXgesvdjSetTolerance(
gesvdj_params,
tol);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* default value of max. sweeps is 100 */
status = cusolverDnXgesvdjSetMaxSweeps(
gesvdj_params,
max_sweeps);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* step 3: copy A to device */
cudaStat1 = cudaMalloc ((void**)&d_A   , sizeof(float)*lda*n);
cudaStat2 = cudaMalloc ((void**)&d_S   , sizeof(float)*n);
cudaStat3 = cudaMalloc ((void**)&d_U   , sizeof(float)*lda*m);
cudaStat4 = cudaMalloc ((void**)&d_V   , sizeof(float)*lda*n);
cudaStat5 = cudaMalloc ((void**)&d_info, sizeof(int));
assert(cudaSuccess == cudaStat1);
assert(cudaSuccess == cudaStat2);
assert(cudaSuccess == cudaStat3);
assert(cudaSuccess == cudaStat4);
assert(cudaSuccess == cudaStat5);
cudaStat1 = cudaMemcpy(d_A, A, sizeof(float)*lda*n, cudaMemcpyHostToDevice);
assert(cudaSuccess == cudaStat1);
/* step 4: query workspace of SVD */
status = cusolverDnSgesvdj_bufferSize(
cusolverH,
jobz, /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
/* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
econ, /* econ = 1 for economy size */
m,    /* nubmer of rows of A, 0 <= m */
n,    /* number of columns of A, 0 <= n  */
d_A,  /* m-by-n */
lda,  /* leading dimension of A */
d_S,  /* min(m,n) */
/* the singular values in descending order */
d_U,  /* m-by-m if econ = 0 */
/* m-by-min(m,n) if econ = 1 */
lda,  /* leading dimension of U, ldu >= max(1,m) */
d_V,  /* n-by-n if econ = 0  */
/* n-by-min(m,n) if econ = 1  */
lda,  /* leading dimension of V, ldv >= max(1,n) */
&lwork,
gesvdj_params);
assert(CUSOLVER_STATUS_SUCCESS == status);
cudaStat1 = cudaMalloc((void**)&d_work , sizeof(float)*lwork);
assert(cudaSuccess == cudaStat1);
/* step 5: compute SVD */
//Iterating over SVD calculation, not part of example
int iters;
for (iters = 10; iters > 0; iters--){
status = cusolverDnSgesvdj(
cusolverH,
jobz,  /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
/* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
econ,  /* econ = 1 for economy size */
m,     /* nubmer of rows of A, 0 <= m */
n,     /* number of columns of A, 0 <= n  */
d_A,   /* m-by-n */
lda,   /* leading dimension of A */
d_S,   /* min(m,n)  */
/* the singular values in descending order */
d_U,   /* m-by-m if econ = 0 */
/* m-by-min(m,n) if econ = 1 */
lda,   /* leading dimension of U, ldu >= max(1,m) */
d_V,   /* n-by-n if econ = 0  */
/* n-by-min(m,n) if econ = 1  */
lda,   /* leading dimension of V, ldv >= max(1,n) */
d_work,
lwork,
d_info,
gesvdj_params);
cudaStat1 = cudaDeviceSynchronize();
assert(CUSOLVER_STATUS_SUCCESS == status);
assert(cudaSuccess == cudaStat1);
cudaStat1 = cudaMemcpy(U, d_U, sizeof(float)*lda*m, cudaMemcpyDeviceToHost);
cudaStat2 = cudaMemcpy(V, d_V, sizeof(float)*lda*n, cudaMemcpyDeviceToHost);
cudaStat3 = cudaMemcpy(S, d_S, sizeof(float)*n    , cudaMemcpyDeviceToHost);
cudaStat4 = cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost);
cudaStat5 = cudaDeviceSynchronize();
assert(cudaSuccess == cudaStat1);
assert(cudaSuccess == cudaStat2);
assert(cudaSuccess == cudaStat3);
assert(cudaSuccess == cudaStat4);
assert(cudaSuccess == cudaStat5);
if ( 0 == info ){
printf("gesvdj converges n");
}else if ( 0 > info ){
printf("%d-th parameter is wrong n", -info);
exit(1);
}else{
printf("WARNING: info = %d : gesvdj does not converge n", info );
}
}
/* step 6: measure error of singular value */
status = cusolverDnXgesvdjGetSweeps(
cusolverH,
gesvdj_params,
&executed_sweeps);
assert(CUSOLVER_STATUS_SUCCESS == status);
status = cusolverDnXgesvdjGetResidual(
cusolverH,
gesvdj_params,
&residual);
assert(CUSOLVER_STATUS_SUCCESS == status);
printf("residual |A - U*S*V**H|_F = %E n", residual );
printf("number of executed sweeps = %d nn", executed_sweeps );
/*  free resources  */
if (A      ) free(A);
if (V      ) free(V);
if (U      ) free(U);
if (d_A    ) cudaFree(d_A);
if (d_S    ) cudaFree(d_S);
if (d_U    ) cudaFree(d_U);
if (d_V    ) cudaFree(d_V);
if (d_info ) cudaFree(d_info);
if (d_work ) cudaFree(d_work);
if (cusolverH    ) cusolverDnDestroy(cusolverH);
if (stream       ) cudaStreamDestroy(stream);
if (gesvdj_params) cusolverDnDestroyGesvdjInfo(gesvdj_params);
cudaDeviceReset();
return 0;
}

感谢您@talonmies在诊断问题方面的帮助。 Cusolver 的 Gesvdj 方法具有经济模式，它将 U 和 V 矩阵存储在更经济的数组中。我为使代码正常工作而进行的修改很简单。

econ = 1
U array size (mxn)
V array size (nxn)
ldv paramater = n

代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <cusolverDn.h>
#include<iostream>
#include<iomanip>
#include<assert.h>
#include <cuda_runtime_api.h>
#include "cuda_runtime.h"
int main(int argc, char*argv[])
{
//Setting which device to run on
const int device = 0;
cudaSetDevice(device);
cusolverDnHandle_t cusolverH = NULL;
cudaStream_t stream = NULL;
gesvdjInfo_t gesvdj_params = NULL;
cusolverStatus_t status = CUSOLVER_STATUS_SUCCESS;
cudaError_t cudaStat1 = cudaSuccess;
cudaError_t cudaStat2 = cudaSuccess;
cudaError_t cudaStat3 = cudaSuccess;
cudaError_t cudaStat4 = cudaSuccess;
cudaError_t cudaStat5 = cudaSuccess;
const long int m = 1048576;
const long int n = 20;
const int lda = m;
// --- Setting the host, Nrows x Ncols matrix
float *A = (float *)malloc(m * n * sizeof(float));
for(long int j = 0; j < m; j++)
for(long int i = 0; i < n; i++)
A[j + i*m] = sqrt((float)(i + j));
float *U = (float *)malloc(m * n * sizeof(float)); /* m-by-m unitary matrix, left singular vectors  */
float *V = (float *)malloc(n * n * sizeof(float)); /* n-by-n unitary matrix, right singular vectors */
float S[n];     /* numerical singular value */ 
float *d_A = NULL;  /* device copy of A */
float *d_S = NULL;  /* singular values */
float *d_U = NULL;  /* left singular vectors */
float *d_V = NULL;  /* right singular vectors */
int *d_info = NULL;  /* error info */
int lwork = 0;       /* size of workspace */
float *d_work = NULL; /* devie workspace for gesvdj */
int info = 0;        /* host copy of error info */
/* configuration of gesvdj  */
const double tol = 1.e-7;
const int max_sweeps = 100;
const cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvectors.
const int econ = 1 ; /* econ = 1 for economy size */
/* numerical results of gesvdj  */
double residual = 0;
int executed_sweeps = 0;
printf("tol = %E, default value is machine zero n", tol);
printf("max. sweeps = %d, default value is 100n", max_sweeps);
printf("econ = %d n", econ);
printf("=====n");
/* step 1: create cusolver handle, bind a stream */
status = cusolverDnCreate(&cusolverH);
assert(CUSOLVER_STATUS_SUCCESS == status);
cudaStat1 = cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);
assert(cudaSuccess == cudaStat1);
status = cusolverDnSetStream(cusolverH, stream);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* step 2: configuration of gesvdj */
status = cusolverDnCreateGesvdjInfo(&gesvdj_params);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* default value of tolerance is machine zero */
status = cusolverDnXgesvdjSetTolerance(
gesvdj_params,
tol);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* default value of max. sweeps is 100 */
status = cusolverDnXgesvdjSetMaxSweeps(
gesvdj_params,
max_sweeps);
assert(CUSOLVER_STATUS_SUCCESS == status);
/* step 3: copy A to device */
cudaStat1 = cudaMalloc ((void**)&d_A   , sizeof(float)*lda*n);
cudaStat2 = cudaMalloc ((void**)&d_S   , sizeof(float)*n);
cudaStat3 = cudaMalloc ((void**)&d_U   , sizeof(float)*lda*n);
cudaStat4 = cudaMalloc ((void**)&d_V   , sizeof(float)*n*n);
cudaStat5 = cudaMalloc ((void**)&d_info, sizeof(int));
assert(cudaSuccess == cudaStat1);
assert(cudaSuccess == cudaStat2);
assert(cudaSuccess == cudaStat3);
assert(cudaSuccess == cudaStat4);
assert(cudaSuccess == cudaStat5);
cudaStat1 = cudaMemcpy(d_A, A, sizeof(float)*lda*n, cudaMemcpyHostToDevice);
assert(cudaSuccess == cudaStat1);
/* step 4: query workspace of SVD */
status = cusolverDnSgesvdj_bufferSize(
cusolverH,
jobz, /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
/* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
econ, /* econ = 1 for economy size */
m,    /* nubmer of rows of A, 0 <= m */
n,    /* number of columns of A, 0 <= n  */
d_A,  /* m-by-n */
lda,  /* leading dimension of A */
d_S,  /* min(m,n) */
/* the singular values in descending order */
d_U,  /* m-by-m if econ = 0 */
/* m-by-min(m,n) if econ = 1 */
lda,  /* leading dimension of U, ldu >= max(1,m) */
d_V,  /* n-by-n if econ = 0  */
/* n-by-min(m,n) if econ = 1  */
n,  /* leading dimension of V, ldv >= max(1,n) */
&lwork,
gesvdj_params);
assert(CUSOLVER_STATUS_SUCCESS == status);
cudaStat1 = cudaMalloc((void**)&d_work , sizeof(float)*lwork);
assert(cudaSuccess == cudaStat1);
/* step 5: compute SVD */
//Iterating over SVD calculation, not part of example
int iters;
for (iters = 10; iters > 0; iters--){ 
status = cusolverDnSgesvdj(
cusolverH,
jobz,  /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
/* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
econ,  /* econ = 1 for economy size */
m,     /* nubmer of rows of A, 0 <= m */
n,     /* number of columns of A, 0 <= n  */
d_A,   /* m-by-n */
lda,   /* leading dimension of A */
d_S,   /* min(m,n)  */
/* the singular values in descending order */
d_U,   /* m-by-m if econ = 0 */
/* m-by-min(m,n) if econ = 1 */
lda,   /* leading dimension of U, ldu >= max(1,m) */
d_V,   /* n-by-n if econ = 0  */
/* n-by-min(m,n) if econ = 1  */
n,   /* leading dimension of V, ldv >= max(1,n) */
d_work,
lwork,
d_info,
gesvdj_params);
cudaStat1 = cudaDeviceSynchronize();
assert(CUSOLVER_STATUS_SUCCESS == status);
assert(cudaSuccess == cudaStat1);
cudaStat1 = cudaMemcpy(U, d_U, sizeof(float)*lda*n, cudaMemcpyDeviceToHost);
cudaStat2 = cudaMemcpy(V, d_V, sizeof(float)*n*n, cudaMemcpyDeviceToHost);
cudaStat3 = cudaMemcpy(S, d_S, sizeof(float)*n    , cudaMemcpyDeviceToHost);
cudaStat4 = cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost);
cudaStat5 = cudaDeviceSynchronize();
assert(cudaSuccess == cudaStat1);
assert(cudaSuccess == cudaStat2);
assert(cudaSuccess == cudaStat3);
assert(cudaSuccess == cudaStat4);
assert(cudaSuccess == cudaStat5);
if ( 0 == info ){
printf("gesvdj converges n");
}else if ( 0 > info ){
printf("%d-th parameter is wrong n", -info);
exit(1);
}else{
printf("WARNING: info = %d : gesvdj does not converge n", info );
}
}
/* step 6: measure error of singular value */
status = cusolverDnXgesvdjGetSweeps(
cusolverH,
gesvdj_params,
&executed_sweeps);
assert(CUSOLVER_STATUS_SUCCESS == status);
status = cusolverDnXgesvdjGetResidual(
cusolverH,
gesvdj_params,
&residual);
assert(CUSOLVER_STATUS_SUCCESS == status);
printf("residual |A - U*S*V**H|_F = %E n", residual );
printf("number of executed sweeps = %d nn", executed_sweeps );
/*  free resources  */
if (A      ) free(A);
if (V      ) free(V);
if (U      ) free(U);
if (d_A    ) cudaFree(d_A);
if (d_S    ) cudaFree(d_S);
if (d_U    ) cudaFree(d_U);
if (d_V    ) cudaFree(d_V);
if (d_info ) cudaFree(d_info);
if (d_work ) cudaFree(d_work);
if (cusolverH    ) cusolverDnDestroy(cusolverH);
if (stream       ) cudaStreamDestroy(stream);
if (gesvdj_params) cusolverDnDestroyGesvdjInfo(gesvdj_params);
cudaDeviceReset();
return 0;
}

相关内容

  • 没有找到相关文章

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium