Advertisement
desdemona

dddd

Apr 14th, 2015
686
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 4.20 KB | None | 0 0
  1. void axpy (int n, float alpha, const float *x,int incx, float *y, int incy)
  2. {
  3.  
  4. }
  5.  
  6. void hello(char *a, int *b)
  7. {
  8.     a[threadIdx.x] += b[threadIdx.x];
  9. }
  10.  
  11.   int block_size = 4;
  12.   int n_blocks = N/block_size + (N%block_size == 0 ? 0:1);
  13.  
  14. 8DnDBBNu
  15.  
  16.     int threadsPerBlock = 256;
  17.     int blocksPerGrid =(N + threadsPerBlock - 1) / threadsPerBlock;
  18.  
  19.     // sparse matrix vector product: d_Ax = A * d_x
  20.     cusparseScsrmv(cusparseHandle,CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, nz, &alpha, descr, d_val, d_row, d_col, d_x, &beta, d_Ax);  // PODMIEN FUNCKJE (ZADANIE-II)
  21.    // vectorMatVec<<<blocksPerGrid, threadsPerBlock>>>(d_Ax, d_x, d_row, d_col, d_val, N);
  22.  
  23.  
  24.     //azpy: d_r = d_r + alpham1 * d_Ax
  25.     //cublasSaxpy(cublasHandle, N, &alpham1, d_Ax, 1, d_r, 1);                              // PODMIEN FUNCKJE (ZADANIE-I)
  26.     vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_r,d_Ax,alpham1,N);
  27.     //dot:  r1 = d_r * d_r
  28.     cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1);                        // PODMIEN FUNCKJE (ZADANIE-III)
  29.  
  30.     k = 1;
  31.         printf("My functions\n"
  32.     while (r1 > tol*tol && k <= max_iter)
  33.     {
  34.         if (k > 1)
  35.         {
  36.             b = r1 / r0;
  37.             //scal: d_p = b * d_p
  38.             //cublasStatus = cublasSscal(cublasHandle, N, &b, d_p, 1);                        // PODMIEN FUNCKJE (ZADANIE-I)
  39.                 vectorScal<<<blocksPerGrid, threadsPerBlock>>>(d_p, b, N);  
  40.             //axpy:  d_p = d_p + alpha * d_r
  41.             cublasStatus =vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_p,d_r, alpham1,N);         // PODMIEN FUNCKJE (ZADANIE-I)
  42.         //    vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_p,d_r, alpha,N);  
  43.         }
  44.         else
  45.         {
  46.             //cpy: d_p = d_r
  47.             //cublasStatus = cublasScopy(cublasHandle, N, d_r, 1, d_p, 1);                    // PODMIEN FUNCKJE (ZADANIE-I)
  48.                 vectorCopy<<<blocksPerGrid, threadsPerBlock>>>(d_p,d_r, N);  
  49.         }
  50.  
  51.         //sparse matrix-vector product: d_Ax = A * d_p
  52.         cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, nz, &alpha, descr, d_val, d_row, d_col, d_p, &beta, d_Ax); // PODMIEN FUNCKJE (ZADANIE-II)
  53.        
  54.         cublasStatus = cublasSdot(cublasHandle, N, d_p, 1, d_Ax, 1, &dot);                  // PODMIEN FUNCKJE (ZADANIE-III)
  55.         a = r1 / dot;
  56.  
  57.         //axpy: d_x = d_x + a*d_p
  58.        // cublasStatus = vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_x,d_p,a,N);                 // PODMIEN FUNCKJE (ZADANIE-I)
  59.         vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_x,d_p,a,N);    
  60.         na = -a;
  61.          
  62.         //axpy:  d_r = d_r + na * d_Ax
  63.         //cublasStatus = vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_r,d_Ax, na ,N);                 // PODMIEN FUNCKJE (ZADANIE-I)
  64.         vectorAxpy<<<blocksPerGrid, threadsPerBlock>>>(d_r,d_Ax, na ,N);  
  65.         r0 = r1;
  66.        
  67.         //dot: r1 = d_r * d_r
  68.         cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1);                    // PODMIEN FUNCKJE (ZADANIE-III)
  69.         cudaThreadSynchronize();
  70.         printf("iteration = %3d, residual = %e\n", k, sqrt(r1));
  71.         k++;
  72.     }
  73.  
  74.     cudaMemcpy(x, d_x, N*sizeof(float), cudaMemcpyDeviceToHost);
  75.  
  76.     float rsum, diff, err = 0.0;
  77.  
  78.     for (int i = 0; i < N; i++)
  79.     {
  80.         rsum = 0.0;
  81.  
  82.         for (int j = I[i]; j < I[i+1]; j++)
  83.         {
  84.             rsum += val[j]*x[J[j]];
  85.         }
  86.  
  87.         diff = fabs(rsum - rhs[i]);
  88.  
  89.         if (diff > err)
  90.         {
  91.             err = diff;
  92.         }
  93.     }
  94.  
  95.     cusparseDestroy(cusparseHandle);
  96.     cublasDestroy(cublasHandle);
  97.  
  98.     free(I);
  99.     free(J);
  100.     free(val);
  101.     free(x);
  102.     free(rhs);
  103.     cudaFree(d_col);
  104.     cudaFree(d_row);
  105.     cudaFree(d_val);
  106.     cudaFree(d_x);
  107.     cudaFree(d_r);
  108.     cudaFree(d_p);
  109.     cudaFree(d_Ax);
  110.  
  111.     cudaDeviceReset();
  112.  
  113.     printf("Test Summary:  Error amount = %e\n", err);
  114.     //exit((k <= max_iter) ? 0 : 1);
  115.  
  116.  
  117. }
  118.  
  119. int main(int argc, char **argv)
  120. {
  121.     //int N = 1e6;//1 << 20;
  122.     //int N = 256 * (1<<10)  -10 ; //1e6;//1 << 20;
  123.     int N = 1e5;
  124.     int M = N;
  125.    
  126.     cgs_basic(argc, argv, N, M);
  127.    
  128.     cgs_TODO(argc, argv, N, M);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement