/[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 5219 by caltinay, Tue Oct 21 06:50:20 2014 UTC revision 5220 by caltinay, Thu Oct 23 03:50:37 2014 UTC
# Line 242  void spmv_cds(const Matrix&  A, Line 242  void spmv_cds(const Matrix&  A,
242      // optimization for special case      // optimization for special case
243      if (block_size == 2) {      if (block_size == 2) {
244          if (A.symmetric) {          if (A.symmetric) {
245                // if there is a main diagonal block, it is the first in offsets
246                // and should be skipped in the first loop below since the main
247                // diagonal is processed in the second loop
248                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              {              {
# Line 250  void spmv_cds(const Matrix&  A, Line 254  void spmv_cds(const Matrix&  A,
254                      ValueType sum1 = initialize(y[row]);                      ValueType sum1 = initialize(y[row]);
255                      ValueType sum2 = initialize(y[row+1]);                      ValueType sum2 = initialize(y[row+1]);
256                      // process subdiagonal blocks                      // process subdiagonal blocks
257                      for(IndexType d = 0; d < num_diagonals; d++)                      for(IndexType d = 0; d < num_diagonals-d0; d++)
258                      {                      {
259                          const IndexType diag = num_diagonals-d-1;                          const IndexType diag = num_diagonals-d-1;
                         // skip main diagonal block which is handled below  
                         if (A.diagonal_offsets[diag] == 0)  
                             continue;  
260                          const IndexType col = row - A.diagonal_offsets[diag]*2;                          const IndexType col = row - A.diagonal_offsets[diag]*2;
261                          if (col >= 0 && col <= num_rows)                          if (col >= 0 && col <= num_rows-2)
262                          {                          {
263                              sum1 = reduce(sum1,combine(A.values(col,  2*diag),  x[col]));                              sum1 = reduce(sum1,combine(A.values(col,  2*diag),  x[col]));
264                              sum2 = reduce(sum2,combine(A.values(col+1,2*diag),  x[col]));                              sum2 = reduce(sum2,combine(A.values(col,  2*diag+1),x[col]));
265                              sum1 = reduce(sum1,combine(A.values(col,  2*diag+1),x[col+1]));                              sum1 = reduce(sum1,combine(A.values(col+1,2*diag),  x[col+1]));
266                              sum2 = reduce(sum2,combine(A.values(col+1,2*diag+1),x[col+1]));                              sum2 = reduce(sum2,combine(A.values(col+1,2*diag+1),x[col+1]));
267                          }                          }
268                      }                      }
# Line 269  void spmv_cds(const Matrix&  A, Line 270  void spmv_cds(const Matrix&  A,
270                      for(IndexType d = 0; d < num_diagonals; d++)                      for(IndexType d = 0; d < num_diagonals; d++)
271                      {                      {
272                          const IndexType col = row + A.diagonal_offsets[d]*2;                          const IndexType col = row + A.diagonal_offsets[d]*2;
273                          if (col >= 0 && col <= num_rows)                          if (col >= 0 && col <= num_rows-2)
274                          {                          {
275                              sum1 = reduce(sum1,combine(A.values(row,  2*d),  x[col]));                              sum1 = reduce(sum1,combine(A.values(row,  2*d),  x[col]));
276                              sum2 = reduce(sum2,combine(A.values(row+1,2*d),  x[col]));                              sum2 = reduce(sum2,combine(A.values(row+1,2*d),  x[col]));
# Line 281  void spmv_cds(const Matrix&  A, Line 282  void spmv_cds(const Matrix&  A,
282                      y[row+1] = sum2;                      y[row+1] = sum2;
283                  }                  }
284              }              }
285          } else {          } else { // A.symmetric
286  #pragma omp parallel for  #pragma omp parallel for
287              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)
288              {              {
# Line 293  void spmv_cds(const Matrix&  A, Line 294  void spmv_cds(const Matrix&  A,
294                      for(IndexType d = 0; d < num_diagonals; d++)                      for(IndexType d = 0; d < num_diagonals; d++)
295                      {                      {
296                          const IndexType col = row + A.diagonal_offsets[d]*2;                          const IndexType col = row + A.diagonal_offsets[d]*2;
297                          if (col >= 0 && col <= num_rows)                          if (col >= 0 && col <= num_rows-2)
298                          {                          {
299                              sum1 = reduce(sum1,combine(A.values(row,  2*d),  x[col]));                              sum1 = reduce(sum1,combine(A.values(row,  2*d),  x[col]));
300                              sum2 = reduce(sum2,combine(A.values(row+1,2*d),  x[col]));                              sum2 = reduce(sum2,combine(A.values(row+1,2*d),  x[col]));
# Line 308  void spmv_cds(const Matrix&  A, Line 309  void spmv_cds(const Matrix&  A,
309          } // A.symmetric          } // A.symmetric
310      } else { // block size      } else { // block size
311          if (A.symmetric) {          if (A.symmetric) {
312                // if there is a main diagonal block, it is the first in offsets
313                // and should be skipped in the first loop below since the main
314                // diagonal is processed in the second loop
315                const IndexType d0 = (A.diagonal_offsets[0] == 0 ? 1 : 0);
316  #pragma omp parallel for  #pragma omp parallel for
317              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)              for(IndexType ch = 0; ch < num_rows; ch+=chunksize)
318              {              {
# Line 316  void spmv_cds(const Matrix&  A, Line 321  void spmv_cds(const Matrix&  A,
321                      y[row] = initialize(y[row]);                      y[row] = initialize(y[row]);
322    
323                      // process subdiagonal blocks                      // process subdiagonal blocks
324                      for (IndexType d = 0; d < num_diagonals; d++)                      for (IndexType d = 0; d < num_diagonals-d0; d++)
325                      {                      {
326                          const IndexType diag = num_diagonals-d-1;                          const IndexType diag = num_diagonals-d-1;
                         // skip main diagonal block which is handled below  
                         if (A.diagonal_offsets[diag] == 0)  
                             continue;  
327                          const IndexType col = block_size*(row/block_size) - A.diagonal_offsets[diag]*block_size;                          const IndexType col = block_size*(row/block_size) - A.diagonal_offsets[diag]*block_size;
328                          if (col >= 0 && col <= num_rows-block_size)                          if (col >= 0 && col <= num_rows-block_size)
329                          {                          {
# Line 344  void spmv_cds(const Matrix&  A, Line 346  void spmv_cds(const Matrix&  A,
346                              // for each column in block                              // for each column in block
347                              for (IndexType i = 0; i < block_size; i++)                              for (IndexType i = 0; i < block_size; i++)
348                              {                              {
349                                  const ValueType& Aij = A.values(row+i, d*block_size+row%block_size);                                  const ValueType& Aij = A.values(row, d*block_size+i);
350                                  const ValueType& xj = x[col + i];                                  const ValueType& xj = x[col + i];
351    
352                                  y[row] = reduce(y[row], combine(Aij, xj));                                  y[row] = reduce(y[row], combine(Aij, xj));

Legend:
Removed from v.5219  
changed lines
  Added in v.5220

  ViewVC Help
Powered by ViewVC 1.1.26