Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreSTR.c
Go to the documentation of this file.
1
15#include <math.h>
16
17#include "fasp.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Declare Private Functions --*/
22/*---------------------------------*/
23
24static inline void fasp_darray_cp_nc3 (const REAL *x, REAL *y);
25static inline void fasp_darray_cp_nc5 (const REAL *x, REAL *y);
26static inline void fasp_darray_cp_nc7 (const REAL *x, REAL *y);
27
28/*---------------------------------*/
29/*-- Public Functions --*/
30/*---------------------------------*/
31
44void fasp_precond_dstr_diag (REAL *r,
45 REAL *z,
46 void *data)
47{
48 const precond_diag_str *diag = (precond_diag_str *)data;
49 const REAL *diagptr = diag->diag.val;
50 const INT nc = diag->nc, nc2 = nc*nc;
51 const INT m = diag->diag.row/nc2;
52
53 INT i;
54 for ( i=0; i<m; ++i ) {
55 fasp_blas_smat_mxv(&(diagptr[i*nc2]),&(r[i*nc]),&(z[i*nc]),nc);
56 }
57}
58
71void fasp_precond_dstr_ilu0 (REAL *r,
72 REAL *z,
73 void *data)
74{
75 INT i, ic, ic2;
76 REAL *zz,*zr,*tc;
77 INT nline, nplane;
78
79 dSTRmat *ILU_data=(dSTRmat *)data;
80 INT m=ILU_data->ngrid;
81 INT nc=ILU_data->nc;
82 INT nc2=nc*nc;
83 INT nx=ILU_data->nx;
84 INT ny=ILU_data->ny;
85 INT nz=ILU_data->nz;
86 INT nxy=ILU_data->nxy;
87 INT size=m*nc;
88
89#if DEBUG_MODE > 0
90 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
91#endif
92
93 if (nx == 1) {
94 nline = ny;
95 nplane = m;
96 }
97 else if (ny == 1) {
98 nline = nx;
99 nplane = m;
100 }
101 else if (nz == 1) {
102 nline = nx;
103 nplane = m;
104 }
105 else {
106 nline = nx;
107 nplane = nxy;
108 }
109
110 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
111
112 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
113
114 zr=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
115
116 // copy residual r to zr, to save r
117 memcpy(zr,r,(size)*sizeof(REAL));
118
119 if (nc == 1) {
120 // forward sweep: solve unit lower matrix equation L*zz=zr
121 zz[0]=zr[0];
122 for (i=1;i<m;++i) {
123 zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
124 if (i>=nline) zz[i]=zz[i]-ILU_data->offdiag[2][i-nline]*zz[i-nline];
125 if (i>=nplane) zz[i]=zz[i]-ILU_data->offdiag[4][i-nplane]*zz[i-nplane];
126 }
127
128 // backward sweep: solve upper matrix equation U*z=zz
129 z[m-1]=zz[m-1]*ILU_data->diag[m-1];
130 for (i=m-2;i>=0;i--) {
131 zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
132 if (i<m-nline) zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nline];
133 if (i<m-nplane) zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nplane];
134 z[i]=zz[i]*ILU_data->diag[i];
135 }
136
137 } // end if (nc == 1)
138
139 else if (nc == 3) {
140 // forward sweep: solve unit lower matrix equation L*zz=zr
141 fasp_darray_cp_nc3(&(zr[0]),&(zz[0]));
142
143 for (i=1;i<m;++i) {
144 ic=i*nc;
145
146 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
147 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
148 if (i>=nline) {
149 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
150 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
151 }
152 if (i>=nplane) {
153 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
154 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
155 }
156 fasp_darray_cp_nc3(&(zr[ic]),&(zz[ic]));
157 } // end for (i=1;i<m;++i)
158
159 // backward sweep: solve upper matrix equation U*z=zz
160 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
161
162 for (i=m-2;i>=0;i--) {
163
164 ic=i*nc;
165 ic2=i*nc2;
166
167 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
168 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
169
170 if (i<m-nline) {
171 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
172 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
173 }
174
175 if (i<m-nplane) {
176 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
177 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
178 }
179
180 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
181 } // end for for (i=m-2;i>=0;i--)
182
183 } // end else if (nc == 3)
184
185 else if (nc == 5) {
186 // forward sweep: solve unit lower matrix equation L*zz=zr
187 fasp_darray_cp_nc5(&(zr[0]),&(zz[0]));
188
189 for (i=1;i<m;++i) {
190 ic=i*nc;
191
192 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
193 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
194 if (i>=nline) {
195 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
196 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
197 }
198 if (i>=nplane) {
199 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
200 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
201 }
202 fasp_darray_cp_nc5(&(zr[ic]),&(zz[ic]));
203 } // end for (i=1;i<m;++i)
204
205 // backward sweep: solve upper matrix equation U*z=zz
206 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
207
208 for (i=m-2;i>=0;i--) {
209
210 ic=i*nc;
211 ic2=i*nc2;
212
213 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
214 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
215
216 if (i<m-nline) {
217 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
218 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
219 }
220
221 if (i<m-nplane) {
222 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
223 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
224 }
225
226 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
227 } // end for for (i=m-2;i>=0;i--)
228
229 } // end else if (nc == 5)
230
231
232 else if (nc == 7) {
233 // forward sweep: solve unit lower matrix equation L*zz=zr
234 fasp_darray_cp_nc7(&(zr[0]),&(zz[0]));
235
236 for (i=1;i<m;++i) {
237 ic=i*nc;
238
239 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
240 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
241 if (i>=nline) {
242 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
243 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
244 }
245 if (i>=nplane) {
246 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
247 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
248 }
249 fasp_darray_cp_nc7(&(zr[ic]),&(zz[ic]));
250 } // end for (i=1;i<m;++i)
251
252 // backward sweep: solve upper matrix equation U*z=zz
253 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
254
255 for (i=m-2;i>=0;i--) {
256
257 ic=i*nc;
258 ic2=i*nc2;
259
260 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
261 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
262
263 if (i<m-nline) {
264 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
265 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
266 }
267
268 if (i<m-nplane) {
269 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
270 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
271 }
272
273 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
274 } // end for for (i=m-2;i>=0;i--)
275
276 } // end else if (nc == 7)
277
278 else {
279 // forward sweep: solve unit lower matrix equation L*zz=zr
280 fasp_darray_cp(nc,&(zr[0]),&(zz[0]));
281 for (i=1;i<m;++i) {
282 ic=i*nc;
283
284 fasp_blas_smat_mxv(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc,nc);
285 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
286
287 if (i>=nline) {
288 fasp_blas_smat_mxv(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc,nc);
289 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
290 }
291
292 if (i>=nplane) {
293 fasp_blas_smat_mxv(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc,nc);
294 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
295 }
296
297 fasp_darray_cp(nc,&(zr[ic]),&(zz[ic]));
298
299 } // end for (i=1; i<m; ++i)
300
301 // backward sweep: solve upper matrix equation U*z=zz
302 fasp_blas_smat_mxv(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]),nc);
303
304 for (i=m-2;i>=0;i--) {
305 ic=i*nc;
306 ic2=i*nc2;
307
308 fasp_blas_smat_mxv(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc,nc);
309 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
310
311 if (i<m-nline) {
312 fasp_blas_smat_mxv(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc,nc);
313 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
314 }
315
316 if (i<m-nplane) {
317 fasp_blas_smat_mxv(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc,nc);
318 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
319 }
320
321 fasp_blas_smat_mxv(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]),nc);
322
323 }// end for (i=m-2;i>=0;i--)
324 } // end else
325
326 fasp_mem_free(zr); zr = NULL;
327 fasp_mem_free(zz); zz = NULL;
328 fasp_mem_free(tc); tc = NULL;
329
330#if DEBUG_MODE > 0
331 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
332#endif
333
334 return;
335}
336
349void fasp_precond_dstr_ilu1 (REAL *r,
350 REAL *z,
351 void *data)
352{
353 REAL *zz,*zr,*tc;
354
355 dSTRmat *ILU_data=(dSTRmat *)data;
356 INT i,ic, ic2;
357 INT m=ILU_data->ngrid;
358 INT nc=ILU_data->nc;
359 INT nc2=nc*nc;
360 INT nx=ILU_data->nx;
361 INT ny=ILU_data->ny;
362 INT nz=ILU_data->nz;
363 INT nxy=ILU_data->nxy;
364 INT size=m*nc;
365 INT nline, nplane;
366
367 if (nx == 1) {
368 nline = ny;
369 nplane = m;
370 }
371 else if (ny == 1) {
372 nline = nx;
373 nplane = m;
374 }
375 else if (nz == 1) {
376 nline = nx;
377 nplane = m;
378 }
379 else {
380 nline = nx;
381 nplane = nxy;
382 }
383
384 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
385
386 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
387
388 zr=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
389
390 // copy residual r to zr, to save r
391 for (i=0;i<size;++i) zr[i]=r[i];
392 if (nc == 1) {
393 // forward sweep: solve unit lower matrix equation L*zz=zr
394 zz[0]=zr[0];
395 for (i=1;i<m;++i) {
396
397 zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
398 if (i>=nline-1)
399 zz[i]=zz[i]-ILU_data->offdiag[2][i-nline+1]*zz[i-nline+1];
400
401 if (i>=nline)
402 zz[i]=zz[i]-ILU_data->offdiag[4][i-nline]*zz[i-nline];
403 if (i>=nplane-nline)
404 zz[i]=zz[i]-ILU_data->offdiag[6][i-nplane+nline]*zz[i-nplane+nline];
405 if (i>=nplane-1)
406 zz[i]=zz[i]-ILU_data->offdiag[8][i-nplane+1]*zz[i-nplane+1];
407 if (i>=nplane)
408 zz[i]=zz[i]-ILU_data->offdiag[10][i-nplane]*zz[i-nplane];
409 }
410
411 // backward sweep: solve upper matrix equation U*z=zz
412
413 z[m-1]=zz[m-1]*ILU_data->diag[m-1];
414 for (i=m-2;i>=0;i--) {
415
416 zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
417 if (i+nline-1<m)
418 zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nline-1];
419 if (i+nline<m)
420 zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nline];
421 if (i+nplane-nline<m)
422 zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nplane-nline];
423 if (i+nplane-1<m)
424 zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nplane-1];
425 if (i+nplane<m)
426 zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nplane];
427
428 z[i]=ILU_data->diag[i]*zz[i];
429
430 }
431
432 } // end if (nc == 1)
433
434 else if (nc == 3) {
435
436 // forward sweep: solve unit lower matrix equation L*zz=zr
437 fasp_darray_cp_nc3(&(zr[0]),&(zz[0]));
438
439 for (i=1;i<m;++i) {
440 ic=i*nc;
441
442 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
443 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
444 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
445
446 if (i>=nline-1) {
447 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
448 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
449 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
450 }
451
452 if (i>=nline) {
453 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
454 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
455 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
456 }
457
458 if (i>=nplane-nline) {
459 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
460 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
461 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
462 }
463
464 if (i>=nplane-1) {
465 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
466 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
467 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
468 }
469
470 if (i>=nplane) {
471 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
472 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
473 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
474 }
475
476 fasp_darray_cp_nc3(&(zr[ic]),&(zz[ic]));
477 }
478
479 // backward sweep: solve upper matrix equation U*z=zz
480
481 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
482 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
483
484 for (i=m-2;i>=0;i--) {
485 ic=i*nc;
486 ic2=ic*nc;
487
488 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
489 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
490 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
491
492 if (i+nline-1<m) {
493 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
494 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
495 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
496 }
497
498 if (i+nline<m) {
499 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
500 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
501 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
502 }
503
504 if (i+nplane-nline<m) {
505 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
506 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
507 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
508 }
509
510 if (i+nplane-1<m) {
511 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
512 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
513 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
514 }
515
516 if (i+nplane<m) {
517 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
518 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
519 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
520 }
521
522 //z[i]=ILU_data->diag[i]*zz[i];
523 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
524 } // end for (i=m-2;i>=0;i--)
525
526 } // end if (nc == 3)
527
528 else if (nc == 5) {
529
530 // forward sweep: solve unit lower matrix equation L*zz=zr
531 fasp_darray_cp_nc5(&(zr[0]),&(zz[0]));
532
533 for (i=1;i<m;++i) {
534 ic=i*nc;
535
536 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
537 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
538 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
539
540 if (i>=nline-1) {
541 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
542 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
543 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
544 }
545
546 if (i>=nline) {
547 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
548 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
549 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
550 }
551
552 if (i>=nplane-nline) {
553 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
554 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
555 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
556 }
557
558 if (i>=nplane-1) {
559 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
560 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
561 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
562 }
563
564 if (i>=nplane) {
565 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
566 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
567 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
568 }
569
570 fasp_darray_cp_nc5(&(zr[ic]),&(zz[ic]));
571 }
572
573 // backward sweep: solve upper matrix equation U*z=zz
574
575 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
576 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
577
578 for (i=m-2;i>=0;i--) {
579 ic=i*nc;
580 ic2=ic*nc;
581
582 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
583 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
584 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
585
586 if (i+nline-1<m) {
587 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
588 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
589 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
590 }
591
592 if (i+nline<m) {
593 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
594 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
595 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
596 }
597
598 if (i+nplane-nline<m) {
599 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
600 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
601 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
602 }
603
604 if (i+nplane-1<m) {
605 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
606 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
607 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
608 }
609
610 if (i+nplane<m) {
611 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
612 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
613 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
614 }
615
616 //z[i]=ILU_data->diag[i]*zz[i];
617 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
618 } // end for (i=m-2;i>=0;i--)
619
620 } // end if (nc == 5)
621
622 else if (nc == 7) {
623
624 // forward sweep: solve unit lower matrix equation L*zz=zr
625 fasp_darray_cp_nc7(&(zr[0]),&(zz[0]));
626
627 for (i=1;i<m;++i) {
628 ic=i*nc;
629
630 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
631 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
632 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
633
634 if (i>=nline-1) {
635 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
636 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
637 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
638 }
639
640 if (i>=nline) {
641 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
642 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
643 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
644 }
645
646 if (i>=nplane-nline) {
647 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
648 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
649 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
650 }
651
652 if (i>=nplane-1) {
653 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
654 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
655 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
656 }
657
658 if (i>=nplane) {
659 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
660 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
661 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
662 }
663
664 fasp_darray_cp_nc7(&(zr[ic]),&(zz[ic]));
665 }
666
667 // backward sweep: solve upper matrix equation U*z=zz
668
669 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
670 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
671
672 for (i=m-2;i>=0;i--) {
673 ic=i*nc;
674 ic2=ic*nc;
675
676 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
677 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
678 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
679
680 if (i+nline-1<m) {
681 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
682 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
683 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
684 }
685
686 if (i+nline<m) {
687 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
688 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
689 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
690 }
691
692 if (i+nplane-nline<m) {
693 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
694 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
695 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
696 }
697
698 if (i+nplane-1<m) {
699 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
700 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
701 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
702 }
703
704 if (i+nplane<m) {
705 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
706 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
707 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
708 }
709
710 //z[i]=ILU_data->diag[i]*zz[i];
711 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
712 } // end for (i=m-2;i>=0;i--)
713
714 } // end if (nc == 7)
715
716 else {
717 // forward sweep: solve unit lower matrix equation L*zz=zr
718 fasp_darray_cp(nc,&(zr[0]),&(zz[0]));
719
720 for (i=1;i<m;++i) {
721 ic=i*nc;
722 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
723 fasp_blas_smat_mxv(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc,nc);
724 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
725
726 if (i>=nline-1) {
727 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
728 fasp_blas_smat_mxv(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc,nc);
729 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
730 }
731
732 if (i>=nline) {
733 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
734 fasp_blas_smat_mxv(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc,nc);
735 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
736 }
737
738 if (i>=nplane-nline) {
739 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
740 fasp_blas_smat_mxv(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc,nc);
741 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
742 }
743
744 if (i>=nplane-1) {
745 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
746 fasp_blas_smat_mxv(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc,nc);
747 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
748 }
749
750 if (i>=nplane) {
751 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
752 fasp_blas_smat_mxv(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc,nc);
753 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
754 }
755 fasp_darray_cp(nc,&(zr[ic]),&(zz[ic]));
756 }
757
758 // backward sweep: solve upper matrix equation U*z=zz
759
760 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
761 fasp_blas_smat_mxv(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]),nc);
762
763 for (i=m-2;i>=0;i--) {
764 ic=i*nc;
765 ic2=ic*nc;
766 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
767 fasp_blas_smat_mxv(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc,nc);
768 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
769
770 if (i+nline-1<m) {
771 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
772 fasp_blas_smat_mxv(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc,nc);
773 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
774 }
775
776 if (i+nline<m) {
777 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
778 fasp_blas_smat_mxv(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc,nc);
779 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
780 }
781
782 if (i+nplane-nline<m) {
783 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
784 fasp_blas_smat_mxv(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc,nc);
785 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
786 }
787
788 if (i+nplane-1<m) {
789 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
790 fasp_blas_smat_mxv(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc,nc);
791 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
792 }
793
794 if (i+nplane<m) {
795 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
796 fasp_blas_smat_mxv(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc,nc);
797 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
798 }
799
800 //z[i]=ILU_data->diag[i]*zz[i];
801 fasp_blas_smat_mxv(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]),nc);
802 }
803 } // end else
804
805 fasp_mem_free(zr); zr = NULL;
806 fasp_mem_free(zz); zz = NULL;
807 fasp_mem_free(tc); tc = NULL;
808
809 return;
810}
811
824void fasp_precond_dstr_ilu0_forward (REAL *r,
825 REAL *z,
826 void *data)
827{
828 INT i, ic;
829 REAL *zz,*zr,*tc;
830 INT nline, nplane;
831
832 dSTRmat *ILU_data=(dSTRmat *)data;
833 INT m=ILU_data->ngrid;
834 INT nc=ILU_data->nc;
835 INT nc2=nc*nc;
836 INT nx=ILU_data->nx;
837 INT ny=ILU_data->ny;
838 INT nz=ILU_data->nz;
839 INT nxy=ILU_data->nxy;
840 INT size=m*nc;
841
842 if (nx == 1) {
843 nline = ny;
844 nplane = m;
845 }
846 else if (ny == 1) {
847 nline = nx;
848 nplane = m;
849 }
850 else if (nz == 1) {
851 nline = nx;
852 nplane = m;
853 }
854 else {
855 nline = nx;
856 nplane = nxy;
857 }
858
859 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
860
861 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
862
863 zr=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
864
865 // copy residual r to zr, to save r
866 memcpy(zr,r,(size)*sizeof(REAL));
867 if (nc == 1) {
868 // forward sweep: solve unit lower matrix equation L*zz=zr
869 zz[0]=zr[0];
870 for (i=1;i<m;++i) {
871 zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
872 if (i>=nline) zz[i]=zz[i]-ILU_data->offdiag[2][i-nline]*zz[i-nline];
873 if (i>=nplane) zz[i]=zz[i]-ILU_data->offdiag[4][i-nplane]*zz[i-nplane];
874 }
875 } // end if (nc == 1)
876
877 else if (nc == 3) {
878 // forward sweep: solve unit lower matrix equation L*zz=zr
879 fasp_darray_cp_nc3(&(zr[0]),&(zz[0]));
880
881 for (i=1;i<m;++i) {
882 ic=i*nc;
883 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
884 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
885 if (i>=nline) {
886 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
887 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
888 }
889 if (i>=nplane) {
890 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
891 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
892 }
893 fasp_darray_cp_nc3(&(zr[ic]),&(zz[ic]));
894 } // end for (i=1;i<m;++i)
895
896 } // end else if (nc == 3)
897
898 else if (nc == 5) {
899 // forward sweep: solve unit lower matrix equation L*zz=zr
900 fasp_darray_cp_nc5(&(zr[0]),&(zz[0]));
901
902 for (i=1;i<m;++i) {
903 ic=i*nc;
904 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
905 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
906 if (i>=nline) {
907 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
908 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
909 }
910 if (i>=nplane) {
911 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
912 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
913 }
914 fasp_darray_cp_nc5(&(zr[ic]),&(zz[ic]));
915 } // end for (i=1;i<m;++i)
916
917 } // end else if (nc == 5)
918
919
920 else if (nc == 7) {
921 // forward sweep: solve unit lower matrix equation L*zz=zr
922 fasp_darray_cp_nc7(&(zr[0]),&(zz[0]));
923
924 for (i=1;i<m;++i) {
925 ic=i*nc;
926 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
927 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
928 if (i>=nline) {
929 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
930 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
931 }
932 if (i>=nplane) {
933 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
934 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
935 }
936 fasp_darray_cp_nc7(&(zr[ic]),&(zz[ic]));
937 } // end for (i=1;i<m;++i)
938
939 } // end else if (nc == 7)
940
941
942 else {
943 // forward sweep: solve unit lower matrix equation L*zz=zr
944 fasp_darray_cp(nc,&(zr[0]),&(zz[0]));
945 for (i=1;i<m;++i) {
946 ic=i*nc;
947 fasp_blas_smat_mxv(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc,nc);
948 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
949
950 if (i>=nline) {
951 fasp_blas_smat_mxv(&(ILU_data->offdiag[2][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc,nc);
952 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
953 }
954
955 if (i>=nplane) {
956 fasp_blas_smat_mxv(&(ILU_data->offdiag[4][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc,nc);
957 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
958 }
959
960 fasp_darray_cp(nc,&(zr[ic]),&(zz[ic]));
961
962 } // end for (i=1; i<m; ++i)
963
964 } // end else
965
966 memcpy(z,zz,(size)*sizeof(REAL));
967
968 fasp_mem_free(zr); zr = NULL;
969 fasp_mem_free(zz); zz = NULL;
970 fasp_mem_free(tc); tc = NULL;
971
972 return;
973}
974
987void fasp_precond_dstr_ilu0_backward (REAL *r,
988 REAL *z,
989 void *data)
990{
991 INT i, ic, ic2;
992 REAL *zz,*tc;
993 INT nline, nplane;
994
995 dSTRmat *ILU_data=(dSTRmat *)data;
996 INT m=ILU_data->ngrid;
997 INT nc=ILU_data->nc;
998 INT nc2=nc*nc;
999 INT nx=ILU_data->nx;
1000 INT ny=ILU_data->ny;
1001 INT nz=ILU_data->nz;
1002 INT nxy=ILU_data->nxy;
1003 INT size=m*nc;
1004
1005 if (nx == 1) {
1006 nline = ny;
1007 nplane = m;
1008 }
1009 else if (ny == 1) {
1010 nline = nx;
1011 nplane = m;
1012 }
1013 else if (nz == 1) {
1014 nline = nx;
1015 nplane = m;
1016 }
1017 else {
1018 nline = nx;
1019 nplane = nxy;
1020 }
1021
1022 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
1023
1024 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
1025
1026 // copy residual r to zr, to save r
1027 memcpy(zz,r,(size)*sizeof(REAL));
1028 if (nc == 1) {
1029 // backward sweep: solve upper matrix equation U*z=zz
1030
1031 z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1032 for (i=m-2;i>=0;i--) {
1033 zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1034 if (i<m-nline) zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nline];
1035 if (i<m-nplane) zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nplane];
1036 z[i]=zz[i]*ILU_data->diag[i];
1037 }
1038
1039 } // end if (nc == 1)
1040
1041 else if (nc == 3) {
1042 // backward sweep: solve upper matrix equation U*z=zz
1043 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1044
1045 for (i=m-2;i>=0;i--) {
1046
1047 ic=i*nc;
1048 ic2=i*nc2;
1049
1050 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1051 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1052
1053 if (i<m-nline) {
1054 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
1055 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1056 }
1057
1058 if (i<m-nplane) {
1059 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
1060 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1061 }
1062
1063 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1064 } // end for for (i=m-2;i>=0;i--)
1065
1066 } // end else if (nc == 3)
1067
1068 else if (nc == 5) {
1069 // backward sweep: solve upper matrix equation U*z=zz
1070 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1071
1072 for (i=m-2;i>=0;i--) {
1073
1074 ic=i*nc;
1075 ic2=i*nc2;
1076
1077 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1078 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1079
1080 if (i<m-nline) {
1081 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
1082 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1083 }
1084
1085 if (i<m-nplane) {
1086 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
1087 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1088 }
1089
1090 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1091 } // end for for (i=m-2;i>=0;i--)
1092
1093 } // end else if (nc == 5)
1094
1095 else if (nc == 7) {
1096 // backward sweep: solve upper matrix equation U*z=zz
1097 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1098
1099 for (i=m-2;i>=0;i--) {
1100
1101 ic=i*nc;
1102 ic2=i*nc2;
1103
1104 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1105 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1106
1107 if (i<m-nline) {
1108 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc);
1109 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1110 }
1111
1112 if (i<m-nplane) {
1113 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc);
1114 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1115 }
1116
1117 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1118 } // end for for (i=m-2;i>=0;i--)
1119
1120 } // end else if (nc == 7)
1121
1122
1123 else
1124 {
1125 // backward sweep: solve upper matrix equation U*z=zz
1126 fasp_blas_smat_mxv(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]),nc);
1127
1128 for (i=m-2;i>=0;i--) {
1129 ic=i*nc;
1130 ic2=i*nc2;
1131
1132 fasp_blas_smat_mxv(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc,nc);
1133 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1134
1135 if (i<m-nline) {
1136 fasp_blas_smat_mxv(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline)*nc]),tc,nc);
1137 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1138 }
1139
1140 if (i<m-nplane) {
1141 fasp_blas_smat_mxv(&(ILU_data->offdiag[5][ic2]),&(z[(i+nplane)*nc]),tc,nc);
1142 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1143 }
1144
1145 fasp_blas_smat_mxv(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]),nc);
1146
1147 }// end for (i=m-2;i>=0;i--)
1148 } // end else
1149
1150 fasp_mem_free(zz); zz = NULL;
1151 fasp_mem_free(tc); tc = NULL;
1152
1153 return;
1154}
1155
1168void fasp_precond_dstr_ilu1_forward (REAL *r,
1169 REAL *z,
1170 void *data)
1171{
1172 REAL *zz,*zr,*tc;
1173
1174 dSTRmat *ILU_data=(dSTRmat *)data;
1175 INT i,ic;
1176 INT m=ILU_data->ngrid;
1177 INT nc=ILU_data->nc;
1178 INT nc2=nc*nc;
1179 INT nx=ILU_data->nx;
1180 INT ny=ILU_data->ny;
1181 INT nz=ILU_data->nz;
1182 INT nxy=ILU_data->nxy;
1183 INT size=m*nc;
1184 INT nline, nplane;
1185
1186 if (nx == 1) {
1187 nline = ny;
1188 nplane = m;
1189 }
1190 else if (ny == 1) {
1191 nline = nx;
1192 nplane = m;
1193 }
1194 else if (nz == 1) {
1195 nline = nx;
1196 nplane = m;
1197 }
1198 else {
1199 nline = nx;
1200 nplane = nxy;
1201 }
1202
1203 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
1204
1205 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
1206
1207 zr=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
1208
1209 // copy residual r to zr, to save r
1210 //for (i=0;i<size;++i) zr[i]=r[i];
1211 memcpy(zr,r,(size)*sizeof(REAL));
1212 if (nc == 1) {
1213 // forward sweep: solve unit lower matrix equation L*zz=zr
1214 zz[0]=zr[0];
1215 for (i=1;i<m;++i) {
1216
1217 zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
1218 if (i>=nline-1)
1219 zz[i]=zz[i]-ILU_data->offdiag[2][i-nline+1]*zz[i-nline+1];
1220
1221 if (i>=nline)
1222 zz[i]=zz[i]-ILU_data->offdiag[4][i-nline]*zz[i-nline];
1223 if (i>=nplane-nline)
1224 zz[i]=zz[i]-ILU_data->offdiag[6][i-nplane+nline]*zz[i-nplane+nline];
1225 if (i>=nplane-1)
1226 zz[i]=zz[i]-ILU_data->offdiag[8][i-nplane+1]*zz[i-nplane+1];
1227 if (i>=nplane)
1228 zz[i]=zz[i]-ILU_data->offdiag[10][i-nplane]*zz[i-nplane];
1229 }
1230
1231 } // end if (nc == 1)
1232
1233 else if (nc == 3) {
1234
1235 // forward sweep: solve unit lower matrix equation L*zz=zr
1236 fasp_darray_cp_nc3(&(zr[0]),&(zz[0]));
1237
1238 for (i=1;i<m;++i) {
1239 ic=i*nc;
1240 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
1241 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
1242 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1243
1244 if (i>=nline-1) {
1245 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
1246 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
1247 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1248 }
1249
1250 if (i>=nline) {
1251 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
1252 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
1253 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1254 }
1255
1256 if (i>=nplane-nline) {
1257 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
1258 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
1259 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1260 }
1261
1262 if (i>=nplane-1) {
1263 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
1264 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
1265 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1266 }
1267
1268 if (i>=nplane) {
1269 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
1270 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
1271 fasp_blas_darray_axpy_nc3(-1,tc,&(zr[ic]));
1272 }
1273
1274 fasp_darray_cp_nc3(&(zr[ic]),&(zz[ic]));
1275 }
1276
1277 } // end if (nc == 3)
1278
1279 else if (nc == 5) {
1280
1281 // forward sweep: solve unit lower matrix equation L*zz=zr
1282 fasp_darray_cp_nc5(&(zr[0]),&(zz[0]));
1283
1284 for (i=1;i<m;++i) {
1285 ic=i*nc;
1286 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
1287 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
1288 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1289
1290 if (i>=nline-1) {
1291 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
1292 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
1293 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1294 }
1295
1296 if (i>=nline) {
1297 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
1298 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
1299 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1300 }
1301
1302 if (i>=nplane-nline) {
1303 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
1304 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
1305 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1306 }
1307
1308 if (i>=nplane-1) {
1309 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
1310 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
1311 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1312 }
1313
1314 if (i>=nplane) {
1315 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
1316 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
1317 fasp_blas_darray_axpy_nc5(-1,tc,&(zr[ic]));
1318 }
1319
1320 fasp_darray_cp_nc5(&(zr[ic]),&(zz[ic]));
1321 }
1322
1323 } // end if (nc == 5)
1324
1325 else if (nc == 7) {
1326
1327 // forward sweep: solve unit lower matrix equation L*zz=zr
1328 fasp_darray_cp_nc7(&(zr[0]),&(zz[0]));
1329
1330 for (i=1;i<m;++i) {
1331 ic=i*nc;
1332 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
1333 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc);
1334 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1335
1336 if (i>=nline-1) {
1337 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
1338 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc);
1339 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1340 }
1341
1342 if (i>=nline) {
1343 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
1344 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc);
1345 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1346 }
1347
1348 if (i>=nplane-nline) {
1349 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
1350 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc);
1351 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1352 }
1353
1354 if (i>=nplane-1) {
1355 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
1356 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc);
1357 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1358 }
1359
1360 if (i>=nplane) {
1361 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
1362 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc);
1363 fasp_blas_darray_axpy_nc7(-1,tc,&(zr[ic]));
1364 }
1365
1366 fasp_darray_cp_nc7(&(zr[ic]),&(zz[ic]));
1367 }
1368
1369 } // end if (nc == 7)
1370
1371 else {
1372 // forward sweep: solve unit lower matrix equation L*zz=zr
1373 fasp_darray_cp(nc,&(zr[0]),&(zz[0]));
1374 for (i=1;i<m;++i) {
1375 ic=i*nc;
1376 //zz[i]=zr[i]-ILU_data->offdiag[0][i-1]*zz[i-1];
1377 fasp_blas_smat_mxv(&(ILU_data->offdiag[0][(i-1)*nc2]),&(zz[(i-1)*nc]),tc,nc);
1378 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1379
1380 if (i>=nline-1) {
1381 //zz[i]=zz[i]-ILU_data->offdiag[2][i-nx+1]*zz[i-nx+1];
1382 fasp_blas_smat_mxv(&(ILU_data->offdiag[2][(i-nline+1)*nc2]),&(zz[(i-nline+1)*nc]),tc,nc);
1383 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1384 }
1385
1386 if (i>=nline) {
1387 //zz[i]=zz[i]-ILU_data->offdiag[4][i-nx]*zz[i-nx];
1388 fasp_blas_smat_mxv(&(ILU_data->offdiag[4][(i-nline)*nc2]),&(zz[(i-nline)*nc]),tc,nc);
1389 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1390 }
1391
1392 if (i>=nplane-nline) {
1393 //zz[i]=zz[i]-ILU_data->offdiag[6][i-nxy+nx]*zz[i-nxy+nx];
1394 fasp_blas_smat_mxv(&(ILU_data->offdiag[6][(i-nplane+nline)*nc2]),&(zz[(i-nplane+nline)*nc]),tc,nc);
1395 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1396 }
1397
1398 if (i>=nplane-1) {
1399 // zz[i]=zz[i]-ILU_data->offdiag[8][i-nxy+1]*zz[i-nxy+1];
1400 fasp_blas_smat_mxv(&(ILU_data->offdiag[8][(i-nplane+1)*nc2]),&(zz[(i-nplane+1)*nc]),tc,nc);
1401 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1402 }
1403
1404 if (i>=nplane) {
1405 //zz[i]=zz[i]-ILU_data->offdiag[10][i-nxy]*zz[i-nxy];
1406 fasp_blas_smat_mxv(&(ILU_data->offdiag[10][(i-nplane)*nc2]),&(zz[(i-nplane)*nc]),tc,nc);
1407 fasp_blas_darray_axpy(nc,-1,tc,&(zr[ic]));
1408 }
1409 fasp_darray_cp(nc,&(zr[ic]),&(zz[ic]));
1410 }
1411 } // end else
1412
1413 memcpy(z,zz,(size)*sizeof(REAL));
1414
1415 fasp_mem_free(zr); zr = NULL;
1416 fasp_mem_free(zz); zz = NULL;
1417 fasp_mem_free(tc); tc = NULL;
1418
1419 return;
1420}
1421
1434void fasp_precond_dstr_ilu1_backward (REAL *r,
1435 REAL *z,
1436 void *data)
1437{
1438 REAL *zz,*tc;
1439
1440 dSTRmat *ILU_data=(dSTRmat *)data;
1441 INT i,ic, ic2;
1442 INT m=ILU_data->ngrid;
1443 INT nc=ILU_data->nc;
1444 INT nc2=nc*nc;
1445 INT nx=ILU_data->nx;
1446 INT ny=ILU_data->ny;
1447 INT nz=ILU_data->nz;
1448 INT nxy=ILU_data->nxy;
1449 INT size=m*nc;
1450 INT nline, nplane;
1451
1452 if (nx == 1) {
1453 nline = ny;
1454 nplane = m;
1455 }
1456 else if (ny == 1) {
1457 nline = nx;
1458 nplane = m;
1459 }
1460 else if (nz == 1) {
1461 nline = nx;
1462 nplane = m;
1463 }
1464 else {
1465 nline = nx;
1466 nplane = nxy;
1467 }
1468
1469 tc=(REAL*)fasp_mem_calloc(nc, sizeof(REAL));
1470
1471 zz=(REAL*)fasp_mem_calloc(size, sizeof(REAL));
1472
1473 // copy residual r to zr, to save r
1474 //for (i=0;i<size;++i) zr[i]=r[i];
1475 memcpy(zz,r,(size)*sizeof(REAL));
1476 if (nc == 1) {
1477 // backward sweep: solve upper matrix equation U*z=zz
1478
1479 z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1480 for (i=m-2;i>=0;i--) {
1481
1482 zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1483 if (i+nline-1<m)
1484 zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nline-1];
1485 if (i+nline<m)
1486 zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nline];
1487 if (i+nplane-nline<m)
1488 zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nplane-nline];
1489 if (i+nplane-1<m)
1490 zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nplane-1];
1491 if (i+nplane<m)
1492 zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nplane];
1493
1494 z[i]=ILU_data->diag[i]*zz[i];
1495
1496 }
1497
1498 } // end if (nc == 1)
1499
1500 else if (nc == 3) {
1501 // backward sweep: solve upper matrix equation U*z=zz
1502
1503 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1504 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1505
1506 for (i=m-2;i>=0;i--) {
1507 ic=i*nc;
1508 ic2=ic*nc;
1509
1510 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1511 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1512 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1513
1514 if (i+nline-1<m) {
1515 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
1516 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
1517 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1518 }
1519
1520 if (i+nline<m) {
1521 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
1522 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
1523 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1524 }
1525
1526 if (i+nplane-nline<m) {
1527 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
1528 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
1529 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1530 }
1531
1532 if (i+nplane-1<m) {
1533 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
1534 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
1535 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1536 }
1537
1538 if (i+nplane<m) {
1539 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
1540 fasp_blas_smat_mxv_nc3(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
1541 fasp_blas_darray_axpy_nc3(-1,tc,&(zz[ic]));
1542 }
1543
1544 //z[i]=ILU_data->diag[i]*zz[i];
1545 fasp_blas_smat_mxv_nc3(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1546 } // end for (i=m-2;i>=0;i--)
1547
1548 } // end if (nc == 3)
1549
1550 else if (nc == 5) {
1551 // backward sweep: solve upper matrix equation U*z=zz
1552
1553 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1554 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1555
1556 for (i=m-2;i>=0;i--) {
1557 ic=i*nc;
1558 ic2=ic*nc;
1559
1560 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1561 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1562 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1563
1564 if (i+nline-1<m) {
1565 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
1566 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
1567 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1568 }
1569
1570 if (i+nline<m) {
1571 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
1572 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
1573 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1574 }
1575
1576 if (i+nplane-nline<m) {
1577 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
1578 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
1579 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1580 }
1581
1582 if (i+nplane-1<m) {
1583 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
1584 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
1585 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1586 }
1587
1588 if (i+nplane<m) {
1589 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
1590 fasp_blas_smat_mxv_nc5(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
1591 fasp_blas_darray_axpy_nc5(-1,tc,&(zz[ic]));
1592 }
1593
1594 //z[i]=ILU_data->diag[i]*zz[i];
1595 fasp_blas_smat_mxv_nc5(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1596 } // end for (i=m-2;i>=0;i--)
1597
1598 } // end if (nc == 5)
1599
1600 else if (nc == 7) {
1601 // backward sweep: solve upper matrix equation U*z=zz
1602
1603 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1604 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]));
1605
1606 for (i=m-2;i>=0;i--) {
1607 ic=i*nc;
1608 ic2=ic*nc;
1609
1610 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1611 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc);
1612 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1613
1614 if (i+nline-1<m) {
1615 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
1616 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc);
1617 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1618 }
1619
1620 if (i+nline<m) {
1621 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
1622 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc);
1623 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1624 }
1625
1626 if (i+nplane-nline<m) {
1627 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
1628 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc);
1629 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1630 }
1631
1632 if (i+nplane-1<m) {
1633 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
1634 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc);
1635 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1636 }
1637
1638 if (i+nplane<m) {
1639 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
1640 fasp_blas_smat_mxv_nc7(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc);
1641 fasp_blas_darray_axpy_nc7(-1,tc,&(zz[ic]));
1642 }
1643
1644 //z[i]=ILU_data->diag[i]*zz[i];
1645 fasp_blas_smat_mxv_nc7(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]));
1646 } // end for (i=m-2;i>=0;i--)
1647
1648 } // end if (nc == 7)
1649
1650 else {
1651 // backward sweep: solve upper matrix equation U*z=zz
1652
1653 // z[m-1]=zz[m-1]*ILU_data->diag[m-1];
1654 fasp_blas_smat_mxv(&(ILU_data->diag[(m-1)*nc2]),&(zz[(m-1)*nc]),&(z[(m-1)*nc]),nc);
1655
1656 for (i=m-2;i>=0;i--) {
1657 ic=i*nc;
1658 ic2=ic*nc;
1659 //zz[i]=zz[i]-ILU_data->offdiag[1][i]*z[i+1];
1660 fasp_blas_smat_mxv(&(ILU_data->offdiag[1][ic2]),&(z[(i+1)*nc]),tc,nc);
1661 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1662
1663 if (i+nline-1<m) {
1664 //zz[i]=zz[i]-ILU_data->offdiag[3][i]*z[i+nx-1];
1665 fasp_blas_smat_mxv(&(ILU_data->offdiag[3][ic2]),&(z[(i+nline-1)*nc]),tc,nc);
1666 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1667 }
1668
1669 if (i+nline<m) {
1670 //zz[i]=zz[i]-ILU_data->offdiag[5][i]*z[i+nx];
1671 fasp_blas_smat_mxv(&(ILU_data->offdiag[5][ic2]),&(z[(i+nline)*nc]),tc,nc);
1672 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1673 }
1674
1675 if (i+nplane-nline<m) {
1676 //zz[i]=zz[i]-ILU_data->offdiag[7][i]*z[i+nxy-nx];
1677 fasp_blas_smat_mxv(&(ILU_data->offdiag[7][ic2]),&(z[(i+nplane-nline)*nc]),tc,nc);
1678 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1679 }
1680
1681 if (i+nplane-1<m) {
1682 //zz[i]=zz[i]-ILU_data->offdiag[9][i]*z[i+nxy-1];
1683 fasp_blas_smat_mxv(&(ILU_data->offdiag[9][ic2]),&(z[(i+nplane-1)*nc]),tc,nc);
1684 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1685 }
1686
1687 if (i+nplane<m) {
1688 //zz[i]=zz[i]-ILU_data->offdiag[11][i]*z[i+nxy];
1689 fasp_blas_smat_mxv(&(ILU_data->offdiag[11][ic2]),&(z[(i+nplane)*nc]),tc,nc);
1690 fasp_blas_darray_axpy(nc,-1,tc,&(zz[ic]));
1691 }
1692 //z[i]=ILU_data->diag[i]*zz[i];
1693 fasp_blas_smat_mxv(&(ILU_data->diag[ic2]),&(zz[ic]),&(z[ic]),nc);
1694 }
1695 } // end else
1696
1697 fasp_mem_free(zz); zz = NULL;
1698 fasp_mem_free(tc); tc = NULL;
1699
1700 return;
1701}
1702
1715void fasp_precond_dstr_blockgs (REAL *r,
1716 REAL *z,
1717 void *data)
1718{
1719 precond_data_str *predata=(precond_data_str *)data;
1720 dSTRmat *A = predata->A_str;
1721 dvector *diaginv = predata->diaginv;
1722 ivector *pivot = predata->pivot;
1723 ivector *order = predata->order;
1724 ivector *neigh = predata->neigh;
1725
1726 INT i;
1727 const INT nc = A->nc;
1728 const INT ngrid = A->ngrid;
1729 const INT n = nc*ngrid; // whole size
1730
1731 dvector zz, rr;
1732 zz.row=rr.row=n; zz.val=z; rr.val=r;
1733 fasp_dvec_set(n,&zz,0.0);
1734
1735 for (i=0; i<1; ++i)
1736 fasp_smoother_dstr_swz(A, &rr, &zz, diaginv, pivot, neigh, order);
1737}
1738
1739/*---------------------------------*/
1740/*-- Private Functions --*/
1741/*---------------------------------*/
1742
1754static inline void fasp_darray_cp_nc3 (const REAL *x,
1755 REAL *y)
1756{
1757 memcpy(y, x, 3*sizeof(REAL));
1758}
1759
1771static inline void fasp_darray_cp_nc5 (const REAL *x,
1772 REAL *y)
1773{
1774 memcpy(y, x, 5*sizeof(REAL));
1775}
1776
1788static inline void fasp_darray_cp_nc7 (const REAL *x,
1789 REAL *y)
1790{
1791 memcpy(y, x, 7*sizeof(REAL));
1792}
1793
1794/*---------------------------------*/
1795/*-- End of File --*/
1796/*---------------------------------*/
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_dvec_set
void fasp_dvec_set(INT n, dvector *x, const REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: AuxVector.c:222
fasp_blas_darray_axpy_nc3
void fasp_blas_darray_axpy_nc3(const REAL a, const REAL *x, REAL *y)
y = a*x + y, length of x and y should be 3
Definition: BlaArray.c:255
fasp_blas_darray_axpy_nc7
void fasp_blas_darray_axpy_nc7(const REAL a, const REAL *x, REAL *y)
y = a*x + y, length of x and y should be 7
Definition: BlaArray.c:327
fasp_blas_darray_axpy_nc5
void fasp_blas_darray_axpy_nc5(const REAL a, const REAL *x, REAL *y)
y = a*x + y, length of x and y should be 5
Definition: BlaArray.c:282
fasp_blas_darray_axpy
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:90
fasp_blas_smat_mxv_nc7
void fasp_blas_smat_mxv_nc7(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 7*7 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:200
fasp_blas_smat_mxv_nc5
void fasp_blas_smat_mxv_nc5(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 5*5 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:176
fasp_blas_smat_mxv
void fasp_blas_smat_mxv(const REAL *a, const REAL *b, REAL *c, const INT n)
Compute the product of a small full matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:238
fasp_blas_smat_mxv_nc3
void fasp_blas_smat_mxv_nc3(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 3*3 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:133
fasp_precond_dstr_ilu0_forward
void fasp_precond_dstr_ilu0_forward(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(0) decomposition: Lz = r.
Definition: PreSTR.c:824
fasp_precond_dstr_diag
void fasp_precond_dstr_diag(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreSTR.c:44
fasp_precond_dstr_ilu1_forward
void fasp_precond_dstr_ilu1_forward(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(1) decomposition: Lz = r.
Definition: PreSTR.c:1168
fasp_precond_dstr_blockgs
void fasp_precond_dstr_blockgs(REAL *r, REAL *z, void *data)
CPR-type preconditioner (STR format)
Definition: PreSTR.c:1715
fasp_precond_dstr_ilu1
void fasp_precond_dstr_ilu1(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(1) decomposition.
Definition: PreSTR.c:349
fasp_precond_dstr_ilu0
void fasp_precond_dstr_ilu0(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(0) decomposition.
Definition: PreSTR.c:71
fasp_precond_dstr_ilu1_backward
void fasp_precond_dstr_ilu1_backward(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(1) decomposition: Uz = r.
Definition: PreSTR.c:1434
fasp_precond_dstr_ilu0_backward
void fasp_precond_dstr_ilu0_backward(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(0) decomposition: Uz = r.
Definition: PreSTR.c:987
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
INT
#define INT
Definition: fasp.h:72
ILU_data
Data for ILU setup.
Definition: fasp.h:651
dSTRmat
Structure matrix of REAL type.
Definition: fasp.h:316
dSTRmat::ngrid
INT ngrid
number of grids
Definition: fasp.h:334
dSTRmat::nc
INT nc
size of each block (number of components)
Definition: fasp.h:331
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
precond_data_str
Data for preconditioners in dSTRmat format.
Definition: fasp.h:987
precond_data_str::neigh
ivector * neigh
array to store neighbor information
Definition: fasp.h:1061
precond_data_str::pivot
ivector * pivot
the pivot for the GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:1049
precond_data_str::A_str
dSTRmat * A_str
store the whole reservoir block in STR format
Definition: fasp.h:1038
precond_data_str::order
ivector * order
order for smoothing
Definition: fasp.h:1058
precond_data_str::diaginv
dvector * diaginv
the inverse of the diagonals for GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:1046
precond_diag_str
Data for diagonal preconditioners in dSTRmat format.
Definition: fasp.h:1079
precond_diag_str::nc
INT nc
number of components
Definition: fasp.h:1082
precond_diag_str::diag
dvector diag
diagonal elements
Definition: fasp.h:1085