/[escript]/trunk/cusplibrary/cusp/detail/host/spmv.h
ViewVC logotype

Diff of /trunk/cusplibrary/cusp/detail/host/spmv.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5243 by caltinay, Thu Nov 6 10:09:23 2014 UTC revision 5244 by caltinay, Thu Nov 6 11:03:58 2014 UTC
# Line 247  void spmv_cds(const Matrix&  A, Line 247  void spmv_cds(const Matrix&  A,
247              // diagonal is processed in the second loop              // diagonal is processed in the second loop
248              const IndexType d0 = (A.diagonal_offsets[0] == 0 ? 1 : 0);              const IndexType d0 = (A.diagonal_offsets[0] == 0 ? 1 : 0);
249  #pragma omp parallel for  #pragma omp parallel for
250              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for (IndexType ch = 0; ch < num_rows; ch+=chunksize)
251              {              {
252                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row+=2)                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row++)
253                  {                  {
254                      ValueType sum1 = initialize(y[row]);                      y[row] = initialize(y[row]);
255                      ValueType sum2 = initialize(y[row+1]);                  }
256                      // process subdiagonal blocks  
257                      for(IndexType d = 0; d < num_diagonals-d0; d++)                  // process subdiagonal blocks
258                    for (IndexType d = 0; d < num_diagonals-d0; d++)
259                    {
260                        const IndexType diag = num_diagonals-d-1;
261                        const IndexType k = -2*A.diagonal_offsets[diag];
262                        for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row+=2)
263                      {                      {
264                          const IndexType diag = num_diagonals-d-1;                          const IndexType col = row + k;
                         const IndexType col = row - A.diagonal_offsets[diag]*2;  
265                          if (col >= 0 && col <= num_rows-2)                          if (col >= 0 && col <= num_rows-2)
266                          {                          {
267                              sum1 = reduce(sum1,combine(A.values(col,  2*diag),  x[col]));                              y[row]   = reduce(y[row],  combine(A.values(col,  2*diag),  x[col]));
268                              sum2 = reduce(sum2,combine(A.values(col,  2*diag+1),x[col]));                              y[row]   = reduce(y[row],  combine(A.values(col+1,2*diag),  x[col+1]));
269                              sum1 = reduce(sum1,combine(A.values(col+1,2*diag),  x[col+1]));                              y[row+1] = reduce(y[row+1],combine(A.values(col,  2*diag+1),x[col]));
270                              sum2 = reduce(sum2,combine(A.values(col+1,2*diag+1),x[col+1]));                              y[row+1] = reduce(y[row+1],combine(A.values(col+1,2*diag+1),x[col+1]));
271                          }                          }
272                      }                      }
273                      // process main and upper diagonal blocks                  }
274                      for(IndexType d = 0; d < num_diagonals; d++)                  // process main and upper diagonal blocks
275                    for (IndexType d = 0; d < num_diagonals; d++)
276                    {
277                        const IndexType k = 2*A.diagonal_offsets[d];
278                        for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row+=2)
279                      {                      {
280                          const IndexType col = row + A.diagonal_offsets[d]*2;                          const IndexType col = row + k;
281                          if (col >= 0 && col <= num_rows-2)                          if (col >= 0 && col <= num_rows-2)
282                          {                          {
283                              sum1 = reduce(sum1,combine(A.values(row,  2*d),  x[col]));                              y[row]   = reduce(y[row],  combine(A.values(row,  2*d),  x[col]));
284                              sum2 = reduce(sum2,combine(A.values(row+1,2*d),  x[col]));                              y[row+1] = reduce(y[row+1],combine(A.values(row+1,2*d),  x[col]));
285                              sum1 = reduce(sum1,combine(A.values(row,  2*d+1),x[col+1]));                              y[row]   = reduce(y[row],  combine(A.values(row,  2*d+1),x[col+1]));
286                              sum2 = reduce(sum2,combine(A.values(row+1,2*d+1),x[col+1]));                              y[row+1] = reduce(y[row+1],combine(A.values(row+1,2*d+1),x[col+1]));
287                          }                          }
288                      }                      }
                     y[row] = sum1;  
                     y[row+1] = sum2;  
289                  }                  }
290              }              }
291          } else { // A.symmetric          } else { // !A.symmetric
292  #pragma omp parallel for  #pragma omp parallel for
293              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for (IndexType ch = 0; ch < num_rows; ch+=chunksize)
294              {              {
295                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row+=2)                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row+=2)
296                  {                  {
297                      ValueType sum1 = initialize(y[row]);                      ValueType sum1 = initialize(y[row]);
298                      ValueType sum2 = initialize(y[row+1]);                      ValueType sum2 = initialize(y[row+1]);
299                      // for each diagonal block                      // for each diagonal block
300                      for(IndexType d = 0; d < num_diagonals; d++)                      for (IndexType d = 0; d < num_diagonals; d++)
301                      {                      {
302                          const IndexType col = row + A.diagonal_offsets[d]*2;                          const IndexType col = row + A.diagonal_offsets[d]*2;
303                          if (col >= 0 && col <= num_rows-2)                          if (col >= 0 && col <= num_rows-2)
# Line 316  void spmv_cds(const Matrix&  A, Line 322  void spmv_cds(const Matrix&  A,
322              const ValueType* values = thrust::raw_pointer_cast(&A.values.values[0]);              const ValueType* values = thrust::raw_pointer_cast(&A.values.values[0]);
323              const IndexType pitch = A.values.pitch;              const IndexType pitch = A.values.pitch;
324  #pragma omp parallel for  #pragma omp parallel for
325              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for (IndexType ch = 0; ch < num_rows; ch+=chunksize)
326              {              {
327                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row++)                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row++)
328                  {                  {
# Line 373  void spmv_cds(const Matrix&  A, Line 379  void spmv_cds(const Matrix&  A,
379                              }                              }
380                          }                          }
381                      }                      }
382                  }                  } // diagonals
383              }              }
384          } else { // A.symmetric          } else { // !A.symmetric
385  #pragma omp parallel for  #pragma omp parallel for
386              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for (IndexType ch = 0; ch < num_rows; ch+=chunksize)
387              {              {
388                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row++)                  for (IndexType row = ch; row<std::min(ch+chunksize,num_rows); row++)
389                  {                  {
# Line 385  void spmv_cds(const Matrix&  A, Line 391  void spmv_cds(const Matrix&  A,
391                  }                  }
392    
393                  // for each diagonal block                  // for each diagonal block
394                  for(IndexType d = 0; d < num_diagonals; d++)                  for (IndexType d = 0; d < num_diagonals; d++)
395                  {                  {
396                      const IndexType k = A.diagonal_offsets[d]*block_size;                      const IndexType k = A.diagonal_offsets[d]*block_size;
397    
398                      for(IndexType row=ch; row<std::min(ch+chunksize,num_rows); row+=block_size)                      for (IndexType row=ch; row<std::min(ch+chunksize,num_rows); row+=block_size)
399                      {                      {
400                          const IndexType col = row + k;                          const IndexType col = row + k;
401                          if (col >= 0 && col <= num_rows-block_size)                          if (col >= 0 && col <= num_rows-block_size)
# Line 408  void spmv_cds(const Matrix&  A, Line 414  void spmv_cds(const Matrix&  A,
414                              }                              }
415                          }                          }
416                      }                      }
417                  }                  } // diagonals
418              }              } // row chunks
419          }          } // A.symmetric
420      }      } // block size
421  }  }
422    
423  template <typename Matrix,  template <typename Matrix,

Legend:
Removed from v.5243  
changed lines
  Added in v.5244

  ViewVC Help
Powered by ViewVC 1.1.26