Guest User

Untitled

a guest
Feb 26th, 2013
171
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <cusp/hyb_matrix.h>
  2. #include <cusp/dia_matrix.h>
  3. #include <cusp/csr_matrix.h>
  4. #include <cusp/gallery/poisson.h>
  5. #include <cusp/krylov/gmres.h>
  6. #include <cusp/io/matrix_market.h>
  7. #include <cusp/print.h>
  8. #include <cstdlib>
  9. #include <iostream>
  10. #include <ctime>
  11.  
  12.  
  13.  
  14. struct event_pair
  15. {
  16.   cudaEvent_t start;
  17.   cudaEvent_t end;
  18. };
  19.  
  20. /*Start clock used for timing kernels*/
  21. inline void start_timer(event_pair * p)
  22. {
  23.   cudaEventCreate(&p->start);
  24.   cudaEventCreate(&p->end);
  25.   cudaEventRecord(p->start, 0);
  26. }
  27.  
  28. /*Stop clock used for timing kernels*/
  29. inline void stop_timer(event_pair * p, char * kernel_name)
  30. {
  31.   cudaEventRecord(p->end, 0);
  32.   cudaEventSynchronize(p->end);
  33.  
  34.   float elapsed_time;
  35.   cudaEventElapsedTime(&elapsed_time, p->start, p->end);
  36.   printf("%s took %.1f ms\n",kernel_name, elapsed_time);
  37.   cudaEventDestroy(p->start);
  38.   cudaEventDestroy(p->end);
  39. }
  40.  
  41. /////////////////////////////////////////////////////////////////
  42. // where to perform the computation
  43. typedef cusp::device_memory MemorySpace;
  44.  
  45. // which floating point type to use
  46. typedef double ValueType;
  47.  
  48. int main(int argc, char** argv)
  49. {
  50.     // create an empty sparse matrix structure (HYB format)
  51.   cusp::dia_matrix<int, ValueType, MemorySpace> A;
  52.  
  53.     // load a matrix stored in MatrixMarket format
  54.     cusp::io::read_matrix_market_file(A, argv[1]);
  55.  
  56.     // allocate storage for solution (x) and right hand side (b)
  57.     cusp::array1d<ValueType, MemorySpace> x(A.num_rows, ValueType(1));
  58.     cusp::array1d<ValueType, MemorySpace> b(A.num_rows);
  59.  
  60.     cusp::multiply(A,x,b);
  61.  
  62.     // set initial guess
  63.     thrust::fill( x.begin(), x.end(), ValueType(2.0) );
  64.  
  65.     // set stopping criteria:
  66.     //  iteration_limit    = 100
  67.     //  relative_tolerance = 1e-6
  68.     int iteration_limit = atoi(argv[2]);
  69.     int relative_tolerance = 1e-20;
  70.     cusp::default_monitor<ValueType> monitor(b, iteration_limit, relative_tolerance);
  71.     int restart = 30;
  72.  
  73.  
  74.         event_pair timer;
  75.         // set initial guess
  76.        thrust::fill( x.begin(), x.end(), ValueType(2.0) );
  77.        start_timer(&timer);
  78.        cusp::krylov::gmres(A,x,b,restart,monitor);
  79.        cudaThreadSynchronize();
  80.        stop_timer(&timer,"GMRES solver");
  81.        
  82.  
  83.   // report solver results
  84.       if (monitor.converged())
  85.       {
  86.           std::cout << "Solver converged to " << monitor.relative_tolerance() << " relative tolerance";
  87.           std::cout << " after " << monitor.iteration_count() << " iterations" << std::endl;
  88.       }
  89.       else
  90.       {
  91.           std::cout << "Solver reached iteration limit " << monitor.iteration_limit() << " before converging";
  92.           std::cout << " to " << monitor.relative_tolerance() << " relative tolerance " << std::endl;
  93.       }
  94.  
  95.  
  96.      
  97.      std::cout << "Size of the matrix A is " << A.num_rows << "  " << A.num_cols << std::endl;  
  98.    
  99.     cusp::print(x);
  100.     return 0;
  101. }
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×