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