找回密码
 立即注册
发表于 2021-6-2 20:54:28
4971
我搭建了一个平台后在github上找了个例子试运行运行不了,我已经加载了cublas和cusparse这两个需要加载的库了.希望有人可以帮忙解答问题,因为我附近的人没人学过这个

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
使用道具 举报 回复
发表于 2021-6-3 09:49:15
我在这个网站上找到的代码想去跑跑,这个代码应该在linux系统中编写的
  1. #include <cuda_runtime.h>
  2. #include <stdio.h>
  3. #include "cusparse.h"
  4. #include "cublas_v2.h"
  5. #include <time.h>
  6. #ifdef _WIN32
  7. #include <windows.h>
  8. #else
  9. #include <sys/time.h>
  10. #endif

  11. #define CLEANUP(s){\
  12.         printf("%s\n", s);\
  13.         if(h_valA) free(h_valA);\
  14.         if(h_csrRowPtrA) free(h_csrRowPtrA);\
  15.         if(h_csrColIndA) free(h_csrColIndA);\
  16.         if(valA) cudaFree(valA);\
  17.         if(csrRowPtrA) cudaFree(csrRowPtrA);\
  18.         if(csrColIndA) cudaFree(csrColIndA);\
  19.         if(h_b) free(h_b);\
  20.         if(b) cudaFree(b);\
  21.         if(r) cudaFree(r);\
  22.         if(p) cudaFree(p);\
  23.         if(Ap) cudaFree(Ap);\
  24.         if(h_x) free(h_x);\
  25.         if(x) cudaFree(x);\
  26.         if(handle1) cublasDestroy(handle1);\
  27.         if(descrA) cusparseDestroyMatDescr(descrA);\
  28.         if(handle2) cusparseDestroy(handle2);\
  29.         cudaDeviceReset();\
  30. }

  31. //用于统计程序时间
  32. double cpuSecond(){
  33.         struct timeval tp;
  34.         gettimeofday(&tp, NULL);
  35.         return ((double)tp.tv_sec + (double)tp.tv_usec * 1e-6);
  36. }

  37. int main(){
  38.         //initialization
  39.         cudaError_t cudaStat1, cudaStat2, cudaStat3;       

  40.         double* h_valA = 0;
  41.         int* h_csrRowPtrA = 0;
  42.         int* h_csrColIndA = 0;
  43.         double* valA = 0;
  44.         int* csrRowPtrA = 0;
  45.         int* csrColIndA = 0;
  46.         double* h_b = 0;
  47.         double* b = 0;
  48.         double* r = 0;
  49.         double* p = 0;
  50.         double* Ap = 0;
  51.         double* h_x = 0;
  52.         double* x = 0;

  53.         cublasStatus_t status1;
  54.         cublasHandle_t handle1 = 0;
  55.         cusparseStatus_t status2;
  56.         cusparseHandle_t handle2 = 0;
  57.         cusparseMatDescr_t descrA = 0;

  58.         int n, nnz;
  59.         double dzero = 0.0;
  60.         double done = 1.0;
  61.         /*
  62.         use CG to solve following linear system:
  63.                 Ax = b;
  64.                 A = |3 0 2|
  65.                         |0 2 0|
  66.                         |2 0 1|
  67.                 b = |3.5|
  68.                         |1.5|
  69.                         |2.0|
  70.         correct result x:
  71.                 x = |0.50|
  72.                         |0.75|
  73.                         |1.00|
  74.            */
  75.         n = 3;
  76.         nnz = 5;
  77.         //which device(gpu) to use
  78.         cudaStat1 = cudaSetDevice(0);
  79.         if(cudaStat1 != cudaSuccess){
  80.                 CLEANUP("Device Set failed");
  81.                 return EXIT_FAILURE;
  82.         }

  83.         //matrix A(csr format) in host
  84.         h_valA = (double*)malloc(nnz*sizeof(double));
  85.         h_csrRowPtrA = (int*)malloc((n+1)*sizeof(int));
  86.         h_csrColIndA = (int*)malloc(nnz*sizeof(int));
  87.         if(!h_valA||!h_csrRowPtrA||!h_csrColIndA){
  88.                 CLEANUP("Host malloc failed(A)");
  89.                 return EXIT_FAILURE;
  90.         }

  91.         h_valA[0] = 3.0;
  92.         h_valA[1] = 2.0;
  93.         h_valA[2] = 2.0;
  94.         h_valA[3] = 2.0;
  95.         h_valA[4] = 1.0;
  96.        
  97.         h_csrRowPtrA[0] = 0;
  98.         h_csrRowPtrA[1] = 2;
  99.         h_csrRowPtrA[2] = 3;
  100.         h_csrRowPtrA[3] = 5;
  101.        
  102.         h_csrColIndA[0] = 0;
  103.         h_csrColIndA[1] = 2;
  104.         h_csrColIndA[2] = 1;
  105.         h_csrColIndA[3] = 0;
  106.         h_csrColIndA[4] = 2;

  107.         //matrix A(csr format) in device
  108.         cudaStat1 = cudaMalloc((void**)&valA, nnz*sizeof(double));
  109.         cudaStat2 = cudaMalloc((void**)&csrRowPtrA, (n+1)*sizeof(int));
  110.         cudaStat3 = cudaMalloc((void**)&csrColIndA, nnz*sizeof(int));
  111.         if(cudaStat1 != cudaSuccess || cudaStat2 != cudaSuccess || cudaStat3 != cudaSuccess){
  112.                 CLEANUP("device malloc failed(A)");
  113.                 return EXIT_FAILURE;
  114.         }
  115.         cudaStat1 = cudaMemcpy(valA, h_valA, (size_t)(nnz*sizeof(double)), cudaMemcpyHostToDevice);
  116.         cudaStat2 = cudaMemcpy(csrRowPtrA, h_csrRowPtrA, (size_t)((n+1)*sizeof(int)), cudaMemcpyHostToDevice);
  117.         cudaStat3 = cudaMemcpy(csrColIndA, h_csrColIndA, (size_t)(nnz*sizeof(int)), cudaMemcpyHostToDevice);
  118.         if(cudaStat1 != cudaSuccess || cudaStat2 != cudaSuccess || cudaStat3 != cudaSuccess){
  119.                 CLEANUP("memcpy from host to device failed(A)");
  120.                 return EXIT_FAILURE;
  121.         }

  122.         //vector b in host
  123.         h_b = (double*)malloc(n*sizeof(double));
  124.         h_b[0] = 3.5;
  125.         h_b[1] = 1.5;
  126.         h_b[2] = 2.0;

  127.         //vector b in device
  128.         cudaStat1 = cudaMalloc((void**)&b, n*sizeof(double));
  129.         if(cudaStat1 != cudaSuccess){
  130.                 CLEANUP("device malloc failed(b)");
  131.                 return EXIT_FAILURE;
  132.         }
  133.         cudaStat1 = cudaMemcpy(b, h_b, (size_t)(n*sizeof(double)), cudaMemcpyHostToDevice);
  134.         if(cudaStat1 != cudaSuccess){
  135.                 CLEANUP("memcpy from host to device failed(b)");
  136.                 return EXIT_FAILURE;
  137.         }

  138.         //vector r in device
  139.         cudaStat1 = cudaMalloc((void**)&r, n*sizeof(double));
  140.         if(cudaStat1 != cudaSuccess){
  141.                 CLEANUP("device malloc failed(r)");
  142.                 return EXIT_FAILURE;
  143.         }

  144.         //vector p in device
  145.         cudaStat1 = cudaMalloc((void**)&p, n*sizeof(double));
  146.         if(cudaStat1 != cudaSuccess){
  147.                 CLEANUP("device malloc failed(p)");
  148.                 return EXIT_FAILURE;
  149.         }

  150.         //vector Ap in device
  151.         cudaStat1 = cudaMalloc((void**)&Ap, n*sizeof(double));
  152.         if(cudaStat1 != cudaSuccess){
  153.                 CLEANUP("device malloc failed(Ap)");
  154.                 return EXIT_FAILURE;
  155.         }

  156.         //vector x in host
  157.         h_x = (double*)malloc(n*sizeof(double));
  158.         for(int i=0; i<n; ++i)
  159.                 h_x[i] = 0.0;

  160.         //vector x in device
  161.         cudaStat1 = cudaMalloc((void**)&x, n*sizeof(double));
  162.         if(cudaStat1 != cudaSuccess){
  163.                 CLEANUP("device malloc faile(x)");
  164.                 return EXIT_FAILURE;
  165.         }
  166.         cudaStat1 = cudaMemcpy(x, h_x, (size_t)(n*sizeof(double)), cudaMemcpyHostToDevice);
  167.         if(cudaStat1 != cudaSuccess){
  168.                 CLEANUP("memcpy from host to device failed(x)");
  169.                 return EXIT_FAILURE;
  170.         }

  171.         //initialize cublas library
  172.         status1 = cublasCreate(&handle1);
  173.         if(status1 != CUBLAS_STATUS_SUCCESS){
  174.                 CLEANUP("CUBLAS initialize failed");
  175.                 return EXIT_FAILURE;
  176.         }
  177.         //initialize cusparse library
  178.         status2 = cusparseCreate(&handle2);
  179.         if(status2 != CUSPARSE_STATUS_SUCCESS){
  180.                 CLEANUP("CUSPARSE Library initialize failed");
  181.                 return EXIT_FAILURE;
  182.         }
  183.         //initialize matrix descriptor
  184.         status2 = cusparseCreateMatDescr(&descrA);
  185.         if(status2 != CUSPARSE_STATUS_SUCCESS){
  186.                 CLEANUP("Matrix descriptor initialize failed");
  187.                 return EXIT_FAILURE;
  188.         }
  189.         cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL);
  190.         cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ZERO);


  191.         //CG
  192.         //1:
  193.         //max iteration number = 2000
  194.         //relative residual = 1e-7
  195.         //compute initial residual r = b - Ax_0
  196.         //x_0 = |0|
  197.         //                |0|
  198.         //                |0|
  199.         //initial p = r
  200.         int maxit = 2000;
  201.         double tol = 1e-7;
  202.         double alpha, beta, rhop, rho;
  203.         //r = b
  204.         status1 = cublasDcopy(handle1, n, b, 1, r, 1);
  205.         if(status1 != CUBLAS_STATUS_SUCCESS){
  206.                 CLEANUP("Vector copy failed(r)");
  207.                 return EXIT_FAILURE;
  208.         }
  209.         //p = b
  210.         status1 = cublasDcopy(handle1, n, b, 1, p, 1);
  211.         if(status1 != CUBLAS_STATUS_SUCCESS){
  212.                 CLEANUP("Vector copy failed(p)");
  213.                 return EXIT_FAILURE;
  214.         }
  215.         status1 = cublasDnrm2(handle1, n, r, 1, &rho);
  216.         if(status1 != CUBLAS_STATUS_SUCCESS){
  217.                 CLEANUP("compute norm2 failed(r)");
  218.                 return EXIT_FAILURE;
  219.         }
  220.         rho = rho*rho;

  221.         //2: repeat until convergence
  222.         for(int i=0; i<maxit; ++i){
  223.                 status2 = cusparseDcsrmv(handle2, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, nnz, &done,
  224.                                                                         descrA, valA, csrRowPtrA, csrColIndA, p, &dzero, Ap);
  225.                 if(status2 != CUSPARSE_STATUS_SUCCESS){
  226.                         CLEANUP("compute mv failed(Ap)");
  227.                         return EXIT_FAILURE;
  228.                 }
  229.                 double tmp;
  230.                 status1 = cublasDdot(handle1, n, p, 1, Ap, 1, &tmp);
  231.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  232.                         CLEANUP("compute dot failed(tmp)");
  233.                         return EXIT_FAILURE;
  234.                 }
  235.                 //alpha = (r^T*r)/(p^T*A*p)
  236.                 alpha = rho/tmp;
  237.                 //x = x + alpha*p
  238.                 status1 = cublasDaxpy(handle1, n, &alpha, p, 1, x, 1);
  239.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  240.                         CLEANUP("compute axpy failed(x)");
  241.                         return EXIT_FAILURE;
  242.                 }
  243.                 double tmp2 = -alpha;
  244.                 //r = r - alpha*A*p
  245.                 status1 = cublasDaxpy(handle1, n, &tmp2, Ap, 1, r, 1);
  246.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  247.                         CLEANUP("compute axpy failed(r)");
  248.                         return EXIT_FAILURE;
  249.                 }

  250.                 rhop = rho;
  251.                 status1 = cublasDnrm2(handle1, n, r, 1, &rho);
  252.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  253.                         CLEANUP("compute norm2 failed(r)");
  254.                         return EXIT_FAILURE;
  255.                 }
  256.                 if(rho < tol){
  257.                         break;
  258.                 }
  259.                 rho = rho*rho;
  260.                
  261.                 //beta = rho/rhop
  262.                 beta = rho/rhop;
  263.                
  264.                 //p = r + beta*p
  265.                 status1 = cublasDscal(handle1, n, &beta, p, 1);
  266.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  267.                         CLEANUP("compute scal failed(p)");
  268.                         return EXIT_FAILURE;
  269.                 }
  270.                 status1 = cublasDaxpy(handle1, n, &done, r, 1, p, 1);
  271.                 if(status1 != CUBLAS_STATUS_SUCCESS){
  272.                         CLEANUP("compute axpy failed(p)");
  273.                         return EXIT_FAILURE;
  274.                 }
  275.         }
  276.        
  277.         //copy x to h_x, and output h_x
  278.         cudaStat1 = cudaMemcpy(h_x, x, (size_t)(n*sizeof(double)), cudaMemcpyDeviceToHost);
  279.         if(cudaStat1 != cudaSuccess){
  280.                 CLEANUP("memcpy from device to host failed(x)");
  281.                 return EXIT_FAILURE;
  282.         }
  283.         for(int i=0; i<n; ++i)
  284.                 printf("%f\n", h_x[i]);

  285.         CLEANUP("Success!");
  286.         return 0;
  287. }
复制代码
使用道具 举报 回复 支持 反对
发新帖
您需要登录后才可以回帖 登录 | 立即注册