Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
ItrSmootherSTR.c
Go to the documentation of this file.
1
15#include <math.h>
16
17#include "fasp.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Declare Private Functions --*/
22/*---------------------------------*/
23
24static void blkcontr2 (INT, INT, INT, INT, REAL *, REAL *, REAL *);
25static void aAxpby (REAL, REAL, INT, REAL *, REAL *, REAL *);
26
27/*---------------------------------*/
28/*-- Public Functions --*/
29/*---------------------------------*/
30
43void fasp_smoother_dstr_jacobi (dSTRmat *A,
44 dvector *b,
45 dvector *u)
46{
47 INT nc = A->nc; // size of each block (number of components)
48 INT ngrid = A->ngrid; // number of grids
49 REAL *diag = A->diag; // Diagonal entries
50 REAL *diaginv = NULL; // Diagonal inverse, same size and storage scheme as A->diag
51
52 INT nc2 = nc*nc;
53 INT size = nc2*ngrid;
54 INT block = 0;
55 INT start = 0;
56
57 if (nc > 1) {
58 // allocate memory
59 diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
60
61 // diaginv = diag;
62 fasp_darray_cp(size,diag,diaginv);
63
64 // generate diaginv
65 for (block = 0; block < ngrid; block ++) {
66 fasp_smat_inv(diaginv+start, nc);
67 start += nc2;
68 }
69 }
70
71 fasp_smoother_dstr_jacobi1(A, b, u, diaginv);
72
73 fasp_mem_free(diaginv); diaginv = NULL;
74}
75
76
92void fasp_smoother_dstr_jacobi1 (dSTRmat *A,
93 dvector *b,
94 dvector *u,
95 REAL *diaginv)
96{
97 // information of A
98 INT ngrid = A->ngrid; // number of grids
99 INT nc = A->nc; // size of each block (number of components)
100 INT nband = A->nband; // number of off-diag band
101 INT *offsets = A->offsets; // offsets of the off-diagals
102 REAL *diag = A->diag; // Diagonal entries
103 REAL **offdiag = A->offdiag; // Off-diagonal entries
104
105 // values of dvector b and u
106 REAL *b_val = b->val;
107 REAL *u_val = u->val;
108
109 // local variables
110 INT block = 0;
111 INT point = 0;
112 INT band = 0;
113 INT width = 0;
114 INT size = nc*ngrid;
115 INT nc2 = nc*nc;
116 INT start = 0;
117 INT column = 0;
118 INT start_data = 0;
119 INT start_DATA = 0;
120 INT start_vecb = 0;
121 INT start_vecu = 0;
122
123 // auxiliary array
124 REAL *b_tmp = NULL;
125
126 // this should be done once and for all!!
127 b_tmp = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
128
129 // b_tmp = b_val
130 fasp_darray_cp(size,b_val,b_tmp);
131
132 // It's not necessary to assign the smoothing order since the results doesn't depend on it
133 if (nc == 1) {
134 for (point = 0; point < ngrid; point ++) {
135 for (band = 0; band < nband; band ++) {
136 width = offsets[band];
137 column = point + width;
138 if (width < 0) {
139 if (column >= 0) {
140 b_tmp[point] -= offdiag[band][column]*u_val[column];
141 }
142 }
143 else {// width > 0
144 if (column < ngrid) {
145 b_tmp[point] -= offdiag[band][point]*u_val[column];
146 }
147 }
148 } // end for band
149 } // end for point
150
151 for (point = 0; point < ngrid; point ++) {
152 // zero-diagonal should be tested previously
153 u_val[point] = b_tmp[point] / diag[point];
154 }
155 } // end if (nc == 1)
156 else if (nc > 1) {
157 for (block = 0; block < ngrid; block ++) {
158 start_DATA = nc2*block;
159 start_vecb = nc*block;
160 for (band = 0; band < nband; band ++) {
161 width = offsets[band];
162 column = block + width;
163 if (width < 0) {
164 if (column >= 0) {
165 start_data = nc2*column;
166 start_vecu = nc*column;
167 blkcontr2( start_data, start_vecu, start_vecb,
168 nc, offdiag[band], u_val, b_tmp );
169 }
170 }
171 else {// width > 0
172 if (column < ngrid) {
173 start_vecu = nc*column;
174 blkcontr2( start_DATA, start_vecu, start_vecb,
175 nc, offdiag[band], u_val, b_tmp );
176 }
177 }
178 } // end for band
179 } // end for block
180
181 for (block = 0; block < ngrid; block ++) {
182 start = nc*block;
183 fasp_blas_smat_mxv(diaginv+nc2*block, b_tmp+start, u_val+start, nc);
184 }
185 } // end else if (nc > 1)
186 else {
187 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
188 return;
189 }
190
191 fasp_mem_free(b_tmp); b_tmp = NULL;
192}
193
217void fasp_smoother_dstr_gs (dSTRmat *A,
218 dvector *b,
219 dvector *u,
220 const INT order,
221 INT *mark)
222{
223 INT nc = A->nc; // size of each block (number of components)
224 INT ngrid = A->ngrid; // number of grids
225 REAL *diag = A->diag; // Diagonal entries
226 REAL *diaginv = NULL; // Diagonal inverse(when nc>1),same size and storage scheme as A->diag
227
228 INT nc2 = nc*nc;
229 INT size = nc2*ngrid;
230 INT block = 0;
231 INT start = 0;
232
233 if (nc > 1) {
234 // allocate memory
235 diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
236
237 // diaginv = diag;
238 fasp_darray_cp(size,diag,diaginv);
239
240 // generate diaginv
241 for (block = 0; block < ngrid; block ++) {
242 fasp_smat_inv(diaginv+start, nc);
243 start += nc2;
244 }
245 }
246
247 fasp_smoother_dstr_gs1(A, b, u, order, mark, diaginv);
248
249 fasp_mem_free(diaginv); diaginv = NULL;
250}
251
277void fasp_smoother_dstr_gs1 (dSTRmat *A,
278 dvector *b,
279 dvector *u,
280 const INT order,
281 INT *mark,
282 REAL *diaginv)
283{
284
285 if (!mark) {
286 if (order == ASCEND) // smooth ascendingly
287 {
288 fasp_smoother_dstr_gs_ascend(A, b, u, diaginv);
289 }
290 else if (order == DESCEND) // smooth descendingly
291 {
292 fasp_smoother_dstr_gs_descend(A, b, u, diaginv);
293 }
294 }
295 else {
296 if (order == USERDEFINED) // smooth according to the order 'mark' defined by user
297 {
298 fasp_smoother_dstr_gs_order(A, b, u, diaginv, mark);
299 }
300 else // smooth according to 'mark', where 'mark' is a CF_marker array
301 {
302 fasp_smoother_dstr_gs_cf(A, b, u, diaginv, mark, order);
303 }
304 }
305}
306
322void fasp_smoother_dstr_gs_ascend (dSTRmat *A,
323 dvector *b,
324 dvector *u,
325 REAL *diaginv)
326{
327 // information of A
328 INT ngrid = A->ngrid; // number of grids
329 INT nc = A->nc; // size of each block (number of components)
330 INT nband = A->nband; // number of off-diag band
331 INT *offsets = A->offsets; // offsets of the off-diagals
332 REAL *diag = A->diag; // Diagonal entries
333 REAL **offdiag = A->offdiag; // Off-diagonal entries
334
335 // values of dvector b and u
336 REAL *b_val = b->val;
337 REAL *u_val = u->val;
338
339 // local variables
340 INT block = 0;
341 INT point = 0;
342 INT band = 0;
343 INT width = 0;
344 INT nc2 = nc*nc;
345 INT ncb = 0;
346 INT column = 0;
347 INT start_data = 0;
348 INT start_DATA = 0;
349 INT start_vecu = 0;
350 REAL rhs = 0.0;
351
352 // auxiliary array(nc*1 vector)
353 REAL *vec_tmp = NULL;
354
355 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
356
357 if (nc == 1) {
358 for (point = 0; point < ngrid; point ++) {
359 rhs = b_val[point];
360 for (band = 0; band < nband; band ++) {
361 width = offsets[band];
362 column = point + width;
363 if (width < 0) {
364 if (column >= 0) {
365 rhs -= offdiag[band][column]*u_val[column];
366 }
367 }
368 else { // width > 0
369 if (column < ngrid) {
370 rhs -= offdiag[band][point]*u_val[column];
371 }
372 }
373 } // end for band
374
375 // zero-diagonal should be tested previously
376 u_val[point] = rhs / diag[point];
377
378 } // end for point
379
380 } // end if (nc == 1)
381
382 else if (nc > 1) {
383 for (block = 0; block < ngrid; block ++) {
384 ncb = nc*block;
385 for (point = 0; point < nc; point ++) {
386 vec_tmp[point] = b_val[ncb+point];
387 }
388 start_DATA = nc2*block;
389 for (band = 0; band < nband; band ++) {
390 width = offsets[band];
391 column = block + width;
392 if (width < 0) {
393 if (column >= 0) {
394 start_data = nc2*column;
395 start_vecu = nc*column;
396 blkcontr2( start_data, start_vecu, 0, nc,
397 offdiag[band], u_val, vec_tmp );
398 }
399 }
400 else { // width > 0
401 if (column < ngrid) {
402 start_vecu = nc*column;
403 blkcontr2( start_DATA, start_vecu, 0, nc,
404 offdiag[band], u_val, vec_tmp );
405 }
406 }
407 } // end for band
408
409 // subblock smoothing
410 fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
411
412 } // end for block
413
414 } // end else if (nc > 1)
415 else {
416 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
417 return;
418 }
419
420 fasp_mem_free(vec_tmp); vec_tmp = NULL;
421}
422
438void fasp_smoother_dstr_gs_descend (dSTRmat *A,
439 dvector *b,
440 dvector *u,
441 REAL *diaginv)
442{
443 // information of A
444 INT ngrid = A->ngrid; // number of grids
445 INT nc = A->nc; // size of each block (number of components)
446 INT nband = A->nband ; // number of off-diag band
447 INT *offsets = A->offsets; // offsets of the off-diagals
448 REAL *diag = A->diag; // Diagonal entries
449 REAL **offdiag = A->offdiag; // Off-diagonal entries
450
451 // values of dvector b and u
452 REAL *b_val = b->val;
453 REAL *u_val = u->val;
454
455 // local variables
456 INT block = 0;
457 INT point = 0;
458 INT band = 0;
459 INT width = 0;
460 INT nc2 = nc*nc;
461 INT ncb = 0;
462 INT column = 0;
463 INT start_data = 0;
464 INT start_DATA = 0;
465 INT start_vecu = 0;
466 REAL rhs = 0.0;
467
468 // auxiliary array(nc*1 vector)
469 REAL *vec_tmp = NULL;
470
471 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
472
473 if (nc == 1) {
474 for (point = ngrid-1; point >= 0; point --) {
475 rhs = b_val[point];
476 for (band = 0; band < nband; band ++) {
477 width = offsets[band];
478 column = point + width;
479 if (width < 0) {
480 if (column >= 0) {
481 rhs -= offdiag[band][column]*u_val[column];
482 }
483 }
484 else { // width > 0
485 if (column < ngrid) {
486 rhs -= offdiag[band][point]*u_val[column];
487 }
488 }
489 } // end for band
490
491 // zero-diagonal should be tested previously
492 u_val[point] = rhs / diag[point];
493
494 } // end for point
495
496 } // end if (nc == 1)
497
498 else if (nc > 1) {
499 for (block = ngrid-1; block >= 0; block --) {
500 ncb = nc*block;
501 for (point = 0; point < nc; point ++) {
502 vec_tmp[point] = b_val[ncb+point];
503 }
504 start_DATA = nc2*block;
505 for (band = 0; band < nband; band ++) {
506 width = offsets[band];
507 column = block + width;
508 if (width < 0) {
509 if (column >= 0) {
510 start_data = nc2*column;
511 start_vecu = nc*column;
512 blkcontr2( start_data, start_vecu, 0, nc,
513 offdiag[band], u_val, vec_tmp );
514 }
515 }
516 else { // width > 0
517 if (column < ngrid) {
518 start_vecu = nc*column;
519 blkcontr2( start_DATA, start_vecu, 0, nc,
520 offdiag[band], u_val, vec_tmp );
521 }
522 }
523 } // end for band
524
525 // subblock smoothing
526 fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
527
528 } // end for block
529
530 } // end else if (nc > 1)
531
532 else {
533 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
534 return;
535 }
536
537 fasp_mem_free(vec_tmp); vec_tmp = NULL;
538}
539
556void fasp_smoother_dstr_gs_order (dSTRmat *A,
557 dvector *b,
558 dvector *u,
559 REAL *diaginv,
560 INT *mark)
561{
562 // information of A
563 INT ngrid = A->ngrid; // number of grids
564 INT nc = A->nc; // size of each block (number of components)
565 INT nband = A->nband; // number of off-diag band
566 INT *offsets = A->offsets; // offsets of the off-diagals
567 REAL *diag = A->diag; // Diagonal entries
568 REAL **offdiag = A->offdiag; // Off-diagonal entries
569
570 // values of dvector b and u
571 REAL *b_val = b->val;
572 REAL *u_val = u->val;
573
574 // local variables
575 INT block = 0;
576 INT point = 0;
577 INT band = 0;
578 INT width = 0;
579 INT nc2 = nc*nc;
580 INT ncb = 0;
581 INT index = 0;
582 INT column = 0;
583 INT start_data = 0;
584 INT start_DATA = 0;
585 INT start_vecu = 0;
586 REAL rhs = 0.0;
587
588 // auxiliary array(nc*1 vector)
589 REAL *vec_tmp = NULL;
590
591 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
592
593 if (nc == 1) {
594 for (index = 0; index < ngrid; index ++) {
595 point = mark[index];
596 rhs = b_val[point];
597 for (band = 0; band < nband; band ++) {
598 width = offsets[band];
599 column = point + width;
600 if (width < 0) {
601 if (column >= 0) {
602 rhs -= offdiag[band][column]*u_val[column];
603 }
604 }
605 else { // width > 0
606 if (column < ngrid) {
607 rhs -= offdiag[band][point]*u_val[column];
608 }
609 }
610 } // end for band
611
612 // zero-diagonal should be tested previously
613 u_val[point] = rhs / diag[point];
614
615 } // end for index
616
617 } // end if (nc == 1)
618
619 else if (nc > 1) {
620 for (index = 0; index < ngrid; index ++) {
621 block = mark[index];
622 ncb = nc*block;
623 for (point = 0; point < nc; point ++) {
624 vec_tmp[point] = b_val[ncb+point];
625 }
626 start_DATA = nc2*block;
627 for (band = 0; band < nband; band ++) {
628 width = offsets[band];
629 column = block + width;
630 if (width < 0) {
631 if (column >= 0) {
632 start_data = nc2*column;
633 start_vecu = nc*column;
634 blkcontr2( start_data, start_vecu, 0, nc,
635 offdiag[band], u_val, vec_tmp );
636 }
637 }
638 else { // width > 0
639 if (column < ngrid) {
640 start_vecu = nc*column;
641 blkcontr2( start_DATA, start_vecu, 0, nc,
642 offdiag[band], u_val, vec_tmp );
643 }
644 }
645 } // end for band
646
647 // subblock smoothing
648 fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
649
650 } // end for index
651
652 } // end else if (nc > 1)
653 else {
654 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
655 return;
656 }
657
658 fasp_mem_free(vec_tmp); vec_tmp = NULL;
659}
660
680void fasp_smoother_dstr_gs_cf (dSTRmat *A,
681 dvector *b,
682 dvector *u,
683 REAL *diaginv,
684 INT *mark,
685 const INT order)
686{
687 // information of A
688 INT ngrid = A->ngrid; // number of grids
689 INT nc = A->nc; // size of each block (number of components)
690 INT nband = A->nband ; // number of off-diag band
691 INT *offsets = A->offsets; // offsets of the off-diagals
692 REAL *diag = A->diag; // Diagonal entries
693 REAL **offdiag = A->offdiag; // Off-diagonal entries
694
695 // values of dvector b and u
696 REAL *b_val = b->val;
697 REAL *u_val = u->val;
698
699 // local variables
700 INT block = 0;
701 INT point = 0;
702 INT band = 0;
703 INT width = 0;
704 INT nc2 = nc*nc;
705 INT ncb = 0;
706 INT column = 0;
707 INT start_data = 0;
708 INT start_DATA = 0;
709 INT start_vecu = 0;
710 INT FIRST = order; // which kind of points to be smoothed firstly?
711 INT SECOND = -order; // which kind of points to be smoothed secondly?
712
713 REAL rhs = 0.0;
714
715 // auxiliary array(nc*1 vector)
716 REAL *vec_tmp = NULL;
717
718 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
719
720 if (nc == 1) {
721 // deal with the points marked FIRST
722 for (point = 0; point < ngrid; point ++) {
723 if (mark[point] == FIRST) {
724 rhs = b_val[point];
725 for (band = 0; band < nband; band ++) {
726 width = offsets[band];
727 column = point + width;
728 if (width < 0) {
729 if (column >= 0) {
730 rhs -= offdiag[band][column]*u_val[column];
731 }
732 }
733 else { // width > 0
734 if (column < ngrid) {
735 rhs -= offdiag[band][point]*u_val[column];
736 }
737 }
738 } // end for band
739
740 // zero-diagonal should be tested previously
741 u_val[point] = rhs / diag[point];
742 } // end if (mark[point] == FIRST)
743 } // end for point
744
745 // deal with the points marked SECOND
746 for (point = 0; point < ngrid; point ++) {
747 if (mark[point] == SECOND) {
748 rhs = b_val[point];
749 for (band = 0; band < nband; band ++) {
750 width = offsets[band];
751 column = point + width;
752 if (width < 0) {
753 if (column >= 0) {
754 rhs -= offdiag[band][column]*u_val[column];
755 }
756 }
757 else { // width > 0
758 if (column < ngrid) {
759 rhs -= offdiag[band][point]*u_val[column];
760 }
761 }
762 } // end for band
763
764 // zero-diagonal should be tested previously
765 u_val[point] = rhs / diag[point];
766 } // end if (mark[point] == SECOND)
767 } // end for point
768
769 } // end if (nc == 1)
770
771 else if (nc > 1) {
772 // deal with the blocks marked FIRST
773 for (block = 0; block < ngrid; block ++) {
774 if (mark[block] == FIRST) {
775 ncb = nc*block;
776 for (point = 0; point < nc; point ++) {
777 vec_tmp[point] = b_val[ncb+point];
778 }
779 start_DATA = nc2*block;
780 for (band = 0; band < nband; band ++) {
781 width = offsets[band];
782 column = block + width;
783 if (width < 0) {
784 if (column >= 0) {
785 start_data = nc2*column;
786 start_vecu = nc*column;
787 blkcontr2( start_data, start_vecu, 0, nc,
788 offdiag[band], u_val, vec_tmp );
789 }
790 }
791 else { // width > 0
792 if (column < ngrid) {
793 start_vecu = nc*column;
794 blkcontr2( start_DATA, start_vecu, 0, nc,
795 offdiag[band], u_val, vec_tmp );
796 }
797 }
798 } // end for band
799
800 // subblock smoothing
801 fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
802 } // end if (mark[block] == FIRST)
803
804 } // end for block
805
806 // deal with the blocks marked SECOND
807 for (block = 0; block < ngrid; block ++) {
808 if (mark[block] == SECOND) {
809 ncb = nc*block;
810 for (point = 0; point < nc; point ++) {
811 vec_tmp[point] = b_val[ncb+point];
812 }
813 start_DATA = nc2*block;
814 for (band = 0; band < nband; band ++) {
815 width = offsets[band];
816 column = block + width;
817 if (width < 0) {
818 if (column >= 0) {
819 start_data = nc2*column;
820 start_vecu = nc*column;
821 blkcontr2( start_data, start_vecu, 0, nc,
822 offdiag[band], u_val, vec_tmp );
823 }
824 }
825 else { // width > 0
826 if (column < ngrid) {
827 start_vecu = nc*column;
828 blkcontr2( start_DATA, start_vecu, 0, nc,
829 offdiag[band], u_val, vec_tmp );
830 }
831 }
832 } // end for band
833
834 // subblock smoothing
835 fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
836 } // end if (mark[block] == SECOND)
837
838 } // end for block
839
840 } // end else if (nc > 1)
841 else {
842 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
843 return;
844 }
845
846 fasp_mem_free(vec_tmp); vec_tmp = NULL;
847}
848
873void fasp_smoother_dstr_sor (dSTRmat *A,
874 dvector *b,
875 dvector *u,
876 const INT order,
877 INT *mark,
878 const REAL weight)
879{
880 INT nc = A->nc; // size of each block (number of components)
881 INT ngrid = A->ngrid; // number of grids
882 REAL *diag = A->diag; // Diagonal entries
883 REAL *diaginv = NULL; // Diagonal inverse(when nc>1),same size and storage scheme as A->diag
884
885 INT nc2 = nc*nc;
886 INT size = nc2*ngrid;
887 INT block = 0;
888 INT start = 0;
889
890 if (nc > 1) {
891 // allocate memory
892 diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
893
894 // diaginv = diag;
895 fasp_darray_cp(size,diag,diaginv);
896
897 // generate diaginv
898 for (block = 0; block < ngrid; block ++) {
899 fasp_smat_inv(diaginv+start, nc);
900 start += nc2;
901 }
902 }
903
904 fasp_smoother_dstr_sor1(A, b, u, order, mark, diaginv, weight);
905
906 fasp_mem_free(diaginv); diaginv = NULL;
907}
908
935void fasp_smoother_dstr_sor1 (dSTRmat *A,
936 dvector *b,
937 dvector *u,
938 const INT order,
939 INT *mark,
940 REAL *diaginv,
941 const REAL weight)
942{
943 if (!mark) {
944 if (order == ASCEND) // smooth ascendingly
945 {
946 fasp_smoother_dstr_sor_ascend(A, b, u, diaginv, weight);
947 }
948 else if (order == DESCEND) // smooth descendingly
949 {
950 fasp_smoother_dstr_sor_descend(A, b, u, diaginv, weight);
951 }
952 }
953 else {
954 if (order == USERDEFINED) // smooth according to the order 'mark' defined by user
955 {
956 fasp_smoother_dstr_sor_order(A, b, u, diaginv, mark, weight);
957 }
958 else // smooth according to 'mark', where 'mark' is a CF_marker array
959 {
960 fasp_smoother_dstr_sor_cf(A, b, u, diaginv, mark, order, weight);
961 }
962 }
963}
964
981void fasp_smoother_dstr_sor_ascend (dSTRmat *A,
982 dvector *b,
983 dvector *u,
984 REAL *diaginv,
985 REAL weight)
986{
987 // information of A
988 INT ngrid = A->ngrid; // number of grids
989 INT nc = A->nc; // size of each block (number of components)
990 INT nband = A->nband ; // number of off-diag band
991 INT *offsets = A->offsets; // offsets of the off-diagals
992 REAL *diag = A->diag; // Diagonal entries
993 REAL **offdiag = A->offdiag; // Off-diagonal entries
994
995 // values of dvector b and u
996 REAL *b_val = b->val;
997 REAL *u_val = u->val;
998
999 // local variables
1000 INT block = 0;
1001 INT point = 0;
1002 INT band = 0;
1003 INT width = 0;
1004 INT nc2 = nc*nc;
1005 INT ncb = 0;
1006 INT column = 0;
1007 INT start_data = 0;
1008 INT start_DATA = 0;
1009 INT start_vecu = 0;
1010 REAL rhs = 0.0;
1011 REAL one_minus_weight = 1.0 - weight;
1012
1013 // auxiliary array(nc*1 vector)
1014 REAL *vec_tmp = NULL;
1015
1016 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1017
1018 if (nc == 1) {
1019 for (point = 0; point < ngrid; point ++) {
1020 rhs = b_val[point];
1021 for (band = 0; band < nband; band ++) {
1022 width = offsets[band];
1023 column = point + width;
1024 if (width < 0) {
1025 if (column >= 0) {
1026 rhs -= offdiag[band][column]*u_val[column];
1027 }
1028 }
1029 else { // width > 0
1030 if (column < ngrid) {
1031 rhs -= offdiag[band][point]*u_val[column];
1032 }
1033 }
1034 } // end for band
1035
1036 // zero-diagonal should be tested previously
1037 u_val[point] = one_minus_weight*u_val[point] +
1038 weight*(rhs / diag[point]);
1039
1040 } // end for point
1041
1042 } // end if (nc == 1)
1043
1044 else if (nc > 1) {
1045 for (block = 0; block < ngrid; block ++) {
1046 ncb = nc*block;
1047 for (point = 0; point < nc; point ++) {
1048 vec_tmp[point] = b_val[ncb+point];
1049 }
1050 start_DATA = nc2*block;
1051 for (band = 0; band < nband; band ++) {
1052 width = offsets[band];
1053 column = block + width;
1054 if (width < 0) {
1055 if (column >= 0) {
1056 start_data = nc2*column;
1057 start_vecu = nc*column;
1058 blkcontr2( start_data, start_vecu, 0, nc,
1059 offdiag[band], u_val, vec_tmp );
1060 }
1061 }
1062 else { // width > 0
1063 if (column < ngrid) {
1064 start_vecu = nc*column;
1065 blkcontr2( start_DATA, start_vecu, 0, nc,
1066 offdiag[band], u_val, vec_tmp );
1067 }
1068 }
1069 } // end for band
1070
1071 // subblock smoothing
1072 aAxpby(weight, one_minus_weight, nc,
1073 diaginv+start_DATA, vec_tmp, u_val+nc*block);
1074
1075 } // end for block
1076
1077 } // end else if (nc > 1)
1078 else {
1079 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
1080 return;
1081 }
1082
1083 fasp_mem_free(vec_tmp); vec_tmp = NULL;
1084}
1085
1102void fasp_smoother_dstr_sor_descend (dSTRmat *A,
1103 dvector *b,
1104 dvector *u,
1105 REAL *diaginv,
1106 REAL weight)
1107{
1108 // information of A
1109 INT ngrid = A->ngrid; // number of grids
1110 INT nc = A->nc; // size of each block (number of components)
1111 INT nband = A->nband ; // number of off-diag band
1112 INT *offsets = A->offsets; // offsets of the off-diagals
1113 REAL *diag = A->diag; // Diagonal entries
1114 REAL **offdiag = A->offdiag; // Off-diagonal entries
1115
1116 // values of dvector b and u
1117 REAL *b_val = b->val;
1118 REAL *u_val = u->val;
1119
1120 // local variables
1121 INT block = 0;
1122 INT point = 0;
1123 INT band = 0;
1124 INT width = 0;
1125 INT nc2 = nc*nc;
1126 INT ncb = 0;
1127 INT column = 0;
1128 INT start_data = 0;
1129 INT start_DATA = 0;
1130 INT start_vecu = 0;
1131 REAL rhs = 0.0;
1132 REAL one_minus_weight = 1.0 - weight;
1133
1134 // auxiliary array(nc*1 vector)
1135 REAL *vec_tmp = NULL;
1136
1137 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1138
1139 if (nc == 1) {
1140 for (point = ngrid-1; point >= 0; point --) {
1141 rhs = b_val[point];
1142 for (band = 0; band < nband; band ++) {
1143 width = offsets[band];
1144 column = point + width;
1145 if (width < 0) {
1146 if (column >= 0) {
1147 rhs -= offdiag[band][column]*u_val[column];
1148 }
1149 }
1150 else { // width > 0
1151 if (column < ngrid) {
1152 rhs -= offdiag[band][point]*u_val[column];
1153 }
1154 }
1155 } // end for band
1156
1157 // zero-diagonal should be tested previously
1158 u_val[point] = one_minus_weight*u_val[point] +
1159 weight*(rhs / diag[point]);
1160
1161 } // end for point
1162
1163 } // end if (nc == 1)
1164
1165 else if (nc > 1) {
1166 for (block = ngrid-1; block >= 0; block --) {
1167 ncb = nc*block;
1168 for (point = 0; point < nc; point ++) {
1169 vec_tmp[point] = b_val[ncb+point];
1170 }
1171 start_DATA = nc2*block;
1172 for (band = 0; band < nband; band ++) {
1173 width = offsets[band];
1174 column = block + width;
1175 if (width < 0) {
1176 if (column >= 0) {
1177 start_data = nc2*column;
1178 start_vecu = nc*column;
1179 blkcontr2( start_data, start_vecu, 0, nc,
1180 offdiag[band], u_val, vec_tmp );
1181 }
1182 }
1183 else { // width > 0
1184 if (column < ngrid) {
1185 start_vecu = nc*column;
1186 blkcontr2( start_DATA, start_vecu, 0, nc,
1187 offdiag[band], u_val, vec_tmp );
1188 }
1189 }
1190 } // end for band
1191
1192 // subblock smoothing
1193 aAxpby(weight, one_minus_weight, nc,
1194 diaginv+start_DATA, vec_tmp, u_val+nc*block);
1195
1196 } // end for block
1197
1198 } // end else if (nc > 1)
1199 else {
1200 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
1201 return;
1202 }
1203
1204 fasp_mem_free(vec_tmp); vec_tmp = NULL;
1205}
1206
1224void fasp_smoother_dstr_sor_order (dSTRmat *A,
1225 dvector *b,
1226 dvector *u,
1227 REAL *diaginv,
1228 INT *mark,
1229 REAL weight)
1230{
1231 // information of A
1232 INT ngrid = A->ngrid; // number of grids
1233 INT nc = A->nc; // size of each block (number of components)
1234 INT nband = A->nband ; // number of off-diag band
1235 INT *offsets = A->offsets; // offsets of the off-diagals
1236 REAL *diag = A->diag; // Diagonal entries
1237 REAL **offdiag = A->offdiag; // Off-diagonal entries
1238
1239 // values of dvector b and u
1240 REAL *b_val = b->val;
1241 REAL *u_val = u->val;
1242
1243 // local variables
1244 INT block = 0;
1245 INT point = 0;
1246 INT band = 0;
1247 INT width = 0;
1248 INT nc2 = nc*nc;
1249 INT ncb = 0;
1250 INT column = 0;
1251 INT index = 0;
1252 INT start_data = 0;
1253 INT start_DATA = 0;
1254 INT start_vecu = 0;
1255 REAL rhs = 0.0;
1256 REAL one_minus_weight = 1.0 - weight;
1257
1258 // auxiliary array(nc*1 vector)
1259 REAL *vec_tmp = NULL;
1260
1261 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1262
1263 if (nc == 1) {
1264 for (index = 0; index < ngrid; index ++) {
1265 point = mark[index];
1266 rhs = b_val[point];
1267 for (band = 0; band < nband; band ++) {
1268 width = offsets[band];
1269 column = point + width;
1270 if (width < 0) {
1271 if (column >= 0) {
1272 rhs -= offdiag[band][column]*u_val[column];
1273 }
1274 }
1275 else { // width > 0
1276 if (column < ngrid) {
1277 rhs -= offdiag[band][point]*u_val[column];
1278 }
1279 }
1280 } // end for band
1281
1282 // zero-diagonal should be tested previously
1283 u_val[point] = one_minus_weight*u_val[point] +
1284 weight*(rhs / diag[point]);
1285
1286 } // end for index
1287
1288 } // end if (nc == 1)
1289
1290 else if (nc > 1) {
1291 for (index = 0; index < ngrid; index ++) {
1292 block = mark[index];
1293 ncb = nc*block;
1294 for (point = 0; point < nc; point ++) {
1295 vec_tmp[point] = b_val[ncb+point];
1296 }
1297 start_DATA = nc2*block;
1298 for (band = 0; band < nband; band ++) {
1299 width = offsets[band];
1300 column = block + width;
1301 if (width < 0) {
1302 if (column >= 0) {
1303 start_data = nc2*column;
1304 start_vecu = nc*column;
1305 blkcontr2( start_data, start_vecu, 0, nc,
1306 offdiag[band], u_val, vec_tmp );
1307 }
1308 }
1309 else { // width > 0
1310 if (column < ngrid) {
1311 start_vecu = nc*column;
1312 blkcontr2( start_DATA, start_vecu, 0, nc,
1313 offdiag[band], u_val, vec_tmp );
1314 }
1315 }
1316 } // end for band
1317
1318 // subblock smoothing
1319 aAxpby(weight, one_minus_weight, nc,
1320 diaginv+start_DATA, vec_tmp, u_val+nc*block);
1321
1322 } // end for index
1323
1324 } // end else if (nc > 1)
1325
1326 else {
1327 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
1328 return;
1329 }
1330
1331 fasp_mem_free(vec_tmp); vec_tmp = NULL;
1332}
1333
1355void fasp_smoother_dstr_sor_cf (dSTRmat *A,
1356 dvector *b,
1357 dvector *u,
1358 REAL *diaginv,
1359 INT *mark,
1360 const INT order,
1361 const REAL weight)
1362{
1363 // information of A
1364 INT ngrid = A->ngrid; // number of grids
1365 INT nc = A->nc; // size of each block (number of components)
1366 INT nband = A->nband ; // number of off-diag band
1367 INT *offsets = A->offsets; // offsets of the off-diagals
1368 REAL *diag = A->diag; // Diagonal entries
1369 REAL **offdiag = A->offdiag; // Off-diagonal entries
1370
1371 // values of dvector b and u
1372 REAL *b_val = b->val;
1373 REAL *u_val = u->val;
1374
1375 // local variables
1376 INT block = 0;
1377 INT point = 0;
1378 INT band = 0;
1379 INT width = 0;
1380 INT nc2 = nc*nc;
1381 INT ncb = 0;
1382 INT column = 0;
1383 INT start_data = 0;
1384 INT start_DATA = 0;
1385 INT start_vecu = 0;
1386 REAL rhs = 0.0;
1387 REAL one_minus_weight = 1.0 - weight;
1388 INT FIRST = order; // which kind of points to be smoothed firstly?
1389 INT SECOND = -order; // which kind of points to be smoothed secondly?
1390
1391 // auxiliary array(nc*1 vector)
1392 REAL *vec_tmp = NULL;
1393
1394 vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1395
1396 if (nc == 1) {
1397 // deal with the points marked FIRST
1398 for (point = 0; point < ngrid; point ++) {
1399 if (mark[point] == FIRST) {
1400 rhs = b_val[point];
1401 for (band = 0; band < nband; band ++) {
1402 width = offsets[band];
1403 column = point + width;
1404 if (width < 0) {
1405 if (column >= 0) {
1406 rhs -= offdiag[band][column]*u_val[column];
1407 }
1408 }
1409 else { // width > 0
1410 if (column < ngrid) {
1411 rhs -= offdiag[band][point]*u_val[column];
1412 }
1413 }
1414 } // end for band
1415
1416 // zero-diagonal should be tested previously
1417 u_val[point] = one_minus_weight*u_val[point] +
1418 weight*(rhs / diag[point]);
1419
1420 } // end if (mark[point] == FIRST)
1421 } // end for point
1422
1423 // deal with the points marked SECOND
1424 for (point = 0; point < ngrid; point ++) {
1425 if (mark[point] == SECOND) {
1426 rhs = b_val[point];
1427 for (band = 0; band < nband; band ++) {
1428 width = offsets[band];
1429 column = point + width;
1430 if (width < 0) {
1431 if (column >= 0) {
1432 rhs -= offdiag[band][column]*u_val[column];
1433 }
1434 }
1435 else { // width > 0
1436 if (column < ngrid) {
1437 rhs -= offdiag[band][point]*u_val[column];
1438 }
1439 }
1440 } // end for band
1441
1442 // zero-diagonal should be tested previously
1443 u_val[point] = rhs / diag[point];
1444 } // end if (mark[point] == SECOND)
1445 } // end for point
1446
1447 } // end if (nc == 1)
1448
1449 else if (nc > 1) {
1450 // deal with the blocks marked FIRST
1451 for (block = 0; block < ngrid; block ++) {
1452 if (mark[block] == FIRST) {
1453 ncb = nc*block;
1454 for (point = 0; point < nc; point ++) {
1455 vec_tmp[point] = b_val[ncb+point];
1456 }
1457 start_DATA = nc2*block;
1458 for (band = 0; band < nband; band ++) {
1459 width = offsets[band];
1460 column = block + width;
1461 if (width < 0) {
1462 if (column >= 0) {
1463 start_data = nc2*column;
1464 start_vecu = nc*column;
1465 blkcontr2( start_data, start_vecu, 0, nc,
1466 offdiag[band], u_val, vec_tmp );
1467 }
1468 }
1469 else { // width > 0
1470 if (column < ngrid) {
1471 start_vecu = nc*column;
1472 blkcontr2( start_DATA, start_vecu, 0, nc,
1473 offdiag[band], u_val, vec_tmp );
1474 }
1475 }
1476 } // end for band
1477
1478 // subblock smoothing
1479 aAxpby(weight, one_minus_weight, nc,
1480 diaginv+start_DATA, vec_tmp, u_val+nc*block);
1481 } // end if (mark[block] == FIRST)
1482
1483 } // end for block
1484
1485 // deal with the blocks marked SECOND
1486 for (block = 0; block < ngrid; block ++) {
1487 if (mark[block] == SECOND) {
1488 ncb = nc*block;
1489 for (point = 0; point < nc; point ++) {
1490 vec_tmp[point] = b_val[ncb+point];
1491 }
1492 start_DATA = nc2*block;
1493 for (band = 0; band < nband; band ++) {
1494 width = offsets[band];
1495 column = block + width;
1496 if (width < 0) {
1497 if (column >= 0) {
1498 start_data = nc2*column;
1499 start_vecu = nc*column;
1500 blkcontr2( start_data, start_vecu, 0, nc,
1501 offdiag[band], u_val, vec_tmp );
1502 }
1503 }
1504 else { // width > 0
1505 if (column < ngrid) {
1506 start_vecu = nc*column;
1507 blkcontr2( start_DATA, start_vecu, 0, nc,
1508 offdiag[band], u_val, vec_tmp );
1509 }
1510 }
1511 } // end for band
1512
1513 // subblock smoothing
1514 aAxpby(weight, one_minus_weight, nc,
1515 diaginv+start_DATA, vec_tmp, u_val+nc*block);
1516 } // end if (mark[block] == SECOND)
1517
1518 } // end for block
1519
1520 } // end else if (nc > 1)
1521 else {
1522 printf("### ERROR: nc is illegal! [%s:%d]\n", __FILE__, __LINE__);
1523 return;
1524 }
1525
1526 fasp_mem_free(vec_tmp); vec_tmp = NULL;
1527}
1528
1543void fasp_generate_diaginv_block (dSTRmat *A,
1544 ivector *neigh,
1545 dvector *diaginv,
1546 ivector *pivot)
1547{
1548 // information about A
1549 const INT nc = A->nc;
1550 const INT ngrid = A->ngrid;
1551 const INT nband = A->nband;
1552
1553 INT *offsets = A->offsets;
1554 REAL *diag = A->diag;
1555 REAL **offdiag = A->offdiag;
1556
1557 // information about neighbors
1558 INT nneigh;
1559 if (!neigh) {
1560 nneigh = 0;
1561 }
1562 else {
1563 nneigh= neigh->row/ngrid;
1564 }
1565
1566 // local variable
1567 INT i, j, k, l, m, n, nbd, p;
1568 INT count;
1569 INT block_size;
1570 INT mem_inv = 0;
1571 INT mem_pivot = 0;
1572
1573 // allocation
1574 REAL *temp = (REAL *)fasp_mem_calloc(((nneigh+1)*nc)*((nneigh+1)*nc)*ngrid, sizeof(REAL));
1575 INT *tmp = (INT *)fasp_mem_calloc(((nneigh+1)*nc)*ngrid, sizeof(INT));
1576
1577 // main loop
1578 for (i=0; i<ngrid; ++i) {
1579 // count number of neighbors of node i
1580 count = 1;
1581 for (l=0; l<nneigh; ++l) {
1582 if (neigh->val[i*nneigh+l] >= 0) count++ ;
1583 }
1584
1585 // prepare the inverse of diagonal block i
1586 block_size = count*nc;
1587
1588 diaginv[i].row = block_size*block_size;
1589 diaginv[i].val = temp + mem_inv;
1590 mem_inv += diaginv[i].row;
1591
1592 pivot[i].row = block_size;
1593 pivot[i].val = tmp + mem_pivot;
1594 mem_pivot += pivot[i].row;
1595
1596 // put the diagonal block corresponding to node i
1597 for (j=0; j<nc; ++j) {
1598 for (k=0; k<nc; ++k) {
1599 diaginv[i].val[j*block_size+k] = diag[i*nc*nc + j*nc + k];
1600 }
1601 }
1602
1603 // put the blocks corresponding to the neighbor of node i
1604 count = 1;
1605 for (l=0; l<nneigh; ++l) {
1606 p = neigh->val[i*nneigh+l];
1607 if (p >= 0){
1608 // put the diagonal block corresponding to this neighbor
1609 for (j=0; j<nc; ++j) {
1610 for (k=0; k<nc; ++k) {
1611 m = count*nc + j; n = count*nc+k;
1612 diaginv[i].val[m*block_size+n] = diag[p*nc*nc+j*nc+k];
1613 }
1614 }
1615
1616 for (nbd=0; nbd<nband; nbd++) {
1617 // put the block corresponding to (i, p)
1618 if ( offsets[nbd] == (p-i) ) {
1619 for (j=0; j<nc; ++j) {
1620 for(k=0; k<nc; ++k) {
1621 m = j; n = count*nc + k;
1622 diaginv[i].val[m*block_size+n] = offdiag[nbd][(p-MAX(p-i, 0))*nc*nc+j*nc+k];
1623 }
1624 }
1625 }
1626
1627 // put the block corresponding to (p, i)
1628 if ( offsets[nbd] == (i-p) ) {
1629 for (j=0; j<nc; ++j) {
1630 for(k=0; k<nc; ++k) {
1631 m = count*nc + j; n = k;
1632 diaginv[i].val[m*block_size+n] = offdiag[nbd][(i-MAX(i-p, 0))*nc*nc+j*nc+k];
1633 }
1634 }
1635 }
1636 }
1637 count++;
1638 } //end if
1639 } // end for (l=0; l<nneigh; ++l)
1640
1641 //fasp_smat_inv(diaginv[i].val, block_size);
1642 fasp_smat_lu_decomp(diaginv[i].val, pivot[i].val, block_size);
1643
1644 } // end of main loop
1645}
1646
1665void fasp_smoother_dstr_swz (dSTRmat *A,
1666 dvector *b,
1667 dvector *u,
1668 dvector *diaginv,
1669 ivector *pivot,
1670 ivector *neigh,
1671 ivector *order)
1672{
1673 // information about A
1674 const INT ngrid = A->ngrid;
1675 const INT nc = A->nc;
1676
1677 // information about neighbors
1678 INT nneigh;
1679 if (!neigh) {
1680 nneigh = 0;
1681 }
1682 else {
1683 nneigh= neigh->row/ngrid;
1684 }
1685
1686 // local variable
1687 INT i, j, k, l, p, ti;
1688
1689 // work space
1690 REAL *temp = (REAL *)fasp_mem_calloc(b->row + (nneigh+1)*nc + (nneigh+1)*nc, sizeof(REAL));
1691 dvector r, e, ri;
1692 r.row = b->row; r.val = temp;
1693 e.row = (nneigh+1)*nc; e.val = temp + b->row;
1694 ri.row = (nneigh+1)*nc; ri.val = temp + b->row + (nneigh+1)*nc;
1695
1696 // initial residual
1697 fasp_dvec_cp(b,&r);fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1698
1699 // main loop
1700 if (!order) {
1701 for (i=0; i<ngrid; ++i) {
1702 //-----------------------------------------------------
1703 // right hand side for A_ii e_i = r_i
1704 // rhs corresponding to node i
1705 for (j=0; j<nc; ++j) {
1706 ri.val[j] = r.val[i*nc + j];
1707 }
1708 // rhs corrsponding to the neighbors of node i
1709 k = 1;
1710 for (l=0; l<nneigh; ++l) {
1711 p=neigh->val[nneigh*i+l];
1712 if ( p>=0 ) {
1713 for (j=0; j<nc; ++j) {
1714 ri.val[k*nc+j] = r.val[p*nc+j];
1715 }
1716
1717 ++k;
1718 } // end if
1719 }
1720
1721 ri.row = k*nc;
1722 //----------------------------------------------------
1723 //----------------------------------------------------
1724 // solve local problem
1725 e.row = k*nc;
1726 //fasp_blas_smat_mxv(diaginv[ti].val, ri.val, k*nc, e.val);
1727 fasp_smat_lu_solve(diaginv[i].val, ri.val, pivot[i].val, e.val, k*nc);
1728 //----------------------------------------------------
1729 //----------------------------------------------------
1730 // update solution
1731 // solution corresponding to node i
1732 for (j=0; j<nc; ++j) {
1733 u->val[i*nc + j] += e.val[j];
1734 }
1735 // solution corresponding to the neighbor of node i
1736 k = 1;
1737 for (l=0; l<nneigh; ++l) {
1738 p=neigh->val[nneigh*i+l];
1739 if ( p>=0 ) {
1740 for (j=0; j<nc; ++j) {
1741 u->val[p*nc+j] += e.val[k*nc+j];
1742 }
1743
1744 ++k;
1745 } // end if
1746 }
1747 //----------------------------------------------------
1748 //----------------------------------------------------
1749 // update residule
1750 fasp_dvec_cp(b,&r); fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1751 }
1752 }
1753 else {
1754 for (i=0; i<ngrid; ++i) {
1755 ti = order->val[i];
1756 //-----------------------------------------------------
1757 // right hand side for A_ii e_i = r_i
1758 // rhs corresponding to node i
1759 for (j=0; j<nc; ++j) {
1760 ri.val[j] = r.val[ti*nc + j];
1761 }
1762 // rhs corrsponding to the neighbors of node i
1763 k = 1;
1764 for (l=0; l<nneigh; ++l) {
1765 p=neigh->val[nneigh*ti+l];
1766 if ( p>=0 ) {
1767 for (j=0; j<nc; ++j) {
1768 ri.val[k*nc+j] = r.val[p*nc+j];
1769 }
1770
1771 ++k;
1772 } // end if
1773 }
1774
1775 ri.row = k*nc;
1776 //----------------------------------------------------
1777 //----------------------------------------------------
1778 // solve local problem
1779 e.row = k*nc;
1780 //fasp_blas_smat_mxv(diaginv[ti].val, ri.val, k*nc, e.val);
1781 fasp_smat_lu_solve(diaginv[ti].val, ri.val, pivot[ti].val, e.val, k*nc);
1782 //----------------------------------------------------
1783 //----------------------------------------------------
1784 // update solution
1785 // solution corresponding to node i
1786 for (j=0; j<nc; ++j) {
1787 u->val[ti*nc + j] += e.val[j];
1788 }
1789 // solution corresponding to the neighbor of node i
1790 k = 1;
1791 for (l=0; l<nneigh; ++l) {
1792 p=neigh->val[nneigh*ti+l];
1793 if ( p>=0 ) {
1794 for (j=0; j<nc; ++j) {
1795 u->val[p*nc+j] += e.val[k*nc+j];
1796 }
1797
1798 ++k;
1799 } // end if
1800 }
1801 //----------------------------------------------------
1802 //----------------------------------------------------
1803 // update residule
1804 fasp_dvec_cp(b,&r); fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1805 }
1806 }// end of main loop
1807}
1808
1809
1810/*---------------------------------*/
1811/*-- Private Functions --*/
1812/*---------------------------------*/
1813
1833static void blkcontr2 (INT start_data,
1834 INT start_vecx,
1835 INT start_vecy,
1836 INT nc,
1837 REAL *data,
1838 REAL *x,
1839 REAL *y)
1840{
1841 INT i,j,k,m;
1842 if (start_vecy == 0) {
1843 for (i = 0; i < nc; i ++) {
1844 k = start_data + i*nc;
1845 for (j = 0; j < nc; j ++) {
1846 y[i] -= data[k+j]*x[start_vecx+j];
1847 }
1848 }
1849 }
1850 else {
1851 for (i = 0; i < nc; i ++) {
1852 k = start_data + i*nc;
1853 m = start_vecy + i;
1854 for (j = 0; j < nc; j ++) {
1855 y[m] -= data[k+j]*x[start_vecx+j];
1856 }
1857 }
1858 }
1859}
1860
1876static void aAxpby (REAL alpha,
1877 REAL beta,
1878 INT size,
1879 REAL *A,
1880 REAL *x,
1881 REAL *y)
1882{
1883 INT i,j;
1884 REAL tmp = 0.0;
1885 if (alpha == 0) {
1886 for (i = 0; i < size; i ++) {
1887 y[i] *= beta;
1888 }
1889 return;
1890 }
1891 tmp = beta / alpha;
1892 // y:=(beta/alpha)y
1893 for (i = 0; i < size; i ++) {
1894 y[i] *= tmp;
1895 }
1896 // y:=y+Ax
1897 for (i = 0; i < size; i ++) {
1898 for (j = 0; j < size; j ++) {
1899 y[i] += A[i*size+j]*x[j];
1900 }
1901 }
1902 // y:=alpha*y
1903 for (i = 0; i < size; i ++) {
1904 y[i] *= alpha;
1905 }
1906}
1907
1908/*---------------------------------*/
1909/*-- End of File --*/
1910/*---------------------------------*/
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_dvec_cp
void fasp_dvec_cp(const dvector *x, dvector *y)
Copy dvector x to dvector y.
Definition: AuxVector.c:334
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_lu_decomp
SHORT fasp_smat_lu_decomp(REAL *A, INT pivot[], const INT n)
LU decomposition of A using Doolittle's method.
Definition: BlaSmallMatLU.c:52
fasp_smat_lu_solve
SHORT fasp_smat_lu_solve(const REAL *A, REAL b[], const INT pivot[], REAL x[], const INT n)
Solving Ax=b using LU decomposition.
Definition: BlaSmallMatLU.c:124
fasp_blas_smat_mxv
void fasp_blas_smat_mxv(const REAL *a, const REAL *b, REAL *c, const INT n)
Compute the product of a small full matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:238
fasp_blas_dstr_aAxpy
void fasp_blas_dstr_aAxpy(const REAL alpha, const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvSTR.c:61
fasp_smoother_dstr_jacobi
void fasp_smoother_dstr_jacobi(dSTRmat *A, dvector *b, dvector *u)
Jacobi method as the smoother.
Definition: ItrSmootherSTR.c:43
fasp_smoother_dstr_sor_cf
void fasp_smoother_dstr_sor_cf(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark, const INT order, const REAL weight)
SOR method as the smoother in the C-F manner.
Definition: ItrSmootherSTR.c:1355
fasp_smoother_dstr_gs_cf
void fasp_smoother_dstr_gs_cf(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark, const INT order)
Gauss method as the smoother in the C-F manner.
Definition: ItrSmootherSTR.c:680
fasp_smoother_dstr_jacobi1
void fasp_smoother_dstr_jacobi1(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv)
Jacobi method as the smoother with diag_inv given.
Definition: ItrSmootherSTR.c:92
fasp_smoother_dstr_sor
void fasp_smoother_dstr_sor(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark, const REAL weight)
SOR method as the smoother.
Definition: ItrSmootherSTR.c:873
fasp_smoother_dstr_sor_ascend
void fasp_smoother_dstr_sor_ascend(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR method as the smoother in the ascending manner.
Definition: ItrSmootherSTR.c:981
fasp_smoother_dstr_gs1
void fasp_smoother_dstr_gs1(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark, REAL *diaginv)
Gauss-Seidel method as the smoother with diag_inv given.
Definition: ItrSmootherSTR.c:277
fasp_smoother_dstr_gs_descend
void fasp_smoother_dstr_gs_descend(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel method as the smoother in the descending manner.
Definition: ItrSmootherSTR.c:438
fasp_smoother_dstr_sor_order
void fasp_smoother_dstr_sor_order(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark, REAL weight)
SOR method as the smoother in the user-defined order.
Definition: ItrSmootherSTR.c:1224
fasp_smoother_dstr_gs_ascend
void fasp_smoother_dstr_gs_ascend(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel method as the smoother in the ascending manner.
Definition: ItrSmootherSTR.c:322
fasp_smoother_dstr_sor1
void fasp_smoother_dstr_sor1(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark, REAL *diaginv, const REAL weight)
SOR method as the smoother.
Definition: ItrSmootherSTR.c:935
fasp_smoother_dstr_gs_order
void fasp_smoother_dstr_gs_order(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark)
Gauss method as the smoother in the user-defined order.
Definition: ItrSmootherSTR.c:556
fasp_generate_diaginv_block
void fasp_generate_diaginv_block(dSTRmat *A, ivector *neigh, dvector *diaginv, ivector *pivot)
Generate inverse of diagonal block for block smoothers.
Definition: ItrSmootherSTR.c:1543
fasp_smoother_dstr_gs
void fasp_smoother_dstr_gs(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark)
Gauss-Seidel method as the smoother.
Definition: ItrSmootherSTR.c:217
fasp_smoother_dstr_sor_descend
void fasp_smoother_dstr_sor_descend(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR method as the smoother in the descending manner.
Definition: ItrSmootherSTR.c:1102
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
ASCEND
#define ASCEND
Definition: fasp_const.h:249
DESCEND
#define DESCEND
Definition: fasp_const.h:250
USERDEFINED
#define USERDEFINED
Type of ordering for smoothers.
Definition: fasp_const.h:246
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::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::offdiag
REAL ** offdiag
off-diagonal entries (dimension is nband * [(ngrid-|offsets|) * nc^2])
Definition: fasp.h:346
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
ivector::row
INT row
number of rows
Definition: fasp.h:371
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374