Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
blas_csrl.c
Go to the documentation of this file.
1 
10 #include "fasp.h"
11 #include "fasp_functs.h"
12 
13 /*---------------------------------*/
14 /*-- Public Functions --*/
15 /*---------------------------------*/
16 
29  REAL *x,
30  REAL *y)
31 {
32  INT dif = A -> dif;
33  INT *nz_diff = A -> nz_diff;
34  INT *rowindex = A -> index;
35  INT *rowstart = A -> start;
36  INT *ja = A -> ja;
37  REAL *a = A -> val;
38 
39  INT i;
40  INT row, col=0;
41  INT len, rowlen;
42  INT firstrow, lastrow;
43 
44  REAL val0, val1;
45 
46  for (len = 0; len < dif; len ++) {
47  firstrow = rowstart[len];
48  lastrow = rowstart[len+1] - 1;
49  rowlen = nz_diff[len];
50 
51  if (lastrow > firstrow ) {
52  //----------------------------------------------------------
53  // Fully-unrolled code for special case (i.g.,rowlen = 5)
54  // Note: you can also set other special case
55  //----------------------------------------------------------
56  if (rowlen == 5) {
57  for (row = firstrow; row < lastrow; row += 2) {
58  val0 = a[col]*x[ja[col]];
59  val1 = a[col+5]*x[ja[col+5]];
60  col ++;
61 
62  val0 += a[col]*x[ja[col]];
63  val1 += a[col+5]*x[ja[col+5]];
64  col ++;
65 
66  val0 += a[col]*x[ja[col]];
67  val1 += a[col+5]*x[ja[col+5]];
68  col ++;
69 
70  val0 += a[col]*x[ja[col]];
71  val1 += a[col+5]*x[ja[col+5]];
72  col ++;
73 
74  val0 += a[col]*x[ja[col]];
75  val1 += a[col+5]*x[ja[col+5]];
76  col ++;
77 
78  y[rowindex[row]] = val0;
79  y[rowindex[row+1]] = val1;
80 
81  col += 5;
82  }
83  }
84  else {
85  //------------------------------------------------------------------
86  // Unroll-and-jammed code for handling two rows at a time
87  //------------------------------------------------------------------
88 
89  for (row = firstrow; row < lastrow; row += 2) {
90  val0 = 0.0;
91  val1 = 0.0;
92  for (i = 0; i < rowlen; i ++) {
93  val0 += a[col]*x[ja[col]];
94  val1 += a[col+rowlen]*x[ja[col+rowlen]];
95  col ++;
96  }
97  y[rowindex[row]] = val0;
98  y[rowindex[row+1]] = val1;
99  col += rowlen;
100  }
101  }
102  firstrow = row;
103  }
104 
105  //-----------------------------------------------------------
106  // Handle leftover rows that can't be handled in bundles
107  // in the unroll-and -jammed loop
108  //-----------------------------------------------------------
109 
110  for (row = firstrow; row <= lastrow; row ++) {
111  val0 = 0.0;
112  for (i = 0; i < rowlen; i ++) {
113  val0 += a[col]*x[ja[col]];
114  col ++;
115  }
116  y[rowindex[row]] = val0;
117  }
118 
119  }
120 
121 }
122 
123 /*---------------------------------*/
124 /*-- End of File --*/
125 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
#define INT
Definition: fasp.h:64
Main header file for FASP.
Sparse matrix of REAL type in CSRL format.
Definition: fasp.h:265
void fasp_blas_dcsrl_mxv(dCSRLmat *A, REAL *x, REAL *y)
Compute y = A*x for a sparse matrix in CSRL format.
Definition: blas_csrl.c:28