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