Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaSpmvBLC.c
Go to the documentation of this file.
1
14#include <time.h>
15
16#include "fasp.h"
17#include "fasp_block.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Public Functions --*/
22/*---------------------------------*/
23
38void fasp_blas_dblc_aAxpy (const REAL alpha,
39 const dBLCmat *A,
40 const REAL *x,
41 REAL *y)
42{
43 // information of A
44 const INT brow = A->brow;
45
46 // local variables
47 register dCSRmat *A11, *A12, *A21, *A22;
48 register dCSRmat *A13, *A23, *A31, *A32, *A33;
49
50 INT row1, col1;
51 INT row2, col2;
52
53 const register REAL *x1, *x2, *x3;
54 register REAL *y1, *y2, *y3;
55
56 INT i,j;
57 INT start_row, start_col;
58
59 switch (brow) {
60
61 case 2:
62 A11 = A->blocks[0];
63 A12 = A->blocks[1];
64 A21 = A->blocks[2];
65 A22 = A->blocks[3];
66
67 row1 = A11->row;
68 col1 = A11->col;
69
70 x1 = x;
71 x2 = &(x[col1]);
72 y1 = y;
73 y2 = &(y[row1]);
74
75 // y1 = alpha*A11*x1 + alpha*A12*x2 + y1
76 if (A11) fasp_blas_dcsr_aAxpy(alpha, A11, x1, y1);
77 if (A12) fasp_blas_dcsr_aAxpy(alpha, A12, x2, y1);
78
79 // y2 = alpha*A21*x1 + alpha*A22*x2 + y2
80 if (A21) fasp_blas_dcsr_aAxpy(alpha, A21, x1, y2);
81 if (A22) fasp_blas_dcsr_aAxpy(alpha, A22, x2, y2);
82
83 break;
84
85 case 3:
86 A11 = A->blocks[0];
87 A12 = A->blocks[1];
88 A13 = A->blocks[2];
89 A21 = A->blocks[3];
90 A22 = A->blocks[4];
91 A23 = A->blocks[5];
92 A31 = A->blocks[6];
93 A32 = A->blocks[7];
94 A33 = A->blocks[8];
95
96 row1 = A11->row;
97 col1 = A11->col;
98 row2 = A22->row;
99 col2 = A22->col;
100
101 x1 = x;
102 x2 = &(x[col1]);
103 x3 = &(x[col1+col2]);
104 y1 = y;
105 y2 = &(y[row1]);
106 y3 = &(y[row1+row2]);
107
108 // y1 = alpha*A11*x1 + alpha*A12*x2 + alpha*A13*x3 + y1
109 if (A11) fasp_blas_dcsr_aAxpy(alpha, A11, x1, y1);
110 if (A12) fasp_blas_dcsr_aAxpy(alpha, A12, x2, y1);
111 if (A13) fasp_blas_dcsr_aAxpy(alpha, A13, x3, y1);
112
113 // y2 = alpha*A21*x1 + alpha*A22*x2 + alpha*A23*x3 + y2
114 if (A21) fasp_blas_dcsr_aAxpy(alpha, A21, x1, y2);
115 if (A22) fasp_blas_dcsr_aAxpy(alpha, A22, x2, y2);
116 if (A23) fasp_blas_dcsr_aAxpy(alpha, A23, x3, y2);
117
118 // y3 = alpha*A31*x1 + alpha*A32*x2 + alpha*A33*x3 + y2
119 if (A31) fasp_blas_dcsr_aAxpy(alpha, A31, x1, y3);
120 if (A32) fasp_blas_dcsr_aAxpy(alpha, A32, x2, y3);
121 if (A33) fasp_blas_dcsr_aAxpy(alpha, A33, x3, y3);
122
123 break;
124
125 default:
126
127 start_row = 0;
128 start_col = 0;
129
130 for (i=0; i<brow; i++) {
131
132 for (j=0; j<brow; j++) {
133
134 if (A->blocks[i*brow+j]) {
135 fasp_blas_dcsr_aAxpy(alpha, A->blocks[i*brow+j],
136 &(x[start_col]), &(y[start_row]));
137 }
138 start_col = start_col + A->blocks[j*brow+j]->col;
139
140 }
141
142 start_row = start_row + A->blocks[i*brow+i]->row;
143 start_col = 0;
144 }
145
146 break;
147
148 } // end of switch
149
150}
151
164void fasp_blas_dblc_mxv (const dBLCmat *A,
165 const REAL *x,
166 REAL *y)
167{
168 // information of A
169 const INT brow = A->brow;
170
171 // local variables
172 register dCSRmat *A11, *A12, *A21, *A22;
173 register dCSRmat *A13, *A23, *A31, *A32, *A33;
174
175 INT row1, col1;
176 INT row2, col2;
177
178 const register REAL *x1, *x2, *x3;
179 register REAL *y1, *y2, *y3;
180
181 INT i,j;
182 INT start_row, start_col;
183
184 switch (brow) {
185
186 case 2:
187 A11 = A->blocks[0];
188 A12 = A->blocks[1];
189 A21 = A->blocks[2];
190 A22 = A->blocks[3];
191
192 row1 = A11->row;
193 col1 = A11->col;
194
195 x1 = x;
196 x2 = &(x[col1]);
197 y1 = y;
198 y2 = &(y[row1]);
199
200 // y1 = A11*x1 + A12*x2
201 if (A11) fasp_blas_dcsr_mxv(A11, x1, y1);
202 if (A12) fasp_blas_dcsr_aAxpy(1.0, A12, x2, y1);
203
204 // y2 = A21*x1 + A22*x2
205 if (A21) fasp_blas_dcsr_mxv(A21, x1, y2);
206 if (A22) fasp_blas_dcsr_aAxpy(1.0, A22, x2, y2);
207
208 break;
209
210 case 3:
211 A11 = A->blocks[0];
212 A12 = A->blocks[1];
213 A13 = A->blocks[2];
214 A21 = A->blocks[3];
215 A22 = A->blocks[4];
216 A23 = A->blocks[5];
217 A31 = A->blocks[6];
218 A32 = A->blocks[7];
219 A33 = A->blocks[8];
220
221 row1 = A11->row;
222 col1 = A11->col;
223 row2 = A22->row;
224 col2 = A22->col;
225
226 x1 = x;
227 x2 = &(x[col1]);
228 x3 = &(x[col1+col2]);
229 y1 = y;
230 y2 = &(y[row1]);
231 y3 = &(y[row1+row2]);
232
233 // y1 = A11*x1 + A12*x2 + A13*x3 + y1
234 if (A11) fasp_blas_dcsr_mxv(A11, x1, y1);
235 if (A12) fasp_blas_dcsr_aAxpy(1.0, A12, x2, y1);
236 if (A13) fasp_blas_dcsr_aAxpy(1.0, A13, x3, y1);
237
238 // y2 = A21*x1 + A22*x2 + A23*x3 + y2
239 if (A21) fasp_blas_dcsr_mxv(A21, x1, y2);
240 if (A22) fasp_blas_dcsr_aAxpy(1.0, A22, x2, y2);
241 if (A23) fasp_blas_dcsr_aAxpy(1.0, A23, x3, y2);
242
243 // y3 = A31*x1 + A32*x2 + A33*x3 + y2
244 if (A31) fasp_blas_dcsr_mxv(A31, x1, y3);
245 if (A32) fasp_blas_dcsr_aAxpy(1.0, A32, x2, y3);
246 if (A33) fasp_blas_dcsr_aAxpy(1.0, A33, x3, y3);
247
248 break;
249
250 default:
251
252 start_row = 0;
253 start_col = 0;
254
255 for (i=0; i<brow; i++) {
256
257 for (j=0; j<brow; j++){
258
259 if (j==0) {
260 if (A->blocks[i*brow+j]){
261 fasp_blas_dcsr_mxv(A->blocks[i*brow+j], &(x[start_col]), &(y[start_row]));
262 }
263 }
264 else {
265 if (A->blocks[i*brow+j]){
266 fasp_blas_dcsr_aAxpy(1.0, A->blocks[i*brow+j], &(x[start_col]), &(y[start_row]));
267 }
268 }
269 start_col = start_col + A->blocks[j*brow+j]->col;
270 }
271
272 start_row = start_row + A->blocks[i*brow+i]->row;
273 start_col = 0;
274 }
275
276 break;
277
278 } // end of switch
279
280}
281
296void fasp_blas_ldblc_aAxpy (const REAL alpha,
297 const dBLCmat *A,
298 const LONGREAL *x,
299 REAL *y)
300{
301 // information of A
302 const INT brow = A->brow;
303
304 // local variables
305 register dCSRmat *A11, *A12, *A21, *A22;
306 register dCSRmat *A13, *A23, *A31, *A32, *A33;
307
308 INT row1, col1;
309 INT row2, col2;
310
311 const register LONGREAL *x1, *x2, *x3;
312 register REAL *y1, *y2, *y3;
313
314 INT i,j;
315 INT start_row, start_col;
316
317 switch (brow) {
318
319 case 2:
320 A11 = A->blocks[0];
321 A12 = A->blocks[1];
322 A21 = A->blocks[2];
323 A22 = A->blocks[3];
324
325 row1 = A11->row;
326 col1 = A11->col;
327
328 x1 = x;
329 x2 = &(x[col1]);
330 y1 = y;
331 y2 = &(y[row1]);
332
333 // y1 = alpha*A11*x1 + alpha*A12*x2 + y1
334 if (A11) fasp_blas_ldcsr_aAxpy(alpha, A11, x1, y1);
335 if (A12) fasp_blas_ldcsr_aAxpy(alpha, A12, x2, y1);
336
337 // y2 = alpha*A21*x1 + alpha*A22*x2 + y2
338 if (A21) fasp_blas_ldcsr_aAxpy(alpha, A21, x1, y2);
339 if (A22) fasp_blas_ldcsr_aAxpy(alpha, A22, x2, y2);
340
341 break;
342
343 case 3:
344 A11 = A->blocks[0];
345 A12 = A->blocks[1];
346 A13 = A->blocks[2];
347 A21 = A->blocks[3];
348 A22 = A->blocks[4];
349 A23 = A->blocks[5];
350 A31 = A->blocks[6];
351 A32 = A->blocks[7];
352 A33 = A->blocks[8];
353
354 row1 = A11->row;
355 col1 = A11->col;
356 row2 = A22->row;
357 col2 = A22->col;
358
359 x1 = x;
360 x2 = &(x[col1]);
361 x3 = &(x[col1+col2]);
362 y1 = y;
363 y2 = &(y[row1]);
364 y3 = &(y[row1+row2]);
365
366 // y1 = alpha*A11*x1 + alpha*A12*x2 + alpha*A13*x3 + y1
367 if (A11) fasp_blas_ldcsr_aAxpy(alpha, A11, x1, y1);
368 if (A12) fasp_blas_ldcsr_aAxpy(alpha, A12, x2, y1);
369 if (A13) fasp_blas_ldcsr_aAxpy(alpha, A13, x3, y1);
370
371 // y2 = alpha*A21*x1 + alpha*A22*x2 + alpha*A23*x3 + y2
372 if (A21) fasp_blas_ldcsr_aAxpy(alpha, A21, x1, y2);
373 if (A22) fasp_blas_ldcsr_aAxpy(alpha, A22, x2, y2);
374 if (A23) fasp_blas_ldcsr_aAxpy(alpha, A23, x3, y2);
375
376 // y3 = alpha*A31*x1 + alpha*A32*x2 + alpha*A33*x3 + y2
377 if (A31) fasp_blas_ldcsr_aAxpy(alpha, A31, x1, y3);
378 if (A32) fasp_blas_ldcsr_aAxpy(alpha, A32, x2, y3);
379 if (A33) fasp_blas_ldcsr_aAxpy(alpha, A33, x3, y3);
380
381 break;
382
383 default:
384
385 start_row = 0;
386 start_col = 0;
387
388 for (i=0; i<brow; i++) {
389
390 for (j=0; j<brow; j++) {
391
392 if (A->blocks[i*brow+j]) {
393 fasp_blas_ldcsr_aAxpy(alpha, A->blocks[i*brow+j],
394 &(x[start_col]), &(y[start_row]));
395 }
396 start_col = start_col + A->blocks[j*brow+j]->col;
397
398 }
399
400 start_row = start_row + A->blocks[i*brow+i]->row;
401 start_col = 0;
402 }
403
404 break;
405
406 } // end of switch
407}
408
409/*---------------------------------*/
410/*-- End of File --*/
411/*---------------------------------*/
fasp_blas_dblc_mxv
void fasp_blas_dblc_mxv(const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvBLC.c:164
fasp_blas_dblc_aAxpy
void fasp_blas_dblc_aAxpy(const REAL alpha, const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvBLC.c:38
fasp_blas_ldblc_aAxpy
void fasp_blas_ldblc_aAxpy(const REAL alpha, const dBLCmat *A, const LONGREAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvBLC.c:296
fasp_blas_ldcsr_aAxpy
void fasp_blas_ldcsr_aAxpy(const REAL alpha, const dCSRmat *A, const LONGREAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:609
fasp_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_blas_dcsr_aAxpy
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:494
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
LONGREAL
#define LONGREAL
Definition: fasp.h:76
INT
#define INT
Definition: fasp.h:72
fasp_block.h
Header file for FASP block matrices.
dBLCmat
Block REAL CSR matrix format.
Definition: fasp_block.h:74
dBLCmat::brow
INT brow
row number of blocks in A, m
Definition: fasp_block.h:77
dBLCmat::blocks
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:83
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154