Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
ilu_setup_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 
29  dSTRmat *LU)
30 {
31  // local variables
32  INT i,i1,ix,ixy,ii;
33  INT *LUoffsets;
34  INT nline, nplane;
35 
36  // information of A
37  INT nc = A->nc;
38  INT nc2 = nc*nc;
39  INT nx = A->nx;
40  INT ny = A->ny;
41  INT nz = A->nz;
42  INT nxy = A->nxy;
43  INT ngrid = A->ngrid;
44  INT nband = A->nband;
45 
46  INT *offsets = A->offsets;
47  REAL *smat=(REAL *)fasp_mem_calloc(nc2,sizeof(REAL));
48  REAL *diag = A->diag;
49  REAL *offdiag0=NULL, *offdiag1=NULL, *offdiag2=NULL;
50  REAL *offdiag3=NULL, *offdiag4=NULL, *offdiag5=NULL;
51 
52  // initialize
53  if (nx == 1) {
54  nline = ny;
55  nplane = ngrid;
56  }
57  else if (ny == 1) {
58  nline = nx;
59  nplane = ngrid;
60  }
61  else if (nz == 1) {
62  nline = nx;
63  nplane = ngrid;
64  }
65  else {
66  nline = nx;
67  nplane = nxy;
68  }
69 
70  // check number of bands
71  if (nband == 4) {
72  LUoffsets=(INT *)fasp_mem_calloc(4,sizeof(INT));
73  LUoffsets[0]=-1; LUoffsets[1]=1; LUoffsets[2]=-nline; LUoffsets[3]=nline;
74  }
75  else if (nband == 6) {
76  LUoffsets=(INT *)fasp_mem_calloc(6,sizeof(INT));
77  LUoffsets[0]=-1; LUoffsets[1]=1; LUoffsets[2]=-nline;
78  LUoffsets[3]=nline; LUoffsets[4]=-nplane; LUoffsets[5]=nplane;
79  }
80  else {
81  printf("%s: number of bands for structured ILU is illegal!\n", __FUNCTION__);
82  return;
83  }
84 
85  // allocate memory to store LU decomposition
86  fasp_dstr_alloc(nx, ny, nz, nxy, ngrid, nband, nc, offsets, LU);
87 
88  // copy diagonal
89  memcpy(LU->diag, diag, (ngrid*nc2)*sizeof(REAL));
90 
91  // check offsets and copy off-diagonals
92  for (i=0; i<nband; ++i) {
93  if (offsets[i] == -1) {
94  offdiag0 = A->offdiag[i];
95  memcpy(LU->offdiag[0],offdiag0,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
96  }
97  else if (offsets[i] == 1) {
98  offdiag1 = A->offdiag[i];
99  memcpy(LU->offdiag[1],offdiag1,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
100  }
101  else if (offsets[i] == -nline) {
102  offdiag2 = A->offdiag[i];
103  memcpy(LU->offdiag[2],offdiag2,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
104  }
105  else if (offsets[i] == nline) {
106  offdiag3 = A->offdiag[i];
107  memcpy(LU->offdiag[3],offdiag3,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
108  }
109  else if (offsets[i] == -nplane) {
110  offdiag4 = A->offdiag[i];
111  memcpy(LU->offdiag[4],offdiag4,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
112  }
113  else if (offsets[i] == nplane) {
114  offdiag5 = A->offdiag[i];
115  memcpy(LU->offdiag[5],offdiag5,((ngrid - ABS(offsets[i]))*nc2)*sizeof(REAL));
116  }
117  else {
118  printf("### ERROR: Offsets for structured ILU is illegal!\n");
119  return;
120  }
121  }
122 
123  // Setup
124  if (nc == 1) {
125  LU->diag[0]=1.0/(LU->diag[0]);
126 
127  for(i=1;i<ngrid;++i) {
128  LU->offdiag[0][i-1]=(offdiag0[i-1])*(LU->diag[i-1]);
129  if (i>=nline)
130  LU->offdiag[2][i-nline]=(offdiag2[i-nline])*(LU->diag[i-nline]);
131  if (i>=nplane)
132  LU->offdiag[4][i-nplane]=(offdiag0[i-nplane])*(LU->diag[i-nplane]);
133 
134  LU->diag[i]=diag[i]-(LU->offdiag[0][i-1])*(LU->offdiag[1][i-1]);
135 
136  if (i>=nline)
137  LU->diag[i]=LU->diag[i]-(LU->offdiag[2][i-nline])*(LU->offdiag[3][i-nline]);
138  if (i>=nplane)
139  LU->diag[i]=LU->diag[i]-(LU->offdiag[4][i-nplane])*(LU->offdiag[5][i-nplane]);
140 
141  LU->diag[i]=1.0/(LU->diag[i]);
142 
143  } // end for (i=1; i<ngrid; ++i)
144 
145  } // end if (nc == 1)
146 
147  else if (nc == 3) {
148 
150 
151  for(i=1;i<ngrid;++i) {
152 
153  i1=(i-1)*9;
154  ix=(i-nline)*9;
155  ixy=(i-nplane)*9;
156  ii=i*9;
157 
158  fasp_blas_smat_mul_nc3(&(offdiag0[i1]),&(LU->diag[i1]),&(LU->offdiag[0][i1]));
159 
160  if (i>=nline)
161  fasp_blas_smat_mul_nc3(&(offdiag2[ix]),&(LU->diag[ix]),&(LU->offdiag[2][ix]));
162  if (i>=nplane)
163  fasp_blas_smat_mul_nc3(&(offdiag4[ixy]),&(LU->diag[ixy]),&(LU->offdiag[4][ixy]));
164 
165  fasp_blas_smat_mul_nc3(&(LU->offdiag[0][i1]),&(LU->offdiag[1][i1]),smat);
166 
167  fasp_blas_array_axpyz_nc3(-1,smat,&(diag[ii]),&(LU->diag[ii]));
168 
169  if (i>=nline) {
170  fasp_blas_smat_mul_nc3(&(LU->offdiag[2][ix]),&(LU->offdiag[3][ix]),smat);
171  fasp_blas_array_axpy_nc3(-1.0,smat,&(LU->diag[ii]));
172  } //end if (i>=nline)
173 
174  if (i>=nplane) {
175  fasp_blas_smat_mul_nc3(&(LU->offdiag[4][ixy]),&(LU->offdiag[5][ixy]),smat);
176  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ii]));
177  } // end if (i>=nplane)
178 
179  fasp_blas_smat_inv_nc3(&(LU->diag[ii]));
180 
181  } // end for(i=1;i<A->ngrid;++i)
182 
183  } // end if (nc == 3)
184 
185  else if (nc == 5) {
186 
188 
189  for(i=1;i<ngrid;++i) {
190 
191  i1=(i-1)*25;
192  ix=(i-nline)*25;
193  ixy=(i-nplane)*25;
194  ii=i*25;
195 
196  fasp_blas_smat_mul_nc5(&(offdiag0[i1]),&(LU->diag[i1]),&(LU->offdiag[0][i1]));
197 
198  if (i>=nline)
199  fasp_blas_smat_mul_nc5(&(offdiag2[ix]),&(LU->diag[ix]),&(LU->offdiag[2][ix]));
200  if (i>=nplane)
201  fasp_blas_smat_mul_nc5(&(offdiag4[ixy]),&(LU->diag[ixy]),&(LU->offdiag[4][ixy]));
202 
203  fasp_blas_smat_mul_nc5(&(LU->offdiag[0][i1]),&(LU->offdiag[1][i1]),smat);
204 
205  fasp_blas_array_axpyz_nc5(-1.0,smat,&(diag[ii]),&(LU->diag[ii]));
206 
207  if (i>=nline) {
208  fasp_blas_smat_mul_nc5(&(LU->offdiag[2][ix]),&(LU->offdiag[3][ix]),smat);
209  fasp_blas_array_axpy_nc5(-1.0,smat,&(LU->diag[ii]));
210  } //end if (i>=nline)
211 
212  if (i>=nplane) {
213  fasp_blas_smat_mul_nc5(&(LU->offdiag[4][ixy]),&(LU->offdiag[5][ixy]),smat);
214  fasp_blas_array_axpy_nc5(-1.0,smat,&(LU->diag[ii]));
215  } // end if (i>=nplane)
216 
217  fasp_blas_smat_inv_nc5(&(LU->diag[ii]));
218 
219  } // end for(i=1;i<A->ngrid;++i)
220 
221  } // end if (nc == 5)
222 
223  else if (nc == 7) {
224 
226 
227  for(i=1;i<ngrid;++i) {
228 
229  i1=(i-1)*49;
230  ix=(i-nline)*49;
231  ixy=(i-nplane)*49;
232  ii=i*49;
233 
234  fasp_blas_smat_mul_nc7(&(offdiag0[i1]),&(LU->diag[i1]),&(LU->offdiag[0][i1]));
235 
236  if (i>=nline)
237  fasp_blas_smat_mul_nc7(&(offdiag2[ix]),&(LU->diag[ix]),&(LU->offdiag[2][ix]));
238  if (i>=nplane)
239  fasp_blas_smat_mul_nc7(&(offdiag4[ixy]),&(LU->diag[ixy]),&(LU->offdiag[4][ixy]));
240 
241  fasp_blas_smat_mul_nc7(&(LU->offdiag[0][i1]),&(LU->offdiag[1][i1]),smat);
242 
243  fasp_blas_array_axpyz_nc7(-1.0,smat,&(diag[ii]),&(LU->diag[ii]));
244 
245  if (i>=nline) {
246  fasp_blas_smat_mul_nc7(&(LU->offdiag[2][ix]),&(LU->offdiag[3][ix]),smat);
247  fasp_blas_array_axpy_nc7(-1.0,smat,&(LU->diag[ii]));
248  } //end if (i>=nline)
249 
250  if (i>=nplane) {
251  fasp_blas_smat_mul_nc7(&(LU->offdiag[4][ixy]),&(LU->offdiag[5][ixy]),smat);
252  fasp_blas_array_axpy_nc7(-1.0,smat,&(LU->diag[ii]));
253  } // end if (i>=nplane)
254 
255  fasp_blas_smat_inv_nc7(&(LU->diag[ii]));
256 
257  } // end for(i=1;i<A->ngrid;++i)
258 
259  } // end if (nc == 7)
260 
261  else {
262  fasp_blas_smat_inv(LU->diag,nc);
263 
264  for(i=1;i<ngrid;++i) {
265 
266  i1=(i-1)*nc2;
267  ix=(i-nline)*nc2;
268  ixy=(i-nplane)*nc2;
269  ii=i*nc2;
270 
271  fasp_blas_smat_mul(&(offdiag0[i1]),&(LU->diag[i1]),&(LU->offdiag[0][i1]),nc);
272 
273  if (i>=nline)
274  fasp_blas_smat_mul(&(offdiag2[ix]),&(LU->diag[ix]),&(LU->offdiag[2][ix]),nc);
275  if (i>=nplane)
276  fasp_blas_smat_mul(&(offdiag4[ixy]),&(LU->diag[ixy]),&(LU->offdiag[4][ixy]),nc);
277 
278  fasp_blas_smat_mul(&(LU->offdiag[0][i1]),&(LU->offdiag[1][i1]),smat,nc);
279 
280  fasp_blas_array_axpyz(nc2,-1,smat,&(diag[ii]),&(LU->diag[ii]));
281 
282  if (i>=nline) {
283  fasp_blas_smat_mul(&(LU->offdiag[2][ix]),&(LU->offdiag[3][ix]),smat,nc);
284  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ii]));
285  } //end if (i>=nline)
286 
287  if (i>=nplane) {
288  fasp_blas_smat_mul(&(LU->offdiag[4][ixy]),&(LU->offdiag[5][ixy]),smat,nc);
289  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ii]));
290  } // end if (i>=nplane)
291 
292  fasp_blas_smat_inv(&(LU->diag[ii]),nc);
293 
294  } // end for(i=1;i<A->ngrid;++i)
295 
296  }
297 
298  fasp_mem_free(smat);
299 
300  return;
301 }
302 
320  dSTRmat *LU)
321 {
322  INT LUnband=12;
323 
324  INT i,j,i1,ix,ixy,ixyx,ix1,ixy1,ic,i1c,ixc,ix1c,ixyc,ixy1c,ixyxc;
325  register REAL *smat,t,*tc;
326 
327  const INT nc=A->nc, nc2=nc*nc;
328  const INT nx=A->nx;
329  const INT ny=A->ny;
330  const INT nz=A->nz;
331  const INT nxy=A->nxy;
332  const INT nband=A->nband;
333  const INT ngrid=A->ngrid;
334  INT nline, nplane;
335 
336  if (nx == 1) {
337  nline = ny;
338  nplane = ngrid;
339  }
340  else if (ny == 1) {
341  nline = nx;
342  nplane = ngrid;
343  }
344  else if (nz == 1) {
345  nline = nx;
346  nplane = ngrid;
347  }
348  else {
349  nline = nx;
350  nplane = nxy;
351  }
352 
353  smat=(REAL *)fasp_mem_calloc(nc2,sizeof(REAL));
354 
355  tc=(REAL *)fasp_mem_calloc(nc2,sizeof(REAL));
356 
357  INT LUoffsets[]={-1,1,1-nline,nline-1,-nline,nline,nline-nplane,nplane-nline,1-nplane,nplane-1,-nplane,nplane};
358 
359  fasp_dstr_alloc(nx,A->ny,A->nz,nxy,ngrid,LUnband,nc,LUoffsets,LU);
360 
361  if (nband == 6) memcpy(LU->offdiag[11],A->offdiag[5],((ngrid-nxy)*nc2)*sizeof(REAL));
362  memcpy(LU->diag,A->diag,nc2*sizeof(REAL));
363 
364  if (nc == 1) {
365  //comput the first row
366  LU->diag[0]=1.0/(LU->diag[0]);
367  LU->offdiag[1][0]=A->offdiag[1][0];
368  LU->offdiag[5][0]=A->offdiag[3][0];
369  LU->offdiag[3][0]=0;
370  LU->offdiag[7][0]=0;
371  LU->offdiag[9][0]=0;
372 
373  for(i=1;i<ngrid;++i) {
374  i1=i-1;ix=i-nline;ixy=i-nplane;ix1=ix+1;ixyx=ixy+nline;ixy1=ixy+1;
375 
376  //comput alpha6[i-nxy]
377  if (ixy>=0)
378  LU->offdiag[10][ixy]=A->offdiag[4][ixy]*LU->diag[ixy];
379 
380  //comput alpha5[ixy1]
381  if (ixy1>=0) {
382  t=0;
383 
384  if (ixy>=0) t-=LU->offdiag[10][ixy]*LU->offdiag[1][ixy];
385 
386  LU->offdiag[8][ixy1]=t*(LU->diag[ixy1]);
387  }
388 
389  //comput alpha4[ixyx]
390  if (ixyx>=0) {
391  t=0;
392 
393  if (ixy>=0) t-=LU->offdiag[10][ixy]*LU->offdiag[5][ixy];
394  if (ixy1>=0) t-=LU->offdiag[8][ixy1]*LU->offdiag[3][ixy1];
395 
396  LU->offdiag[6][ixyx]=t*(LU->diag[ixyx]);
397  }
398 
399  //comput alpha3[ix]
400  if (ix>=0) {
401  t=A->offdiag[2][ix];
402 
403  if (ixy>=0) t-=LU->offdiag[10][ixy]*LU->offdiag[7][ixy];
404 
405  LU->offdiag[4][ix]=t*(LU->diag[ix]);
406  }
407 
408  //comput alpha2[i-nx+1]
409  if (ix1>=0) {
410  t=0;
411 
412  if (ix>=0) t-=LU->offdiag[4][ix]*LU->offdiag[1][ix];
413  if (ixy1>=0) t-=LU->offdiag[8][ixy1]*LU->offdiag[7][ixy1];
414 
415  LU->offdiag[2][ix1]=t*(LU->diag[ix1]);
416  }
417 
418  //comput alpha1[i-1]
419  t=A->offdiag[0][i1];
420 
421  if (ix>=0) t-=LU->offdiag[4][ix]*LU->offdiag[3][ix];
422  if (ixy>=0) t-=LU->offdiag[10][ixy]*LU->offdiag[9][ixy];
423 
424  LU->offdiag[0][i1]=t*(LU->diag[i1]);
425 
426  //comput beta1[i]
427  if (i+1<ngrid) {
428  t=A->offdiag[1][i];
429 
430  if (ix1>=0) t-=LU->offdiag[2][ix1]*LU->offdiag[5][ix1];
431  if (ixy1>=0) t-=LU->offdiag[8][ixy1]*LU->offdiag[11][ixy1];
432 
433  LU->offdiag[1][i]=t;
434  }
435 
436  //comput beta2[i]
437  if (i+nline-1<ngrid) {
438  t=-LU->offdiag[0][i1]*LU->offdiag[5][i1];
439 
440  if (ixyx>=0) t-=LU->offdiag[6][ixyx]*LU->offdiag[9][ixyx];
441 
442  LU->offdiag[3][i]=t;
443  }
444 
445  //comput beta3[i]
446  if (i+nline<ngrid) {
447  t=A->offdiag[3][i];
448 
449  if (ixyx>=0) t-=LU->offdiag[6][ixyx]*LU->offdiag[11][ixyx];
450 
451  LU->offdiag[5][i]=t;
452  }
453 
454 
455  //comput beta4[i]
456  if (i+nplane-nline<ngrid) {
457  t=0;
458 
459  if (ix1>=0) t-=LU->offdiag[2][ix1]*LU->offdiag[9][ix1];
460  if (ix>=0) t-=LU->offdiag[4][ix]*LU->offdiag[11][ix];
461 
462  LU->offdiag[7][i]=t;
463  }
464 
465  //comput beta5[i]
466  if (i+nplane-1<ngrid) LU->offdiag[9][i]=-LU->offdiag[0][i1]*LU->offdiag[11][i1];
467 
468  //comput d[i]
469  LU->diag[i]=A->diag[i]-(LU->offdiag[0][i1])*(LU->offdiag[1][i1]);
470 
471  if (ix1>=0) LU->diag[i]-=(LU->offdiag[2][ix1])*(LU->offdiag[3][ix1]);
472  if (ix>=0) LU->diag[i]-=(LU->offdiag[4][ix])*(LU->offdiag[5][ix]);
473  if (ixyx>=0) LU->diag[i]-=(LU->offdiag[6][ixyx])*(LU->offdiag[7][ixyx]);
474  if (ixy1>=0) LU->diag[i]-=(LU->offdiag[8][ixy1])*(LU->offdiag[9][ixy1]);
475  if (ixy>=0) LU->diag[i]-=(LU->offdiag[10][ixy])*(LU->offdiag[11][ixy]);
476 
477  LU->diag[i]=1.0/(LU->diag[i]);
478 
479  } // end for (i=1; i<ngrid; ++i)
480 
481  } // end if (nc == 1)
482 
483  else if (nc == 3) {
484  //comput the first row
486  memcpy(LU->offdiag[1],A->offdiag[1],9*sizeof(REAL));
487  memcpy(LU->offdiag[5],A->offdiag[3],9*sizeof(REAL));
488 
489  for(i=1;i<ngrid;++i) {
490  i1=i-1;ix=i-nline;ixy=i-nplane;ix1=ix+1;ixyx=ixy+nline;ixy1=ixy+1;
491  ic=i*nc2;i1c=i1*nc2;ixc=ix*nc2;ix1c=ix1*nc2;ixyc=ixy*nc2;ixy1c=ixy1*nc2;ixyxc=ixyx*nc2;
492 
493  //comput alpha6[i-nxy]
494  if (ixy>=0) fasp_blas_smat_mul_nc3(&(A->offdiag[4][ixyc]),&(LU->diag[ixyc]),&(LU->offdiag[10][ixyc]));
495 
496  //comput alpha5[ixy1]
497  if (ixy1>=0) {
498  for (j=0;j<9;++j) tc[j]=0;
499 
500  if (ixy>=0) {
501  fasp_blas_smat_mul_nc3(&(LU->offdiag[10][ixyc]),&(LU->offdiag[1][ixyc]),smat);
502  fasp_blas_array_axpy_nc3(-1,smat,tc);
503  }
504 
505  fasp_blas_smat_mul_nc3(tc,&(LU->diag[ixy1c]),&(LU->offdiag[8][ixy1c]));
506  }
507 
508  //comput alpha4[ixyx]
509  if (ixyx>=0) {
510  for (j=0;j<9;++j) tc[j]=0;
511 
512  if (ixy>=0) {
513  fasp_blas_smat_mul_nc3(&(LU->offdiag[10][ixyc]),&(LU->offdiag[5][ixyc]),smat);
514  fasp_blas_array_axpy_nc3(-1,smat,tc);
515  }
516 
517  if (ixy1>=0) {
518  fasp_blas_smat_mul_nc3(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[3][ixy1c]),smat);
519  fasp_blas_array_axpy_nc3(-1,smat,tc);
520  }
521 
522  fasp_blas_smat_mul_nc3(tc,&(LU->diag[ixyxc]),&(LU->offdiag[6][ixyxc]));
523  }
524 
525  //comput alpha3[ix]
526  if (ix>=0) {
527 
528  memcpy(tc,&(A->offdiag[2][ixc]),9*sizeof(REAL));
529 
530  if (ixy>=0) {
531  fasp_blas_smat_mul_nc3(&(LU->offdiag[10][ixyc]),&(LU->offdiag[7][ixyc]),smat);
532  fasp_blas_array_axpy_nc3(-1,smat,tc);
533  }
534 
535  fasp_blas_smat_mul_nc3(tc,&(LU->diag[ixc]),&(LU->offdiag[4][ixc]));
536  }
537 
538  //comput alpha2[i-nx+1]
539  if (ix1>=0) {
540 
541  for (j=0;j<9;++j) tc[j]=0;
542 
543  if (ix>=0) {
544  fasp_blas_smat_mul_nc3(&(LU->offdiag[4][ixc]),&(LU->offdiag[1][ixc]),smat);
545  fasp_blas_array_axpy_nc3(-1,smat,tc);
546  }
547 
548  if (ixy1>=0) {
549  fasp_blas_smat_mul_nc3(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[7][ixy1c]),smat);
550  fasp_blas_array_axpy_nc3(-1,smat,tc);
551  }
552 
553  fasp_blas_smat_mul_nc3(tc,&(LU->diag[ix1c]),&(LU->offdiag[2][ix1c]));
554 
555  } // end if (ix1 >= 0)
556 
557  //comput alpha1[i-1]
558 
559  memcpy(tc,&(A->offdiag[0][i1c]),9*sizeof(REAL));
560 
561  if (ix>=0) {
562  fasp_blas_smat_mul_nc3(&(LU->offdiag[4][ixc]),&(LU->offdiag[3][ixc]),smat);
563  fasp_blas_array_axpy_nc3(-1,smat,tc);
564  }
565 
566  if (ixy>=0) {
567  fasp_blas_smat_mul_nc3(&(LU->offdiag[10][ixyc]),&(LU->offdiag[9][ixyc]),smat);
568  fasp_blas_array_axpy_nc3(-1,smat,tc);
569  }
570 
571  fasp_blas_smat_mul_nc3(tc,&(LU->diag[i1c]),&(LU->offdiag[0][i1c]));
572 
573  //comput beta1[i]
574  if (i+1<ngrid) {
575 
576  memcpy(&(LU->offdiag[1][ic]),&(A->offdiag[1][ic]),9*sizeof(REAL));
577 
578  if (ix1>=0) {
579  fasp_blas_smat_mul_nc3(&(LU->offdiag[2][ix1c]),&(LU->offdiag[5][ix1c]),smat);
580  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[1][ic]));
581  }
582 
583  if (ixy1>=0) {
584  fasp_blas_smat_mul_nc3(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[11][ixy1c]),smat);
585  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[1][ic]));
586  }
587 
588  }
589 
590  //comput beta2[i]
591  if (i+nline-1<ngrid) {
592 
593  {
594  fasp_blas_smat_mul_nc3(&(LU->offdiag[0][i1c]),&(LU->offdiag[5][i1c]),smat);
595  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[3][ic]));
596  }
597 
598  if (ixyx>=0) {
599  fasp_blas_smat_mul_nc3(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[9][ixyxc]),smat);
600  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[3][ic]));
601  }
602 
603  }
604 
605  //comput beta3[i]
606  if (i+nline<ngrid) {
607 
608  memcpy(&(LU->offdiag[5][ic]),&(A->offdiag[3][ic]),9*sizeof(REAL));
609 
610  if (ixyx>=0) {
611  fasp_blas_smat_mul_nc3(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[11][ixyxc]),smat);
612  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[5][ic]));
613  }
614 
615  }
616 
617  //comput beta4[i]
618  if (i+nplane-nline<ngrid) {
619 
620  if (ix1>=0) {
621  fasp_blas_smat_mul_nc3(&(LU->offdiag[2][ix1c]),&(LU->offdiag[9][ix1c]),smat);
622  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[7][ic]));
623  }
624 
625  if (ix>=0) {
626  fasp_blas_smat_mul_nc3(&(LU->offdiag[4][ixc]),&(LU->offdiag[11][ixc]),smat);
627  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[7][ic]));
628  }
629 
630  }
631 
632  //comput beta5[i]
633  if (i+nplane-1<ngrid) {
634  fasp_blas_smat_mul_nc3(&(LU->offdiag[0][i1c]),&(LU->offdiag[11][i1c]),smat);
635  fasp_blas_array_axpy_nc3(-1,smat,&(LU->offdiag[9][ic]));
636  }
637 
638  //comput d[i]
639  {
640  fasp_blas_smat_mul_nc3(&(LU->offdiag[0][i1c]),&(LU->offdiag[1][i1c]),smat);
641  fasp_blas_array_axpyz_nc3(-1,smat,&(A->diag[ic]),&(LU->diag[ic]));
642  }
643 
644  if (ix1>=0) {
645  fasp_blas_smat_mul_nc3(&(LU->offdiag[2][ix1c]),&(LU->offdiag[3][ix1c]),smat);
646  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ic]));
647  }
648 
649  if (ix>=0) {
650  fasp_blas_smat_mul_nc3(&(LU->offdiag[4][ixc]),&(LU->offdiag[5][ixc]),smat);
651  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ic]));
652  }
653 
654  if (ixyx>=0) {
655  fasp_blas_smat_mul_nc3(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[7][ixyxc]),smat);
656  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ic]));
657  }
658 
659  if (ixy1>=0) {
660  fasp_blas_smat_mul_nc3(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[9][ixy1c]),smat);
661  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ic]));
662  }
663 
664  if (ixy>=0) {
665  fasp_blas_smat_mul_nc3(&(LU->offdiag[10][ixyc]),&(LU->offdiag[11][ixyc]),smat);
666  fasp_blas_array_axpy_nc3(-1,smat,&(LU->diag[ic]));
667  }
668 
669  fasp_blas_smat_inv_nc3(&(LU->diag[ic]));
670 
671  } // end for(i=1;i<ngrid;++i)
672 
673  } // end if (nc == 3)
674 
675  else if (nc == 5) {
676  //comput the first row
677  //fasp_blas_smat_inv_nc5(LU->diag);
678  fasp_blas_smat_inv(LU->diag,5);
679  memcpy(LU->offdiag[1],A->offdiag[1], 25*sizeof(REAL));
680  memcpy(LU->offdiag[5],A->offdiag[3], 25*sizeof(REAL));
681 
682  for(i=1;i<ngrid;++i) {
683  i1=i-1;ix=i-nline;ixy=i-nplane;ix1=ix+1;ixyx=ixy+nline;ixy1=ixy+1;
684  ic=i*nc2;i1c=i1*nc2;ixc=ix*nc2;ix1c=ix1*nc2;ixyc=ixy*nc2;ixy1c=ixy1*nc2;ixyxc=ixyx*nc2;
685 
686  //comput alpha6[i-nxy]
687  if (ixy>=0) fasp_blas_smat_mul_nc5(&(A->offdiag[4][ixyc]),&(LU->diag[ixyc]),&(LU->offdiag[10][ixyc]));
688 
689  //comput alpha5[ixy1]
690  if (ixy1>=0) {
691  for (j=0;j<25;++j) tc[j]=0;
692 
693  if (ixy>=0) {
694  fasp_blas_smat_mul_nc5(&(LU->offdiag[10][ixyc]),&(LU->offdiag[1][ixyc]),smat);
695  fasp_blas_array_axpy_nc5(-1.0,smat,tc);
696  }
697 
698  fasp_blas_smat_mul_nc5(tc,&(LU->diag[ixy1c]),&(LU->offdiag[8][ixy1c]));
699  }
700 
701  //comput alpha4[ixyx]
702  if (ixyx>=0) {
703  for (j=0;j<25;++j) tc[j]=0;
704 
705  if (ixy>=0) {
706  fasp_blas_smat_mul_nc5(&(LU->offdiag[10][ixyc]),&(LU->offdiag[5][ixyc]),smat);
707  fasp_blas_array_axpy_nc5(-1,smat,tc);
708  }
709 
710  if (ixy1>=0) {
711  fasp_blas_smat_mul_nc5(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[3][ixy1c]),smat);
712  fasp_blas_array_axpy_nc5(-1,smat,tc);
713  }
714 
715  fasp_blas_smat_mul_nc5(tc,&(LU->diag[ixyxc]),&(LU->offdiag[6][ixyxc]));
716  }
717 
718  //comput alpha3[ix]
719  if (ix>=0) {
720 
721  memcpy(tc,&(A->offdiag[2][ixc]),25*sizeof(REAL));
722 
723  if (ixy>=0) {
724  fasp_blas_smat_mul_nc5(&(LU->offdiag[10][ixyc]),&(LU->offdiag[7][ixyc]),smat);
725  fasp_blas_array_axpy_nc5(-1,smat,tc);
726  }
727 
728  fasp_blas_smat_mul_nc5(tc,&(LU->diag[ixc]),&(LU->offdiag[4][ixc]));
729  }
730 
731  //comput alpha2[i-nx+1]
732  if (ix1>=0) {
733 
734  for (j=0;j<25;++j) tc[j]=0;
735 
736  if (ix>=0) {
737  fasp_blas_smat_mul_nc5(&(LU->offdiag[4][ixc]),&(LU->offdiag[1][ixc]),smat);
738  fasp_blas_array_axpy_nc5(-1,smat,tc);
739  }
740 
741  if (ixy1>=0) {
742  fasp_blas_smat_mul_nc5(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[7][ixy1c]),smat);
743  fasp_blas_array_axpy_nc5(-1,smat,tc);
744  }
745 
746  fasp_blas_smat_mul_nc5(tc,&(LU->diag[ix1c]),&(LU->offdiag[2][ix1c]));
747 
748  } // end if (ix1 >= 0)
749 
750  //comput alpha1[i-1]
751 
752  memcpy(tc,&(A->offdiag[0][i1c]), 25*sizeof(REAL));
753 
754  if (ix>=0) {
755  fasp_blas_smat_mul_nc5(&(LU->offdiag[4][ixc]),&(LU->offdiag[3][ixc]),smat);
756  fasp_blas_array_axpy_nc5(-1,smat,tc);
757  }
758 
759  if (ixy>=0) {
760  fasp_blas_smat_mul_nc5(&(LU->offdiag[10][ixyc]),&(LU->offdiag[9][ixyc]),smat);
761  fasp_blas_array_axpy_nc5(-1,smat,tc);
762  }
763 
764  fasp_blas_smat_mul_nc5(tc,&(LU->diag[i1c]),&(LU->offdiag[0][i1c]));
765 
766  //comput beta1[i]
767  if (i+1<ngrid) {
768 
769  memcpy(&(LU->offdiag[1][ic]),&(A->offdiag[1][ic]), 25*sizeof(REAL));
770 
771  if (ix1>=0) {
772  fasp_blas_smat_mul_nc5(&(LU->offdiag[2][ix1c]),&(LU->offdiag[5][ix1c]),smat);
773  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[1][ic]));
774  }
775 
776  if (ixy1>=0) {
777  fasp_blas_smat_mul_nc5(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[11][ixy1c]),smat);
778  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[1][ic]));
779  }
780 
781  }
782 
783  //comput beta2[i]
784  if (i+nline-1<ngrid) {
785 
786  {
787  fasp_blas_smat_mul_nc5(&(LU->offdiag[0][i1c]),&(LU->offdiag[5][i1c]),smat);
788  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[3][ic]));
789  }
790 
791  if (ixyx>=0) {
792  fasp_blas_smat_mul_nc5(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[9][ixyxc]),smat);
793  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[3][ic]));
794  }
795 
796  }
797 
798  //comput beta3[i]
799  if (i+nline<ngrid) {
800 
801  memcpy(&(LU->offdiag[5][ic]),&(A->offdiag[3][ic]), 25*sizeof(REAL));
802 
803  if (ixyx>=0) {
804  fasp_blas_smat_mul_nc5(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[11][ixyxc]),smat);
805  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[5][ic]));
806  }
807 
808  }
809 
810  //comput beta4[i]
811  if (i+nplane-nline<ngrid) {
812 
813  if (ix1>=0) {
814  fasp_blas_smat_mul_nc5(&(LU->offdiag[2][ix1c]),&(LU->offdiag[9][ix1c]),smat);
815  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[7][ic]));
816  }
817 
818  if (ix>=0) {
819  fasp_blas_smat_mul_nc5(&(LU->offdiag[4][ixc]),&(LU->offdiag[11][ixc]),smat);
820  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[7][ic]));
821  }
822 
823  }
824 
825  //comput beta5[i]
826  if (i+nplane-1<ngrid) {
827  fasp_blas_smat_mul_nc5(&(LU->offdiag[0][i1c]),&(LU->offdiag[11][i1c]),smat);
828  fasp_blas_array_axpy_nc5(-1,smat,&(LU->offdiag[9][ic]));
829  }
830 
831  //comput d[i]
832  {
833  fasp_blas_smat_mul_nc5(&(LU->offdiag[0][i1c]),&(LU->offdiag[1][i1c]),smat);
834  fasp_blas_array_axpyz_nc5(-1,smat,&(A->diag[ic]),&(LU->diag[ic]));
835  }
836 
837  if (ix1>=0) {
838  fasp_blas_smat_mul_nc5(&(LU->offdiag[2][ix1c]),&(LU->offdiag[3][ix1c]),smat);
839  fasp_blas_array_axpy_nc5(-1,smat,&(LU->diag[ic]));
840  }
841 
842  if (ix>=0) {
843  fasp_blas_smat_mul_nc5(&(LU->offdiag[4][ixc]),&(LU->offdiag[5][ixc]),smat);
844  fasp_blas_array_axpy_nc5(-1,smat,&(LU->diag[ic]));
845  }
846 
847  if (ixyx>=0) {
848  fasp_blas_smat_mul_nc5(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[7][ixyxc]),smat);
849  fasp_blas_array_axpy_nc5(-1,smat,&(LU->diag[ic]));
850  }
851 
852  if (ixy1>=0) {
853  fasp_blas_smat_mul_nc5(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[9][ixy1c]),smat);
854  fasp_blas_array_axpy_nc5(-1,smat,&(LU->diag[ic]));
855  }
856 
857  if (ixy>=0) {
858  fasp_blas_smat_mul_nc5(&(LU->offdiag[10][ixyc]),&(LU->offdiag[11][ixyc]),smat);
859  fasp_blas_array_axpy_nc5(-1,smat,&(LU->diag[ic]));
860  }
861 
862  //fasp_blas_smat_inv_nc5(&(LU->diag[ic]));
863  fasp_blas_smat_inv(&(LU->diag[ic]), 5);
864 
865  } // end for(i=1;i<ngrid;++i)
866 
867  } // end if (nc == 5)
868 
869  else if (nc == 7) {
870  //comput the first row
871  //fasp_blas_smat_inv_nc5(LU->diag);
872  fasp_blas_smat_inv(LU->diag,7);
873  memcpy(LU->offdiag[1],A->offdiag[1], 49*sizeof(REAL));
874  memcpy(LU->offdiag[5],A->offdiag[3], 49*sizeof(REAL));
875 
876  for(i=1;i<ngrid;++i) {
877  i1=i-1;ix=i-nline;ixy=i-nplane;ix1=ix+1;ixyx=ixy+nline;ixy1=ixy+1;
878  ic=i*nc2;i1c=i1*nc2;ixc=ix*nc2;ix1c=ix1*nc2;ixyc=ixy*nc2;ixy1c=ixy1*nc2;ixyxc=ixyx*nc2;
879 
880  //comput alpha6[i-nxy]
881  if (ixy>=0) fasp_blas_smat_mul_nc7(&(A->offdiag[4][ixyc]),&(LU->diag[ixyc]),&(LU->offdiag[10][ixyc]));
882 
883  //comput alpha5[ixy1]
884  if (ixy1>=0) {
885  for (j=0;j<49;++j) tc[j]=0;
886 
887  if (ixy>=0) {
888  fasp_blas_smat_mul_nc7(&(LU->offdiag[10][ixyc]),&(LU->offdiag[1][ixyc]),smat);
889  fasp_blas_array_axpy_nc7(-1.0,smat,tc);
890  }
891 
892  fasp_blas_smat_mul_nc7(tc,&(LU->diag[ixy1c]),&(LU->offdiag[8][ixy1c]));
893  }
894 
895  //comput alpha4[ixyx]
896  if (ixyx>=0) {
897  for (j=0;j<49;++j) tc[j]=0;
898 
899  if (ixy>=0) {
900  fasp_blas_smat_mul_nc7(&(LU->offdiag[10][ixyc]),&(LU->offdiag[5][ixyc]),smat);
901  fasp_blas_array_axpy_nc7(-1,smat,tc);
902  }
903 
904  if (ixy1>=0) {
905  fasp_blas_smat_mul_nc7(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[3][ixy1c]),smat);
906  fasp_blas_array_axpy_nc7(-1,smat,tc);
907  }
908 
909  fasp_blas_smat_mul_nc7(tc,&(LU->diag[ixyxc]),&(LU->offdiag[6][ixyxc]));
910  }
911 
912  //comput alpha3[ix]
913  if (ix>=0) {
914 
915  memcpy(tc,&(A->offdiag[2][ixc]),49*sizeof(REAL));
916 
917  if (ixy>=0) {
918  fasp_blas_smat_mul_nc7(&(LU->offdiag[10][ixyc]),&(LU->offdiag[7][ixyc]),smat);
919  fasp_blas_array_axpy_nc7(-1,smat,tc);
920  }
921 
922  fasp_blas_smat_mul_nc7(tc,&(LU->diag[ixc]),&(LU->offdiag[4][ixc]));
923  }
924 
925  //comput alpha2[i-nx+1]
926  if (ix1>=0) {
927 
928  for (j=0;j<49;++j) tc[j]=0;
929 
930  if (ix>=0) {
931  fasp_blas_smat_mul_nc7(&(LU->offdiag[4][ixc]),&(LU->offdiag[1][ixc]),smat);
932  fasp_blas_array_axpy_nc7(-1,smat,tc);
933  }
934 
935  if (ixy1>=0) {
936  fasp_blas_smat_mul_nc7(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[7][ixy1c]),smat);
937  fasp_blas_array_axpy_nc7(-1,smat,tc);
938  }
939 
940  fasp_blas_smat_mul_nc7(tc,&(LU->diag[ix1c]),&(LU->offdiag[2][ix1c]));
941 
942  } // end if (ix1 >= 0)
943 
944  //comput alpha1[i-1]
945 
946  memcpy(tc,&(A->offdiag[0][i1c]), 49*sizeof(REAL));
947 
948  if (ix>=0) {
949  fasp_blas_smat_mul_nc7(&(LU->offdiag[4][ixc]),&(LU->offdiag[3][ixc]),smat);
950  fasp_blas_array_axpy_nc7(-1,smat,tc);
951  }
952 
953  if (ixy>=0) {
954  fasp_blas_smat_mul_nc7(&(LU->offdiag[10][ixyc]),&(LU->offdiag[9][ixyc]),smat);
955  fasp_blas_array_axpy_nc7(-1,smat,tc);
956  }
957 
958  fasp_blas_smat_mul_nc7(tc,&(LU->diag[i1c]),&(LU->offdiag[0][i1c]));
959 
960  //comput beta1[i]
961  if (i+1<ngrid) {
962 
963  memcpy(&(LU->offdiag[1][ic]),&(A->offdiag[1][ic]), 49*sizeof(REAL));
964 
965  if (ix1>=0) {
966  fasp_blas_smat_mul_nc7(&(LU->offdiag[2][ix1c]),&(LU->offdiag[5][ix1c]),smat);
967  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[1][ic]));
968  }
969 
970  if (ixy1>=0) {
971  fasp_blas_smat_mul_nc7(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[11][ixy1c]),smat);
972  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[1][ic]));
973  }
974 
975  }
976 
977  //comput beta2[i]
978  if (i+nline-1<ngrid) {
979 
980  {
981  fasp_blas_smat_mul_nc7(&(LU->offdiag[0][i1c]),&(LU->offdiag[5][i1c]),smat);
982  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[3][ic]));
983  }
984 
985  if (ixyx>=0) {
986  fasp_blas_smat_mul_nc7(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[9][ixyxc]),smat);
987  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[3][ic]));
988  }
989 
990  }
991 
992  //comput beta3[i]
993  if (i+nline<ngrid) {
994 
995  memcpy(&(LU->offdiag[5][ic]),&(A->offdiag[3][ic]), 49*sizeof(REAL));
996 
997  if (ixyx>=0) {
998  fasp_blas_smat_mul_nc7(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[11][ixyxc]),smat);
999  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[5][ic]));
1000  }
1001 
1002  }
1003 
1004  //comput beta4[i]
1005  if (i+nplane-nline<ngrid) {
1006 
1007  if (ix1>=0) {
1008  fasp_blas_smat_mul_nc7(&(LU->offdiag[2][ix1c]),&(LU->offdiag[9][ix1c]),smat);
1009  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[7][ic]));
1010  }
1011 
1012  if (ix>=0) {
1013  fasp_blas_smat_mul_nc7(&(LU->offdiag[4][ixc]),&(LU->offdiag[11][ixc]),smat);
1014  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[7][ic]));
1015  }
1016 
1017  }
1018 
1019  //comput beta5[i]
1020  if (i+nplane-1<ngrid) {
1021  fasp_blas_smat_mul_nc7(&(LU->offdiag[0][i1c]),&(LU->offdiag[11][i1c]),smat);
1022  fasp_blas_array_axpy_nc7(-1,smat,&(LU->offdiag[9][ic]));
1023  }
1024 
1025  //comput d[i]
1026  {
1027  fasp_blas_smat_mul_nc7(&(LU->offdiag[0][i1c]),&(LU->offdiag[1][i1c]),smat);
1028  fasp_blas_array_axpyz_nc7(-1,smat,&(A->diag[ic]),&(LU->diag[ic]));
1029  }
1030 
1031  if (ix1>=0) {
1032  fasp_blas_smat_mul_nc7(&(LU->offdiag[2][ix1c]),&(LU->offdiag[3][ix1c]),smat);
1033  fasp_blas_array_axpy_nc7(-1,smat,&(LU->diag[ic]));
1034  }
1035 
1036  if (ix>=0) {
1037  fasp_blas_smat_mul_nc7(&(LU->offdiag[4][ixc]),&(LU->offdiag[5][ixc]),smat);
1038  fasp_blas_array_axpy_nc7(-1,smat,&(LU->diag[ic]));
1039  }
1040 
1041  if (ixyx>=0) {
1042  fasp_blas_smat_mul_nc7(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[7][ixyxc]),smat);
1043  fasp_blas_array_axpy_nc7(-1,smat,&(LU->diag[ic]));
1044  }
1045 
1046  if (ixy1>=0) {
1047  fasp_blas_smat_mul_nc7(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[9][ixy1c]),smat);
1048  fasp_blas_array_axpy_nc7(-1,smat,&(LU->diag[ic]));
1049  }
1050 
1051  if (ixy>=0) {
1052  fasp_blas_smat_mul_nc7(&(LU->offdiag[10][ixyc]),&(LU->offdiag[11][ixyc]),smat);
1053  fasp_blas_array_axpy_nc7(-1,smat,&(LU->diag[ic]));
1054  }
1055 
1056  //fasp_blas_smat_inv_nc5(&(LU->diag[ic]));
1057  fasp_blas_smat_inv(&(LU->diag[ic]), 7);
1058 
1059  } // end for(i=1;i<ngrid;++i)
1060 
1061  } // end if (nc == 7)
1062 
1063  else {
1064  //comput the first row
1065  fasp_blas_smat_inv(LU->diag,nc);
1066  memcpy(LU->offdiag[1],A->offdiag[1],nc2*sizeof(REAL));
1067  memcpy(LU->offdiag[5],A->offdiag[3],nc2*sizeof(REAL));
1068 
1069  for(i=1;i<ngrid;++i) {
1070 
1071  i1=i-1;ix=i-nline;ixy=i-nplane;ix1=ix+1;ixyx=ixy+nline;ixy1=ixy+1;
1072  ic=i*nc2;i1c=i1*nc2;ixc=ix*nc2;ix1c=ix1*nc2;ixyc=ixy*nc2;ixy1c=ixy1*nc2;ixyxc=ixyx*nc2;
1073  //comput alpha6[i-nxy]
1074  if (ixy>=0)
1075  fasp_blas_smat_mul(&(A->offdiag[4][ixyc]),&(LU->diag[ixyc]),&(LU->offdiag[10][ixyc]),nc);
1076 
1077  //comput alpha5[ixy1]
1078  if (ixy1>=0) {
1079  for (j=0;j<nc2;++j) tc[j]=0;
1080  if (ixy>=0) {
1081  fasp_blas_smat_mul(&(LU->offdiag[10][ixyc]),&(LU->offdiag[1][ixyc]),smat,nc);
1082  fasp_blas_array_axpy(nc2,-1,smat,tc);
1083  }
1084 
1085  fasp_blas_smat_mul(tc,&(LU->diag[ixy1c]),&(LU->offdiag[8][ixy1c]),nc);
1086  }
1087 
1088  //comput alpha4[ixyx]
1089  if (ixyx>=0) {
1090  for (j=0;j<nc2;++j) tc[j]=0;
1091  if (ixy>=0) {
1092  fasp_blas_smat_mul(&(LU->offdiag[10][ixyc]),&(LU->offdiag[5][ixyc]),smat,nc);
1093  fasp_blas_array_axpy(nc2,-1,smat,tc);
1094  }
1095  if (ixy1>=0) {
1096  fasp_blas_smat_mul(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[3][ixy1c]),smat,nc);
1097  fasp_blas_array_axpy(nc2,-1,smat,tc);
1098  }
1099 
1100  fasp_blas_smat_mul(tc,&(LU->diag[ixyxc]),&(LU->offdiag[6][ixyxc]),nc);
1101  }
1102 
1103  //comput alpha3[ix]
1104  if (ix>=0) {
1105 
1106  memcpy(tc,&(A->offdiag[2][ixc]),nc2*sizeof(REAL));
1107  if (ixy>=0) {
1108  fasp_blas_smat_mul(&(LU->offdiag[10][ixyc]),&(LU->offdiag[7][ixyc]),smat,nc);
1109  fasp_blas_array_axpy(nc2,-1,smat,tc);
1110  }
1111 
1112  fasp_blas_smat_mul(tc,&(LU->diag[ixc]),&(LU->offdiag[4][ixc]),nc);
1113  }
1114 
1115  //comput alpha2[i-nx+1]
1116  if (ix1>=0) {
1117 
1118  for (j=0;j<nc2;++j) tc[j]=0;
1119 
1120  if (ix>=0) {
1121  fasp_blas_smat_mul(&(LU->offdiag[4][ixc]),&(LU->offdiag[1][ixc]),smat,nc);
1122  fasp_blas_array_axpy(nc2,-1,smat,tc);
1123  }
1124 
1125  if (ixy1>=0) {
1126  fasp_blas_smat_mul(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[7][ixy1c]),smat,nc);
1127  fasp_blas_array_axpy(nc2,-1,smat,tc);
1128  }
1129 
1130  fasp_blas_smat_mul(tc,&(LU->diag[ix1c]),&(LU->offdiag[2][ix1c]),nc);
1131  }
1132 
1133  //comput alpha1[i-1]
1134 
1135  memcpy(tc,&(A->offdiag[0][i1c]),nc2*sizeof(REAL));
1136  if (ix>=0) {
1137  fasp_blas_smat_mul(&(LU->offdiag[4][ixc]),&(LU->offdiag[3][ixc]),smat,nc);
1138  fasp_blas_array_axpy(nc2,-1,smat,tc);
1139  }
1140  if (ixy>=0) {
1141  fasp_blas_smat_mul(&(LU->offdiag[10][ixyc]),&(LU->offdiag[9][ixyc]),smat,nc);
1142  fasp_blas_array_axpy(nc2,-1,smat,tc);
1143  }
1144 
1145  fasp_blas_smat_mul(tc,&(LU->diag[i1c]),&(LU->offdiag[0][i1c]),nc);
1146 
1147  //comput beta1[i]
1148  if (i+1<ngrid) {
1149 
1150  memcpy(&(LU->offdiag[1][ic]),&(A->offdiag[1][ic]),nc2*sizeof(REAL));
1151  if (ix1>=0) {
1152  fasp_blas_smat_mul(&(LU->offdiag[2][ix1c]),&(LU->offdiag[5][ix1c]),smat,nc);
1153  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[1][ic]));
1154  }
1155  if (ixy1>=0) {
1156  fasp_blas_smat_mul(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[11][ixy1c]),smat,nc);
1157  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[1][ic]));
1158  }
1159 
1160  }
1161 
1162  //comput beta2[i]
1163  if (i+nline-1<ngrid) {
1164 
1165  {
1166  fasp_blas_smat_mul(&(LU->offdiag[0][i1c]),&(LU->offdiag[5][i1c]),smat,nc);
1167  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[3][ic]));
1168  }
1169 
1170  if (ixyx>=0) {
1171  fasp_blas_smat_mul(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[9][ixyxc]),smat,nc);
1172  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[3][ic]));
1173  }
1174 
1175  }
1176 
1177  //comput beta3[i]
1178  if (i+nline<ngrid) {
1179 
1180  memcpy(&(LU->offdiag[5][ic]),&(A->offdiag[3][ic]),nc2*sizeof(REAL));
1181  if (ixyx>=0) {
1182  fasp_blas_smat_mul(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[11][ixyxc]),smat,nc);
1183  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[5][ic]));
1184  }
1185 
1186  }
1187 
1188  //comput beta4[i]
1189  if (i+nplane-nline<ngrid) {
1190 
1191  if (ix1>=0) {
1192  fasp_blas_smat_mul(&(LU->offdiag[2][ix1c]),&(LU->offdiag[9][ix1c]),smat,nc);
1193  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[7][ic]));
1194  }
1195 
1196  if (ix>=0) {
1197  fasp_blas_smat_mul(&(LU->offdiag[4][ixc]),&(LU->offdiag[11][ixc]),smat,nc);
1198  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[7][ic]));
1199  }
1200  }
1201 
1202  //comput beta5[i]
1203  if (i+nplane-1<ngrid) {
1204  fasp_blas_smat_mul(&(LU->offdiag[0][i1c]),&(LU->offdiag[11][i1c]),smat,nc);
1205  fasp_blas_array_axpy(nc2,-1,smat,&(LU->offdiag[9][ic]));
1206  }
1207 
1208  //comput d[i]
1209  {
1210  fasp_blas_smat_mul(&(LU->offdiag[0][i1c]),&(LU->offdiag[1][i1c]),smat,nc);
1211  fasp_blas_array_axpyz(nc2,-1,smat,&(A->diag[ic]),&(LU->diag[ic]));
1212  }
1213 
1214  if (ix1>=0) {
1215  fasp_blas_smat_mul(&(LU->offdiag[2][ix1c]),&(LU->offdiag[3][ix1c]),smat,nc);
1216  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ic]));
1217  }
1218 
1219  if (ix>=0) {
1220  fasp_blas_smat_mul(&(LU->offdiag[4][ixc]),&(LU->offdiag[5][ixc]),smat,nc);
1221  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ic]));
1222  }
1223 
1224  if (ixyx>=0) {
1225  fasp_blas_smat_mul(&(LU->offdiag[6][ixyxc]),&(LU->offdiag[7][ixyxc]),smat,nc);
1226  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ic]));
1227  }
1228 
1229 
1230  if (ixy1>=0) {
1231  fasp_blas_smat_mul(&(LU->offdiag[8][ixy1c]),&(LU->offdiag[9][ixy1c]),smat,nc);
1232  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ic]));
1233  }
1234 
1235  if (ixy>=0) {
1236  fasp_blas_smat_mul(&(LU->offdiag[10][ixyc]),&(LU->offdiag[11][ixyc]),smat,nc);
1237  fasp_blas_array_axpy(nc2,-1,smat,&(LU->diag[ic]));
1238  }
1239 
1240  fasp_blas_smat_inv(&(LU->diag[ic]),nc);
1241 
1242  }
1243 
1244  } // end else
1245 
1246  fasp_mem_free(smat);
1247  fasp_mem_free(tc);
1248 
1249  return;
1250 }
1251 
1252 /*---------------------------------*/
1253 /*-- End of File --*/
1254 /*---------------------------------*/
INT nx
number of grids in x direction
Definition: fasp.h:307
#define REAL
Definition: fasp.h:67
void fasp_blas_array_axpyz_nc3(const REAL a, REAL *x, REAL *y, REAL *z)
z = a*x + y
Definition: blas_smat.c:527
INT ngrid
number of grids
Definition: fasp.h:322
INT fasp_blas_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: smat.c:554
INT nc
size of each block (number of components)
Definition: fasp.h:319
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
void fasp_blas_smat_mul_nc5(REAL *a, REAL *b, REAL *c)
Compute the matrix product of two 5*5 matrices a and b, stored in c.
Definition: blas_smat.c:299
INT * offsets
offsets of the off-diagonals (length is nband)
Definition: fasp.h:331
Structure matrix of REAL type.
Definition: fasp.h:304
INT nz
number of grids in z direction
Definition: fasp.h:313
void fasp_ilu_dstr_setup1(dSTRmat *A, dSTRmat *LU)
Get ILU(1) decoposition of a structured matrix A.
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void fasp_blas_smat_mul(REAL *a, REAL *b, REAL *c, const INT n)
Compute the matrix product of two small full matrices a and b, stored in c.
Definition: blas_smat.c:448
#define INT
Definition: fasp.h:64
INT nband
number of off-diag bands
Definition: fasp.h:328
void fasp_blas_smat_mul_nc3(REAL *a, REAL *b, REAL *c)
Compute the matrix product of two 3*3 matrices a and b, stored in c.
Definition: blas_smat.c:262
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
Main header file for FASP.
void fasp_ilu_dstr_setup0(dSTRmat *A, dSTRmat *LU)
Get ILU(0) decomposition of a structured matrix A.
Definition: ilu_setup_str.c:28
void fasp_blas_smat_inv_nc7(REAL *a)
Compute the inverse matrix of a 7*7 matrix a.
Definition: smat.c:389
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_smat_inv_nc5(REAL *a)
Compute the inverse matrix of a 5*5 full matrix A (in place)
Definition: smat.c:173
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
REAL * diag
diagonal entries (length is ngrid*(nc^2))
Definition: fasp.h:325
void fasp_blas_smat_inv_nc3(REAL *a)
Compute the inverse matrix of a 3*3 full matrix A (in place)
Definition: smat.c:61
void fasp_blas_smat_mul_nc7(REAL *a, REAL *b, REAL *c)
Compute the matrix product of two 7*7 matrices a and b, stored in c.
Definition: blas_smat.c:358
#define ABS(a)
Definition: fasp.h:74
void fasp_dstr_alloc(const INT nx, const INT ny, const INT nz, const INT nxy, const INT ngrid, const INT nband, const INT nc, INT *offsets, dSTRmat *A)
Allocate STR sparse matrix memory space.
Definition: sparse_str.c:109
INT nxy
number of grids on x-y plane
Definition: fasp.h:316
void fasp_blas_array_axpyz_nc5(const REAL a, REAL *x, REAL *y, REAL *z)
z = a*x + y
Definition: blas_smat.c:560
void fasp_blas_array_axpyz_nc7(const REAL a, REAL *x, REAL *y, REAL *z)
z = a*x + y
Definition: blas_smat.c:611
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
void fasp_blas_array_axpyz(const INT n, const REAL a, REAL *x, REAL *y, REAL *z)
z = a*x + y
Definition: blas_array.c:167
INT ny
number of grids in y direction
Definition: fasp.h:310