Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
blas_bcsr.c
Go to the documentation of this file.
1 
6 #include <time.h>
7 
8 #include "fasp.h"
9 #include "fasp_block.h"
10 #include "fasp_functs.h"
11 
12 /*---------------------------------*/
13 /*-- Public Functions --*/
14 /*---------------------------------*/
15 
30 void fasp_blas_bdcsr_aAxpy (const REAL alpha,
31  block_dCSRmat *A,
32  REAL *x,
33  REAL *y)
34 {
35  // information of A
36  INT brow = A->brow;
37 
38  // local variables
39  register dCSRmat *A11, *A12, *A21, *A22;
40  register dCSRmat *A13, *A23, *A31, *A32, *A33;
41 
42  unsigned INT row1, col1;
43  unsigned INT row2, col2;
44 
45  register REAL *x1, *x2, *y1, *y2;
46  register REAL *x3, *y3;
47 
48  INT i,j;
49  INT start_row;
50  INT start_col;
51 
52  switch (brow) {
53 
54  case 2:
55  A11 = A->blocks[0];
56  A12 = A->blocks[1];
57  A21 = A->blocks[2];
58  A22 = A->blocks[3];
59 
60  row1 = A11->row;
61  col1 = A11->col;
62 
63  x1 = x;
64  x2 = &(x[col1]);
65  y1 = y;
66  y2 = &(y[row1]);
67 
68  // y1 = alpha*A11*x1 + alpha*A12*x2 + y1
69  if (A11) fasp_blas_dcsr_aAxpy(alpha, A11, x1, y1);
70  if (A12) fasp_blas_dcsr_aAxpy(alpha, A12, x2, y1);
71 
72  // y2 = alpha*A21*x1 + alpha*A22*x2 + y2
73  if (A21) fasp_blas_dcsr_aAxpy(alpha, A21, x1, y2);
74  if (A22) fasp_blas_dcsr_aAxpy(alpha, A22, x2, y2);
75 
76  break;
77 
78  case 3:
79  A11 = A->blocks[0];
80  A12 = A->blocks[1];
81  A13 = A->blocks[2];
82  A21 = A->blocks[3];
83  A22 = A->blocks[4];
84  A23 = A->blocks[5];
85  A31 = A->blocks[6];
86  A32 = A->blocks[7];
87  A33 = A->blocks[8];
88 
89  row1 = A11->row;
90  col1 = A11->col;
91  row2 = A22->row;
92  col2 = A22->col;
93 
94  x1 = x;
95  x2 = &(x[col1]);
96  x3 = &(x[col1+col2]);
97  y1 = y;
98  y2 = &(y[row1]);
99  y3 = &(y[row1+row2]);
100 
101  // y1 = alpha*A11*x1 + alpha*A12*x2 + alpha*A13*x3 + y1
102  if (A11) fasp_blas_dcsr_aAxpy(alpha, A11, x1, y1);
103  if (A12) fasp_blas_dcsr_aAxpy(alpha, A12, x2, y1);
104  if (A13) fasp_blas_dcsr_aAxpy(alpha, A13, x3, y1);
105 
106  // y2 = alpha*A21*x1 + alpha*A22*x2 + alpha*A23*x3 + y2
107  if (A21) fasp_blas_dcsr_aAxpy(alpha, A21, x1, y2);
108  if (A22) fasp_blas_dcsr_aAxpy(alpha, A22, x2, y2);
109  if (A23) fasp_blas_dcsr_aAxpy(alpha, A23, x3, y2);
110 
111  // y3 = alpha*A31*x1 + alpha*A32*x2 + alpha*A33*x3 + y2
112  if (A31) fasp_blas_dcsr_aAxpy(alpha, A31, x1, y3);
113  if (A32) fasp_blas_dcsr_aAxpy(alpha, A32, x2, y3);
114  if (A33) fasp_blas_dcsr_aAxpy(alpha, A33, x3, y3);
115 
116  break;
117 
118  default:
119 
120  start_row = 0;
121  start_col = 0;
122 
123  for (i=0; i<brow; i++) {
124 
125  for (j=0; j<brow; j++){
126 
127  if (A->blocks[i*brow+j]){
128  fasp_blas_dcsr_aAxpy(alpha, A->blocks[i*brow+j], &(x[start_col]), &(y[start_row]));
129  }
130  start_col = start_col + A->blocks[j*brow+j]->col;
131  }
132 
133  start_row = start_row + A->blocks[i*brow+i]->row;
134  start_col = 0;
135  }
136 
137  break;
138 
139  } // end of switch
140 
141 }
142 
156  REAL *x,
157  REAL *y)
158 {
159  // information of A
160  INT brow = A->brow;
161 
162  // local variables
163  register dCSRmat *A11, *A12, *A21, *A22;
164  register dCSRmat *A13, *A23, *A31, *A32, *A33;
165 
166  unsigned INT row1, col1;
167  unsigned INT row2, col2;
168 
169  register REAL *x1, *x2, *y1, *y2;
170  register REAL *x3, *y3;
171 
172  INT i,j;
173  INT start_row;
174  INT start_col;
175 
176  switch (brow) {
177 
178  case 2:
179  A11 = A->blocks[0];
180  A12 = A->blocks[1];
181  A21 = A->blocks[2];
182  A22 = A->blocks[3];
183 
184  row1 = A11->row;
185  col1 = A11->col;
186 
187  x1 = x;
188  x2 = &(x[col1]);
189  y1 = y;
190  y2 = &(y[row1]);
191 
192  // y1 = A11*x1 + A12*x2
193  if (A11) fasp_blas_dcsr_mxv(A11, x1, y1);
194  if (A12) fasp_blas_dcsr_aAxpy(1.0, A12, x2, y1);
195 
196  // y2 = A21*x1 + A22*x2
197  if (A21) fasp_blas_dcsr_mxv(A21, x1, y2);
198  if (A22) fasp_blas_dcsr_aAxpy(1.0, A22, x2, y2);
199 
200  break;
201 
202  case 3:
203  A11 = A->blocks[0];
204  A12 = A->blocks[1];
205  A13 = A->blocks[2];
206  A21 = A->blocks[3];
207  A22 = A->blocks[4];
208  A23 = A->blocks[5];
209  A31 = A->blocks[6];
210  A32 = A->blocks[7];
211  A33 = A->blocks[8];
212 
213  row1 = A11->row;
214  col1 = A11->col;
215  row2 = A22->row;
216  col2 = A22->col;
217 
218  x1 = x;
219  x2 = &(x[col1]);
220  x3 = &(x[col1+col2]);
221  y1 = y;
222  y2 = &(y[row1]);
223  y3 = &(y[row1+row2]);
224 
225  // y1 = A11*x1 + A12*x2 + A13*x3 + y1
226  if (A11) fasp_blas_dcsr_mxv(A11, x1, y1);
227  if (A12) fasp_blas_dcsr_aAxpy(1.0, A12, x2, y1);
228  if (A13) fasp_blas_dcsr_aAxpy(1.0, A13, x3, y1);
229 
230  // y2 = A21*x1 + A22*x2 + A23*x3 + y2
231  if (A21) fasp_blas_dcsr_mxv(A21, x1, y2);
232  if (A22) fasp_blas_dcsr_aAxpy(1.0, A22, x2, y2);
233  if (A23) fasp_blas_dcsr_aAxpy(1.0, A23, x3, y2);
234 
235  // y3 = A31*x1 + A32*x2 + A33*x3 + y2
236  if (A31) fasp_blas_dcsr_mxv(A31, x1, y3);
237  if (A32) fasp_blas_dcsr_aAxpy(1.0, A32, x2, y3);
238  if (A33) fasp_blas_dcsr_aAxpy(1.0, A33, x3, y3);
239 
240  break;
241 
242  default:
243 
244  start_row = 0;
245  start_col = 0;
246 
247  for (i=0; i<brow; i++) {
248 
249  for (j=0; j<brow; j++){
250 
251  if (j==0) {
252  if (A->blocks[i*brow+j]){
253  fasp_blas_dcsr_mxv(A->blocks[i*brow+j], &(x[start_col]), &(y[start_row]));
254  }
255  }
256  else {
257  if (A->blocks[i*brow+j]){
258  fasp_blas_dcsr_aAxpy(1.0, A->blocks[i*brow+j], &(x[start_col]), &(y[start_row]));
259  }
260  }
261  start_col = start_col + A->blocks[j*brow+j]->col;
262  }
263 
264  start_row = start_row + A->blocks[i*brow+i]->row;
265  start_col = 0;
266  }
267 
268  break;
269 
270  } // end of switch
271 
272 }
273 
288 void fasp_blas_bdbsr_aAxpy (const REAL alpha,
289  block_BSR *A,
290  REAL *x,
291  REAL *y)
292 {
293  register dBSRmat *Arr = &(A->ResRes);
294  register dCSRmat *Arw = &(A->ResWel);
295  register dCSRmat *Awr = &(A->WelRes);
296  register dCSRmat *Aww = &(A->WelWel);
297 
298  unsigned INT Nr = Arw->row;
299 
300  register REAL *xr = x;
301  register REAL *xw = &(x[Nr]);
302  register REAL *yr = y;
303  register REAL *yw = &(y[Nr]);
304 
305  // yr = alpha*Arr*xr + alpha*Arw*xw + yr
306  fasp_blas_dbsr_aAxpy(alpha, Arr, xr, yr);
307  fasp_blas_dcsr_aAxpy(alpha, Arw, xw, yr);
308 
309  // yw = alpha*Awr*xr + alpha*Aww*xw + yw
310  fasp_blas_dcsr_aAxpy(alpha, Awr, xr, yw);
311  fasp_blas_dcsr_aAxpy(alpha, Aww, xw, yw);
312 }
313 
327  REAL *x,
328  REAL *y)
329 {
330  register dBSRmat *Arr = &(A->ResRes);
331  register dCSRmat *Arw = &(A->ResWel);
332  register dCSRmat *Awr = &(A->WelRes);
333  register dCSRmat *Aww = &(A->WelWel);
334 
335  unsigned INT Nr = Arw->row;
336 
337  register REAL *xr = x;
338  register REAL *xw = &(x[Nr]);
339  register REAL *yr = y;
340  register REAL *yw = &(y[Nr]);
341 
342  // yr = alpha*Arr*xr + alpha*Arw*xw + yr
343  fasp_blas_dbsr_mxv(Arr, xr, yr);
344  fasp_blas_dcsr_aAxpy(1.0, Arw, xw, yr);
345 
346  // yw = alpha*Awr*xr + alpha*Aww*xw + yw
347  fasp_blas_dcsr_mxv(Awr, xr, yw);
348  fasp_blas_dcsr_aAxpy(1.0, Aww, xw, yw);
349 }
350 
351 /*---------------------------------*/
352 /*-- End of File --*/
353 /*---------------------------------*/
dCSRmat WelWel
well-well block
Definition: fasp_block.h:184
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:93
dBSRmat ResRes
reservoir-reservoir block
Definition: fasp_block.h:175
#define REAL
Definition: fasp.h:67
void fasp_blas_bdcsr_mxv(block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:155
dCSRmat WelRes
well-reservoir block
Definition: fasp_block.h:181
Block REAL matrix format for reservoir simulation.
Definition: fasp_block.h:172
void fasp_blas_bdcsr_aAxpy(const REAL alpha, block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_bcsr.c:30
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
#define INT
Definition: fasp.h:64
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
void fasp_blas_bdbsr_aAxpy(const REAL alpha, block_BSR *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_bcsr.c:288
void fasp_blas_bdbsr_mxv(block_BSR *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:326
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT brow
row number of blocks in A, m
Definition: fasp_block.h:87
Block REAL CSR matrix format.
Definition: fasp_block.h:84
void fasp_blas_dbsr_mxv(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x.
Definition: blas_bsr.c:895
void fasp_blas_dcsr_aAxpy(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_csr.c:479
Header file for FASP block matrices.
dCSRmat ResWel
reservoir-well block
Definition: fasp_block.h:178
void fasp_blas_dbsr_aAxpy(const REAL alpha, dBSRmat *A, REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: blas_bsr.c:337
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225