Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
smoother_bsr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 
15 /*---------------------------------*/
16 /*-- Public Functions --*/
17 /*---------------------------------*/
18 
34  dvector *b,
35  dvector *u)
36 {
37  // members of A
38  const INT ROW = A->ROW;
39  const INT nb = A->nb;
40  const INT nb2 = nb*nb;
41  const INT size = ROW*nb2;
42  const INT *IA = A->IA;
43  const INT *JA = A->JA;
44  const REAL *val = A->val;
45 
46  // local variables
47  INT i,k;
48  REAL *diaginv = NULL;
49 
50  INT nthreads = 1, use_openmp = FALSE;
51 
52 #ifdef _OPENMP
53  if ( ROW > OPENMP_HOLDS ) {
54  use_openmp = TRUE;
55  nthreads = FASP_GET_NUM_THREADS();
56  }
57 #endif
58 
59  // allocate memory
60  diaginv = (REAL *)fasp_mem_calloc(size, sizeof(REAL));
61 
62  // get all the diagonal sub-blocks
63  if(use_openmp) {
64  INT mybegin,myend,myid;
65 #ifdef _OPENMP
66 #pragma omp parallel for private(myid,mybegin,myend,i,k)
67 #endif
68  for(myid=0; myid<nthreads; myid++) {
69  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
70  for(i=mybegin; i<myend; i++) {
71  for(k=IA[i]; k<IA[i+1]; ++k)
72  if(JA[k] == i)
73  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
74  }
75  }
76  }
77  else {
78  for (i = 0; i < ROW; ++i) {
79  for (k = IA[i]; k < IA[i+1]; ++k) {
80  if (JA[k] == i)
81  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
82  }
83  }
84  }
85 
86  // compute the inverses of all the diagonal sub-blocks
87  if (nb > 1) {
88  if(use_openmp) {
89  INT mybegin,myend,myid;
90 #ifdef _OPENMP
91 #pragma omp parallel for private(myid,mybegin,myend,i)
92 #endif
93  for(myid=0; myid<nthreads; myid++) {
94  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
95  for(i=mybegin; i<myend; i++) {
96  fasp_blas_smat_inv(diaginv+i*nb2, nb);
97  }
98  }
99  }
100  else {
101  for (i = 0; i < ROW; ++i) {
102  fasp_blas_smat_inv(diaginv+i*nb2, nb);
103  }
104  }
105  }
106  else {
107  if(use_openmp) {
108  INT mybegin, myend, myid;
109 #ifdef _OPENMP
110 #pragma omp parallel for private(myid,mybegin,myend,i)
111 #endif
112  for(myid=0; myid<nthreads; myid++) {
113  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
114  for(i=mybegin; i<myend; i++) {
115  diaginv[i] = 1.0 / diaginv[i];
116  }
117  }
118  }
119  else {
120  for (i = 0; i < ROW; ++i) {
121  // zero-diagonal should be tested previously
122  diaginv[i] = 1.0 / diaginv[i];
123  }
124  }
125  }
126 
127  fasp_smoother_dbsr_jacobi1(A, b, u, diaginv);
128  fasp_mem_free(diaginv);
129 }
130 
149  dvector *b,
150  dvector *u,
151  REAL *diaginv)
152 {
153  // members of A
154  const INT ROW = A->ROW;
155  const INT nb = A->nb;
156  const INT nb2 = nb*nb;
157  // const INT size = ROW*nb2;
158  const INT *IA = A->IA;
159  const INT *JA = A->JA;
160  const REAL *val = A->val;
161 
162  // local variables
163  INT i,k;
164 
165  INT nthreads = 1, use_openmp = FALSE;
166 
167 #ifdef _OPENMP
168  if ( ROW > OPENMP_HOLDS ) {
169  use_openmp = TRUE;
170  nthreads = FASP_GET_NUM_THREADS();
171  }
172 #endif
173 
174  // get all the diagonal sub-blocks
175  if(use_openmp) {
176  INT mybegin,myend,myid;
177 #ifdef _OPENMP
178 #pragma omp parallel for private(myid,mybegin,myend,i,k)
179 #endif
180  for(myid=0; myid<nthreads; myid++) {
181  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
182  for(i=mybegin; i<myend; i++) {
183  for(k=IA[i]; k<IA[i+1]; ++k)
184  if(JA[k] == i)
185  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
186  }
187  }
188  }
189  else {
190  for (i = 0; i < ROW; ++i) {
191  for (k = IA[i]; k < IA[i+1]; ++k) {
192  if (JA[k] == i)
193  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
194  }
195  }
196  }
197 
198  // compute the inverses of all the diagonal sub-blocks
199  if (nb > 1) {
200  if(use_openmp) {
201  INT mybegin,myend,myid;
202 #ifdef _OPENMP
203 #pragma omp parallel for private(myid,mybegin,myend,i)
204 #endif
205  for(myid=0; myid<nthreads; myid++) {
206  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
207  for(i=mybegin; i<myend; i++) {
208  fasp_blas_smat_inv(diaginv+i*nb2, nb);
209  }
210  }
211  }
212  else {
213  for (i = 0; i < ROW; ++i) {
214  fasp_blas_smat_inv(diaginv+i*nb2, nb);
215  }
216  }
217  }
218  else {
219  if(use_openmp) {
220  INT mybegin, myend, myid;
221 #ifdef _OPENMP
222 #pragma omp parallel for private(myid,mybegin,myend,i)
223 #endif
224  for(myid=0; myid<nthreads; myid++) {
225  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
226  for(i=mybegin; i<myend; i++) {
227  diaginv[i] = 1.0 / diaginv[i];
228  }
229  }
230  }
231  else {
232  for (i = 0; i < ROW; ++i) {
233  // zero-diagonal should be tested previously
234  diaginv[i] = 1.0 / diaginv[i];
235  }
236  }
237  }
238 
239 }
240 
258  dvector *b,
259  dvector *u,
260  REAL *diaginv)
261 {
262  // members of A
263  const INT ROW = A->ROW;
264  const INT nb = A->nb;
265  const INT nb2 = nb*nb;
266  const INT size = ROW*nb;
267  const INT *IA = A->IA;
268  const INT *JA = A->JA;
269  REAL *val = A->val;
270 
271  INT nthreads = 1, use_openmp = FALSE;
272 
273 #ifdef _OPENMP
274  if ( ROW > OPENMP_HOLDS ) {
275  use_openmp = TRUE;
276  nthreads = FASP_GET_NUM_THREADS();
277  }
278 #endif
279 
280  // values of dvector b and u
281  REAL *b_val = b->val;
282  REAL *u_val = u->val;
283 
284  // auxiliary array
285  REAL *b_tmp = NULL;
286 
287  // local variables
288  INT i,j,k;
289  INT pb;
290 
291  // b_tmp = b_val
292  b_tmp = (REAL *)fasp_mem_calloc(size, sizeof(REAL));
293  memcpy(b_tmp, b_val, size*sizeof(REAL));
294 
295  // It's not necessary to assign the smoothing order since the result doesn't depend on it
296  if (nb == 1) {
297  if(use_openmp) {
298  INT mybegin, myend, myid;
299 #ifdef _OPENMP
300 #pragma omp parallel for private(myid,mybegin,myend,i,j,k)
301 #endif
302  for (myid=0; myid<nthreads; myid++) {
303  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
304  for (i=mybegin; i<myend; i++) {
305  for (k = IA[i]; k < IA[i+1]; ++k) {
306  j = JA[k];
307  if (j != i)
308  b_tmp[i] -= val[k]*u_val[j];
309  }
310  }
311  }
312 #ifdef _OPENMP
313 #pragma omp parallel for private(myid,mybegin,myend,i)
314 #endif
315  for(myid=0; myid<nthreads; myid++) {
316  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
317  for(i=mybegin; i<myend; i++) {
318  u_val[i] = b_tmp[i]*diaginv[i];
319  }
320  }
321  }
322  else {
323  for (i = 0; i < ROW; ++i) {
324  for (k = IA[i]; k < IA[i+1]; ++k) {
325  j = JA[k];
326  if (j != i)
327  b_tmp[i] -= val[k]*u_val[j];
328  }
329  }
330  for (i = 0; i < ROW; ++i) {
331  u_val[i] = b_tmp[i]*diaginv[i];
332  }
333  }
334 
335  fasp_mem_free(b_tmp);
336  }
337  else if (nb > 1) {
338  if(use_openmp) {
339  INT mybegin, myend, myid;
340 #ifdef _OPENMP
341 #pragma omp parallel for private(myid,mybegin,myend,i,pb,k,j)
342 #endif
343  for(myid=0; myid<nthreads; myid++) {
344  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
345  for (i=mybegin; i<myend; i++) {
346  pb = i*nb;
347  for (k = IA[i]; k < IA[i+1]; ++k) {
348  j = JA[k];
349  if (j != i)
350  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp+pb, nb);
351  }
352  }
353  }
354 #ifdef _OPENMP
355 #pragma omp parallel for private(myid,mybegin,myend,i,pb)
356 #endif
357  for(myid=0; myid<nthreads; myid++) {
358  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
359  for (i=mybegin; i<myend; i++) {
360  pb = i*nb;
361  fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp+pb, u_val+pb, nb);
362  }
363  }
364  }
365  else {
366  for (i = 0; i < ROW; ++i) {
367  pb = i*nb;
368  for (k = IA[i]; k < IA[i+1]; ++k) {
369  j = JA[k];
370  if (j != i)
371  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp+pb, nb);
372  }
373  }
374 
375  for (i = 0; i < ROW; ++i) {
376  pb = i*nb;
377  fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp+pb, u_val+pb, nb);
378  }
379 
380  }
381  fasp_mem_free(b_tmp);
382  }
383  else {
384  printf("### ERROR: nb is illegal! %s : %d\n", __FILE__, __LINE__);
385  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
386  }
387 
388 }
389 
412  dvector *b,
413  dvector *u,
414  INT order,
415  INT *mark )
416 {
417  // members of A
418  const INT ROW = A->ROW;
419  const INT nb = A->nb;
420  const INT nb2 = nb*nb;
421  const INT size = ROW*nb2;
422  const INT *IA = A->IA;
423  const INT *JA = A->JA;
424  const REAL *val = A->val;
425 
426  // local variables
427  INT i,k;
428  REAL *diaginv = NULL;
429 
430  INT nthreads = 1, use_openmp = FALSE;
431 
432 #ifdef _OPENMP
433  if ( ROW > OPENMP_HOLDS ) {
434  use_openmp = TRUE;
435  nthreads = FASP_GET_NUM_THREADS();
436  }
437 #endif
438 
439  // allocate memory
440  diaginv = (REAL *)fasp_mem_calloc(size, sizeof(REAL));
441 
442  // get all the diagonal sub-blocks
443  if(use_openmp) {
444  INT mybegin,myend,myid;
445 #ifdef _OPENMP
446 #pragma omp parallel for private(myid,mybegin,myend,i,k)
447 #endif
448  for(myid=0; myid<nthreads; myid++) {
449  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
450  for(i=mybegin; i<myend; i++) {
451  for(k=IA[i]; k<IA[i+1]; ++k)
452  if(JA[k] == i)
453  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
454  }
455  }
456  }
457  else {
458  for (i = 0; i < ROW; ++i) {
459  for (k = IA[i]; k < IA[i+1]; ++k) {
460  if (JA[k] == i)
461  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
462  }
463  }
464  }
465 
466  // compute the inverses of all the diagonal sub-blocks
467  if (nb > 1) {
468  if(use_openmp) {
469  INT mybegin,myend,myid;
470 #ifdef _OPENMP
471 #pragma omp parallel for private(myid,mybegin,myend,i)
472 #endif
473  for(myid=0; myid<nthreads; myid++) {
474  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
475  for(i=mybegin; i<myend; i++) {
476  fasp_blas_smat_inv(diaginv+i*nb2, nb);
477  }
478  }
479  }
480  else {
481  for (i = 0; i < ROW; ++i) {
482  fasp_blas_smat_inv(diaginv+i*nb2, nb);
483  }
484  }
485  }
486  else {
487  if(use_openmp) {
488  INT mybegin, myend, myid;
489 #ifdef _OPENMP
490 #pragma omp parallel for private(myid,mybegin,myend,i)
491 #endif
492  for(myid=0; myid<nthreads; myid++) {
493  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
494  for(i=mybegin; i<myend; i++) {
495  diaginv[i] = 1.0 / diaginv[i];
496  }
497  }
498  }
499  else {
500  for (i = 0; i < ROW; ++i) {
501  // zero-diagonal should be tested previously
502  diaginv[i] = 1.0 / diaginv[i];
503  }
504  }
505  }
506 
507  fasp_smoother_dbsr_gs1(A, b, u, order, mark, diaginv);
508  fasp_mem_free(diaginv);
509 }
510 
532  dvector *b,
533  dvector *u,
534  INT order,
535  INT *mark,
536  REAL *diaginv)
537 {
538  if (!mark) {
539  if (order == ASCEND) // smooth ascendingly
540  {
541  fasp_smoother_dbsr_gs_ascend(A, b, u, diaginv);
542  }
543  else if (order == DESCEND) // smooth descendingly
544  {
545  fasp_smoother_dbsr_gs_descend(A, b, u, diaginv);
546  }
547  }
548  // smooth according to the order 'mark' defined by user
549  else {
550  fasp_smoother_dbsr_gs_order1(A, b, u, diaginv, mark);
551  }
552 }
553 
569  dvector *b,
570  dvector *u,
571  REAL *diaginv)
572 {
573  // members of A
574  const INT ROW = A->ROW;
575  const INT nb = A->nb;
576  const INT nb2 = nb*nb;
577  const INT *IA = A->IA;
578  const INT *JA = A->JA;
579  REAL *val = A->val;
580 
581  // values of dvector b and u
582  REAL *b_val = b->val;
583  REAL *u_val = u->val;
584 
585  // local variables
586  INT i,j,k;
587  INT pb;
588  REAL rhs = 0.0;
589 
590  if (nb == 1) {
591  for (i = 0; i < ROW; ++i) {
592  rhs = b_val[i];
593  for (k = IA[i]; k < IA[i+1]; ++k) {
594  j = JA[k];
595  if (j != i)
596  rhs -= val[k]*u_val[j];
597  }
598  u_val[i] = rhs*diaginv[i];
599  }
600  }
601  else if (nb > 1) {
602  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
603 
604  for (i = 0; i < ROW; ++i) {
605  pb = i*nb;
606  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
607  for (k = IA[i]; k < IA[i+1]; ++k) {
608  j = JA[k];
609  if (j != i)
610  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
611  }
612  fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp, u_val+pb, nb);
613  }
614 
615  fasp_mem_free(b_tmp);
616  }
617  else {
618  printf("### ERROR: nb is illegal! %s : %d\n", __FILE__, __LINE__);
619  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
620  }
621 
622 }
623 
642  dvector *b,
643  dvector *u)
644 {
645  // members of A
646  const INT ROW = A->ROW;
647  const INT nb = A->nb;
648  const INT nb2 = nb*nb;
649  const INT *IA = A->IA;
650  const INT *JA = A->JA;
651  REAL *val = A->val;
652 
653  // values of dvector b and u
654  REAL *b_val = b->val;
655  REAL *u_val = u->val;
656 
657  // local variables
658  INT i,j,k;
659  INT pb;
660  REAL rhs = 0.0;
661 
662  if (nb == 1) {
663  for (i = 0; i < ROW; ++i) {
664  rhs = b_val[i];
665  for (k = IA[i]; k < IA[i+1]; ++k) {
666  j = JA[k];
667  if (j != i)
668  rhs -= val[k]*u_val[j];
669  }
670  //u_val[i] = rhs*diaginv[i];
671  u_val[i] = rhs;
672  }
673  }
674  else if (nb > 1) {
675  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
676 
677  for (i = 0; i < ROW; ++i) {
678  pb = i*nb;
679  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
680  for (k = IA[i]; k < IA[i+1]; ++k) {
681  j = JA[k];
682  if (j != i)
683  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
684  }
685  //fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp, u_val+pb, nb);
686  memcpy(u_val+pb, b_tmp, nb*sizeof(REAL));
687  }
688 
689  fasp_mem_free(b_tmp);
690  }
691  else {
692  printf("### ERROR: nb is illegal! %s : %d\n", __FILE__, __LINE__);
693  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
694  }
695 
696 }
697 
713  dvector *b,
714  dvector *u,
715  REAL *diaginv )
716 {
717  // members of A
718  const INT ROW = A->ROW;
719  const INT nb = A->nb;
720  const INT nb2 = nb*nb;
721  const INT *IA = A->IA;
722  const INT *JA = A->JA;
723  REAL *val = A->val;
724 
725  // values of dvector b and u
726  REAL *b_val = b->val;
727  REAL *u_val = u->val;
728 
729  // local variables
730  INT i,j,k;
731  INT pb;
732  REAL rhs = 0.0;
733 
734  if (nb == 1) {
735  for (i = ROW-1; i >= 0; i--) {
736  rhs = b_val[i];
737  for (k = IA[i]; k < IA[i+1]; ++k) {
738  j = JA[k];
739  if (j != i)
740  rhs -= val[k]*u_val[j];
741  }
742  u_val[i] = rhs*diaginv[i];
743  }
744  }
745  else if (nb > 1) {
746  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
747 
748  for (i = ROW-1; i >= 0; i--) {
749  pb = i*nb;
750  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
751  for (k = IA[i]; k < IA[i+1]; ++k) {
752  j = JA[k];
753  if (j != i)
754  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
755  }
756  fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp, u_val+pb, nb);
757  }
758 
759  fasp_mem_free(b_tmp);
760  }
761  else {
762  printf("### ERROR: nb is illegal! %s : %d\n", __FILE__, __LINE__);
763  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
764  }
765 
766 }
767 
787  dvector *b,
788  dvector *u)
789 {
790  // members of A
791  const INT ROW = A->ROW;
792  const INT nb = A->nb;
793  const INT nb2 = nb*nb;
794  const INT *IA = A->IA;
795  const INT *JA = A->JA;
796  REAL *val = A->val;
797 
798  // values of dvector b and u
799  REAL *b_val = b->val;
800  REAL *u_val = u->val;
801 
802  // local variables
803  INT i,j,k;
804  INT pb;
805  REAL rhs = 0.0;
806 
807  if (nb == 1) {
808  for (i = ROW-1; i >= 0; i--) {
809  rhs = b_val[i];
810  for (k = IA[i]; k < IA[i+1]; ++k) {
811  j = JA[k];
812  if (j != i)
813  rhs -= val[k]*u_val[j];
814  }
815  //u_val[i] = rhs*diaginv[i];
816  u_val[i] = rhs;
817  }
818  }
819  else if (nb > 1) {
820  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
821 
822  for (i = ROW-1; i >= 0; i--) {
823  pb = i*nb;
824  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
825  for (k = IA[i]; k < IA[i+1]; ++k) {
826  j = JA[k];
827  if (j != i)
828  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
829  }
830  //fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp, u_val+pb, nb);
831  memcpy(u_val+pb, b_tmp, nb*sizeof(REAL));
832  }
833 
834  fasp_mem_free(b_tmp);
835  }
836  else {
837  printf("### ERROR: nb is illegal! %s : %d\n", __FILE__, __LINE__);
838  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
839  }
840 
841 }
842 
859  dvector *b,
860  dvector *u,
861  REAL *diaginv,
862  INT *mark)
863 {
864  // members of A
865  const INT ROW = A->ROW;
866  const INT nb = A->nb;
867  const INT nb2 = nb*nb;
868  const INT *IA = A->IA;
869  const INT *JA = A->JA;
870  REAL *val = A->val;
871 
872  // values of dvector b and u
873  REAL *b_val = b->val;
874  REAL *u_val = u->val;
875 
876  // local variables
877  INT i,j,k;
878  INT I,pb;
879  REAL rhs = 0.0;
880 
881  if (nb == 1) {
882  for (I = 0; I < ROW; ++I) {
883  i = mark[I];
884  rhs = b_val[i];
885  for (k = IA[i]; k < IA[i+1]; ++k) {
886  j = JA[k];
887  if (j != i)
888  rhs -= val[k]*u_val[j];
889  }
890  u_val[i] = rhs*diaginv[i];
891  }
892  }
893  else if (nb > 1) {
894  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
895 
896  for (I = 0; I < ROW; ++I) {
897  i = mark[I];
898  pb = i*nb;
899  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
900  for (k = IA[i]; k < IA[i+1]; ++k) {
901  j = JA[k];
902  if (j != i)
903  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
904  }
905  fasp_blas_smat_mxv(diaginv+nb2*i, b_tmp, u_val+pb, nb);
906  }
907 
908  fasp_mem_free(b_tmp);
909  }
910  else {
911  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
912  }
913 
914 }
915 
937  dvector *b,
938  dvector *u,
939  INT *mark,
940  REAL *work)
941 {
942  // members of A
943  const INT ROW = A->ROW;
944  const INT nb = A->nb;
945  const INT nb2 = nb*nb;
946  const INT *IA = A->IA;
947  const INT *JA = A->JA;
948  REAL *val = A->val;
949 
950  // values of dvector b and u
951  REAL *b_val = b->val;
952  REAL *u_val = u->val;
953 
954  // auxiliary array
955  REAL *b_tmp = work;
956 
957  // local variables
958  INT i,j,k,I,pb;
959  REAL rhs = 0.0;
960 
961  if (nb == 1) {
962  for (I = 0; I < ROW; ++I) {
963  i = mark[I];
964  rhs = b_val[i];
965  for (k = IA[i]; k < IA[i+1]; ++k) {
966  j = JA[k];
967  if (j != i)
968  rhs -= val[k]*u_val[j];
969  }
970  u_val[i] = rhs;
971  }
972  }
973  else if (nb > 1) {
974  for (I = 0; I < ROW; ++I) {
975  i = mark[I];
976  pb = i*nb;
977  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
978  for (k = IA[i]; k < IA[i+1]; ++k) {
979  j = JA[k];
980  if (j != i)
981  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
982  }
983  memcpy(u_val+pb, b_tmp, nb*sizeof(REAL));
984  }
985  }
986  else {
987  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
988  }
989 }
990 
1014  dvector *b,
1015  dvector *u,
1016  INT order,
1017  INT *mark,
1018  REAL weight)
1019 {
1020  // members of A
1021  const INT ROW = A->ROW;
1022  const INT nb = A->nb;
1023  const INT nb2 = nb*nb;
1024  const INT size = ROW*nb2;
1025  const INT *IA = A->IA;
1026  const INT *JA = A->JA;
1027  const REAL *val = A->val;
1028 
1029  // local variables
1030  INT i,k;
1031  REAL *diaginv = NULL;
1032 
1033  INT nthreads = 1, use_openmp = FALSE;
1034 
1035 #ifdef _OPENMP
1036  if ( ROW > OPENMP_HOLDS ) {
1037  use_openmp = TRUE;
1038  nthreads = FASP_GET_NUM_THREADS();
1039  }
1040 #endif
1041 
1042  // allocate memory
1043  diaginv = (REAL *)fasp_mem_calloc(size, sizeof(REAL));
1044 
1045  // get all the diagonal sub-blocks
1046  if(use_openmp) {
1047  INT mybegin,myend,myid;
1048 #ifdef _OPENMP
1049 #pragma omp parallel for private(myid,mybegin,myend,i,k)
1050 #endif
1051  for(myid=0; myid<nthreads; myid++) {
1052  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1053  for(i=mybegin; i<myend; i++) {
1054  for(k=IA[i]; k<IA[i+1]; ++k)
1055  if(JA[k] == i)
1056  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
1057  }
1058  }
1059  }
1060  else {
1061  for (i = 0; i < ROW; ++i) {
1062  for (k = IA[i]; k < IA[i+1]; ++k) {
1063  if (JA[k] == i)
1064  memcpy(diaginv+i*nb2, val+k*nb2, nb2*sizeof(REAL));
1065  }
1066  }
1067  }
1068 
1069  // compute the inverses of all the diagonal sub-blocks
1070  if (nb > 1) {
1071  if(use_openmp) {
1072  INT mybegin,myend,myid;
1073 #ifdef _OPENMP
1074 #pragma omp parallel for private(myid,mybegin,myend,i)
1075 #endif
1076  for(myid=0; myid<nthreads; myid++) {
1077  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1078  for(i=mybegin; i<myend; i++) {
1079  fasp_blas_smat_inv(diaginv+i*nb2, nb);
1080  }
1081  }
1082  }
1083  else {
1084  for (i = 0; i < ROW; ++i) {
1085  fasp_blas_smat_inv(diaginv+i*nb2, nb);
1086  }
1087  }
1088  }
1089  else {
1090  if(use_openmp) {
1091  INT mybegin, myend, myid;
1092 #ifdef _OPENMP
1093 #pragma omp parallel for private(myid,mybegin,myend,i)
1094 #endif
1095  for(myid=0; myid<nthreads; myid++) {
1096  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1097  for(i=mybegin; i<myend; i++) {
1098  diaginv[i] = 1.0 / diaginv[i];
1099  }
1100  }
1101  }
1102  else {
1103  for (i = 0; i < ROW; ++i) {
1104  // zero-diagonal should be tested previously
1105  diaginv[i] = 1.0 / diaginv[i];
1106  }
1107  }
1108  }
1109 
1110  fasp_smoother_dbsr_sor1(A, b, u, order, mark, diaginv, weight);
1111  fasp_mem_free(diaginv);
1112 }
1113 
1136  dvector *b,
1137  dvector *u,
1138  INT order,
1139  INT *mark,
1140  REAL *diaginv,
1141  REAL weight)
1142 {
1143  if (!mark) {
1144  if (order == ASCEND) // smooth ascendingly
1145  {
1146  fasp_smoother_dbsr_sor_ascend(A, b, u, diaginv, weight);
1147  }
1148  else if (order == DESCEND) // smooth descendingly
1149  {
1150  fasp_smoother_dbsr_sor_descend(A, b, u, diaginv, weight);
1151  }
1152  }
1153  // smooth according to the order 'mark' defined by user
1154  else {
1155  fasp_smoother_dbsr_sor_order(A, b, u, diaginv, mark, weight);
1156  }
1157 }
1158 
1177  dvector *b,
1178  dvector *u,
1179  REAL *diaginv,
1180  REAL weight)
1181 {
1182  // members of A
1183  INT ROW = A->ROW;
1184  INT nb = A->nb;
1185  INT *IA = A->IA;
1186  INT *JA = A->JA;
1187  REAL *val = A->val;
1188 
1189  // values of dvector b and u
1190  REAL *b_val = b->val;
1191  REAL *u_val = u->val;
1192 
1193  // local variables
1194  INT nb2 = nb*nb;
1195  INT i,j,k;
1196  INT pb;
1197  REAL rhs = 0.0;
1198  REAL one_minus_weight = 1.0 - weight;
1199 
1200 
1201 #ifdef _OPENMP
1202  // variables for OpenMP
1203  INT myid, mybegin, myend;
1204  INT nthreads = FASP_GET_NUM_THREADS();
1205 #endif
1206 
1207  if (nb == 1) {
1208 #ifdef _OPENMP
1209  if (ROW > OPENMP_HOLDS) {
1210 #pragma omp parallel for private(myid, mybegin, myend, i, rhs, k, j)
1211  for (myid = 0; myid < nthreads; myid++) {
1212  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1213  for (i = mybegin; i < myend; i++) {
1214  rhs = b_val[i];
1215  for (k = IA[i]; k < IA[i+1]; ++k) {
1216  j = JA[k];
1217  if (j != i)
1218  rhs -= val[k]*u_val[j];
1219  }
1220  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1221  }
1222  }
1223  }
1224  else {
1225 #endif
1226  for (i = 0; i < ROW; ++i) {
1227  rhs = b_val[i];
1228  for (k = IA[i]; k < IA[i+1]; ++k) {
1229  j = JA[k];
1230  if (j != i)
1231  rhs -= val[k]*u_val[j];
1232  }
1233  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1234  }
1235 #ifdef _OPENMP
1236  }
1237 #endif
1238  }
1239  else if (nb > 1) {
1240  //REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1241 #ifdef _OPENMP
1242  if (ROW > OPENMP_HOLDS) {
1243  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb*nthreads, sizeof(REAL));
1244  //#pragma omp parallel for private(myid, mybegin, myend, i, pb, b_tmp, k, j)
1245 #pragma omp parallel for private(myid, mybegin, myend, i, pb, k, j)
1246  for (myid = 0; myid < nthreads; myid++) {
1247  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1248  for (i = mybegin; i < myend; i++) {
1249  pb = i*nb;
1250  //memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1251  memcpy(b_tmp+myid*nb, b_val+pb, nb*sizeof(REAL));
1252  for (k = IA[i]; k < IA[i+1]; ++k) {
1253  j = JA[k];
1254  if (j != i)
1255  //fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp+myid*nb, nb);
1256  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1257  }
1258  //fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1259  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp+myid*nb, one_minus_weight, u_val+pb, nb);
1260  }
1261  }
1262  fasp_mem_free(b_tmp);
1263  }
1264  else {
1265 #endif
1266  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1267  for (i = 0; i < ROW; ++i) {
1268  pb = i*nb;
1269  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1270  for (k = IA[i]; k < IA[i+1]; ++k) {
1271  j = JA[k];
1272  if (j != i)
1273  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1274  }
1275  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1276  }
1277  fasp_mem_free(b_tmp);
1278 #ifdef _OPENMP
1279  }
1280 #endif
1281  }
1282  else {
1283  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1284  }
1285 
1286 }
1287 
1306  dvector *b,
1307  dvector *u,
1308  REAL *diaginv,
1309  REAL weight)
1310 {
1311  // members of A
1312  const INT ROW = A->ROW;
1313  const INT nb = A->nb;
1314  const INT nb2 = nb*nb;
1315  const INT *IA = A->IA;
1316  const INT *JA = A->JA;
1317  REAL *val = A->val;
1318  const REAL one_minus_weight = 1.0 - weight;
1319 
1320  // values of dvector b and u
1321  REAL *b_val = b->val;
1322  REAL *u_val = u->val;
1323 
1324  // local variables
1325  INT i,j,k;
1326  INT pb;
1327  REAL rhs = 0.0;
1328 
1329 #ifdef _OPENMP
1330  // variables for OpenMP
1331  INT myid, mybegin, myend;
1332  INT nthreads = FASP_GET_NUM_THREADS();
1333 #endif
1334 
1335  if (nb == 1) {
1336 #ifdef _OPENMP
1337  if (ROW > OPENMP_HOLDS) {
1338 #pragma omp parallel for private(myid, mybegin, myend, i, rhs, k, j)
1339  for (myid = 0; myid < nthreads; myid++) {
1340  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1341  mybegin = ROW-1-mybegin; myend = ROW-1-myend;
1342  for (i = mybegin; i > myend; i--) {
1343  rhs = b_val[i];
1344  for (k = IA[i]; k < IA[i+1]; ++k) {
1345  j = JA[k];
1346  if (j != i)
1347  rhs -= val[k]*u_val[j];
1348  }
1349  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1350  }
1351  }
1352  }
1353  else {
1354 #endif
1355  for (i = ROW-1; i >= 0; i--) {
1356  rhs = b_val[i];
1357  for (k = IA[i]; k < IA[i+1]; ++k) {
1358  j = JA[k];
1359  if (j != i)
1360  rhs -= val[k]*u_val[j];
1361  }
1362  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1363  }
1364 #ifdef _OPENMP
1365  }
1366 #endif
1367  }
1368  else if (nb > 1) {
1369  //REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1370 #ifdef _OPENMP
1371  if (ROW > OPENMP_HOLDS) {
1372  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb*nthreads, sizeof(REAL));
1373  //#pragma omp parallel for private(myid, mybegin, myend, i, pb, b_tmp, k, j)
1374 #pragma omp parallel for private(myid, mybegin, myend, i, pb, k, j)
1375  for (myid = 0; myid < nthreads; myid++) {
1376  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1377  mybegin = ROW-1-mybegin; myend = ROW-1-myend;
1378  for (i = mybegin; i > myend; i--) {
1379  pb = i*nb;
1380  //memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1381  memcpy(b_tmp+myid*nb, b_val+pb, nb*sizeof(REAL));
1382  for (k = IA[i]; k < IA[i+1]; ++k) {
1383  j = JA[k];
1384  if (j != i)
1385  //fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1386  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp+myid*nb, nb);
1387  }
1388  //fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1389  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp+myid*nb, one_minus_weight, u_val+pb, nb);
1390  }
1391  }
1392  fasp_mem_free(b_tmp);
1393  }
1394  else {
1395 #endif
1396  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1397  for (i = ROW-1; i >= 0; i--) {
1398  pb = i*nb;
1399  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1400  for (k = IA[i]; k < IA[i+1]; ++k) {
1401  j = JA[k];
1402  if (j != i)
1403  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1404  }
1405  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1406  }
1407  fasp_mem_free(b_tmp);
1408 #ifdef _OPENMP
1409  }
1410 #endif
1411  }
1412  else {
1413  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1414  }
1415 
1416 }
1417 
1437  dvector *b,
1438  dvector *u,
1439  REAL *diaginv,
1440  INT *mark,
1441  REAL weight)
1442 {
1443  // members of A
1444  const INT ROW = A->ROW;
1445  const INT nb = A->nb;
1446  const INT nb2 = nb*nb;
1447  const INT *IA = A->IA;
1448  const INT *JA = A->JA;
1449  REAL *val = A->val;
1450  const REAL one_minus_weight = 1.0 - weight;
1451 
1452  // values of dvector b and u
1453  REAL *b_val = b->val;
1454  REAL *u_val = u->val;
1455 
1456  // local variables
1457  INT i,j,k;
1458  INT I,pb;
1459  REAL rhs = 0.0;
1460 
1461 #ifdef _OPENMP
1462  // variables for OpenMP
1463  INT myid, mybegin, myend;
1464  INT nthreads = FASP_GET_NUM_THREADS();
1465 #endif
1466 
1467  if (nb == 1) {
1468 #ifdef _OPENMP
1469  if (ROW > OPENMP_HOLDS) {
1470 #pragma omp parallel for private(myid, mybegin, myend, I, i, rhs, k, j)
1471  for (myid = 0; myid < nthreads; myid++) {
1472  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1473  for (I = mybegin; I < myend; ++I) {
1474  i = mark[I];
1475  rhs = b_val[i];
1476  for (k = IA[i]; k < IA[i+1]; ++k) {
1477  j = JA[k];
1478  if (j != i)
1479  rhs -= val[k]*u_val[j];
1480  }
1481  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1482  }
1483  }
1484  }
1485  else {
1486 #endif
1487  for (I = 0; I < ROW; ++I) {
1488  i = mark[I];
1489  rhs = b_val[i];
1490  for (k = IA[i]; k < IA[i+1]; ++k) {
1491  j = JA[k];
1492  if (j != i)
1493  rhs -= val[k]*u_val[j];
1494  }
1495  u_val[i] = one_minus_weight*u_val[i] + weight*(rhs*diaginv[i]);
1496  }
1497 #ifdef _OPENMP
1498  }
1499 #endif
1500  }
1501  else if (nb > 1) {
1502  //REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1503 #ifdef _OPENMP
1504  if (ROW > OPENMP_HOLDS) {
1505  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb*nthreads, sizeof(REAL));
1506  //#pragma omp parallel for private(myid, mybegin, myend, I, i, pb, b_tmp, k, j)
1507 #pragma omp parallel for private(myid, mybegin, myend, I, i, pb, k, j)
1508  for (myid = 0; myid < nthreads; myid++) {
1509  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1510  for (I = mybegin; I < myend; ++I) {
1511  i = mark[I];
1512  pb = i*nb;
1513  //memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1514  memcpy(b_tmp+myid*nb, b_val+pb, nb*sizeof(REAL));
1515  for (k = IA[i]; k < IA[i+1]; ++k) {
1516  j = JA[k];
1517  if (j != i)
1518  //fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1519  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp+myid*nb, nb);
1520  }
1521  //fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1522  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp+myid*nb, one_minus_weight, u_val+pb, nb);
1523  }
1524  }
1525  fasp_mem_free(b_tmp);
1526  }
1527  else {
1528 #endif
1529  REAL *b_tmp = (REAL *)fasp_mem_calloc(nb, sizeof(REAL));
1530  for (I = 0; I < ROW; ++I) {
1531  i = mark[I];
1532  pb = i*nb;
1533  memcpy(b_tmp, b_val+pb, nb*sizeof(REAL));
1534  for (k = IA[i]; k < IA[i+1]; ++k) {
1535  j = JA[k];
1536  if (j != i)
1537  fasp_blas_smat_ymAx(val+k*nb2, u_val+j*nb, b_tmp, nb);
1538  }
1539  fasp_blas_smat_aAxpby(weight, diaginv+nb2*i, b_tmp, one_minus_weight, u_val+pb, nb);
1540  }
1541  fasp_mem_free(b_tmp);
1542 #ifdef _OPENMP
1543  }
1544 #endif
1545  }
1546  else {
1547  fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1548  }
1549 
1550 }
1551 
1567  dvector *b,
1568  dvector *x,
1569  void *data)
1570 {
1571  ILU_data *iludata=(ILU_data *)data;
1572  const INT nb=iludata->nb,m=A->ROW*nb, memneed=5*m;
1573 
1574  REAL *xval = x->val, *bval = b->val;
1575  REAL *zr = iludata->work + 3*m;
1576  REAL *z = zr + m;
1577 
1578  if (iludata->nwork<memneed) goto MEMERR;
1579 
1581  fasp_array_cp(m,bval,zr); fasp_blas_dbsr_aAxpy(-1.0,A,xval,zr);
1582 
1584  fasp_precond_dbsr_ilu(zr,z,iludata);
1585 
1587  fasp_blas_array_axpy(m,1,z,xval);
1588 
1589  return;
1590 
1591 MEMERR:
1592  printf("### ERROR: ILU needs %d memory, only %d available! %s : %d\n",
1593  memneed, iludata->nwork, __FILE__, __LINE__);
1594  fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1595 
1596 }
1597 
1598 /*---------------------------------*/
1599 /*-- End of File --*/
1600 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
void fasp_smoother_dbsr_ilu(dBSRmat *A, dvector *b, dvector *x, void *data)
ILU method as the smoother in solving Au=b with multigrid method.
void fasp_smoother_dbsr_gs1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv)
Gauss-Seidel relaxation.
Definition: smoother_bsr.c:531
#define REAL
Definition: fasp.h:67
void fasp_blas_smat_ymAx(REAL *A, REAL *x, REAL *y, const INT n)
Compute y := y - Ax, where 'A' is a n*n dense matrix.
Definition: blas_smat.c:1207
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 nb
dimension of each sub-block
Definition: fasp_block.h:56
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
REAL * val
actual vector entries
Definition: fasp.h:348
void fasp_smoother_dbsr_sor_order(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark, REAL weight)
SOR relaxation in the user-defined order.
void fasp_blas_smat_aAxpby(const REAL alpha, REAL *A, REAL *x, const REAL beta, REAL *y, const INT n)
Compute y:=alpha*A*x + beta*y.
Definition: blas_smat.c:1308
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
Vector with n entries of REAL type.
Definition: fasp.h:342
void fasp_smoother_dbsr_jacobi(dBSRmat *A, dvector *b, dvector *u)
Jacobi relaxation.
Definition: smoother_bsr.c:33
#define INT
Definition: fasp.h:64
Data for ILU setup.
Definition: fasp.h:400
void fasp_smoother_dbsr_gs_ascend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel relaxation in the ascending order.
Definition: smoother_bsr.c:568
#define DESCEND
Definition: fasp_const.h:232
void fasp_smoother_dbsr_gs_order2(dBSRmat *A, dvector *b, dvector *u, INT *mark, REAL *work)
Gauss-Seidel relaxation in the user-defined order.
Definition: smoother_bsr.c:936
INT nwork
work space size
Definition: fasp.h:421
#define ASCEND
Definition: fasp_const.h:231
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
REAL * val
Definition: fasp_block.h:67
void fasp_smoother_dbsr_jacobi1(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Jacobi relaxation.
Definition: smoother_bsr.c:257
#define OPENMP_HOLDS
Definition: fasp_const.h:248
Main header file for FASP.
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:37
#define ERROR_NUM_BLOCKS
Definition: fasp_const.h:34
void fasp_smoother_dbsr_jacobi_setup(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Setup for jacobi relaxation, fetch the diagonal sub-block matrixes and make them inverse first...
Definition: smoother_bsr.c:148
void fasp_smoother_dbsr_sor_ascend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR relaxation in the ascending order.
void fasp_precond_dbsr_ilu(REAL *r, REAL *z, void *data)
ILU preconditioner.
Definition: precond_bsr.c:306
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
void fasp_smoother_dbsr_gs_descend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel relaxation in the descending order.
Definition: smoother_bsr.c:712
void fasp_smoother_dbsr_sor(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL weight)
SOR relaxation.
void fasp_smoother_dbsr_gs(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark)
Gauss-Seidel relaxation.
Definition: smoother_bsr.c:411
void fasp_smoother_dbsr_gs_order1(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark)
Gauss-Seidel relaxation in the user-defined order.
Definition: smoother_bsr.c:858
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:47
INT nb
block size for BSR type only
Definition: fasp.h:418
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
INT * JA
Definition: fasp_block.h:74
REAL * work
work space
Definition: fasp.h:424
#define FALSE
Definition: fasp_const.h:68
void fasp_smoother_dbsr_gs_ascend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the ascending order.
Definition: smoother_bsr.c:641
INT * IA
integer array of row pointers, the size is ROW+1
Definition: fasp_block.h:70
void fasp_smoother_dbsr_sor1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv, REAL weight)
SOR relaxation.
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
void fasp_smoother_dbsr_gs_descend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the descending order.
Definition: smoother_bsr.c:786
void fasp_blas_dbsr_aAxpy(const REAL alpha, dBSRmat *A, REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: blas_bsr.c:337
void fasp_smoother_dbsr_sor_descend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR relaxation in the descending order.