Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
smoother_str.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #include "fasp.h"
9 #include "fasp_functs.h"
10 
11 static void blkcontr2_str(INT start_data, INT start_vecx, INT start_vecy, INT nc,
12  REAL *data, REAL *x, REAL *y);
13 static void aAxpby(REAL alpha, REAL beta, INT size, REAL *A, REAL *x, REAL *y);
14 
15 /*---------------------------------*/
16 /*-- Public Functions --*/
17 /*---------------------------------*/
18 
32  dvector *b,
33  dvector *u)
34 {
35  INT nc = A->nc; // size of each block (number of components)
36  INT ngrid = A->ngrid; // number of grids
37  REAL *diag = A->diag; // Diagonal entries
38  REAL *diaginv = NULL; // Diagonal inverse, same size and storage scheme as A->diag
39 
40  INT nc2 = nc*nc;
41  INT size = nc2*ngrid;
42  INT block = 0;
43  INT start = 0;
44 
45  if (nc > 1) {
46  // allocate memory
47  diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
48 
49  // diaginv = diag;
50  fasp_array_cp(size,diag,diaginv);
51 
52  // generate diaginv
53  for (block = 0; block < ngrid; block ++) {
54  fasp_blas_smat_inv(diaginv+start, nc);
55  start += nc2;
56  }
57  }
58 
59  fasp_smoother_dstr_jacobi1(A, b, u, diaginv);
60 
61  fasp_mem_free(diaginv);
62 }
63 
64 
80  dvector *b,
81  dvector *u,
82  REAL *diaginv)
83 {
84  // information of A
85  INT ngrid = A->ngrid; // number of grids
86  INT nc = A->nc; // size of each block (number of components)
87  INT nband = A->nband; // number of off-diag band
88  INT *offsets = A->offsets; // offsets of the off-diagals
89  REAL *diag = A->diag; // Diagonal entries
90  REAL **offdiag = A->offdiag; // Off-diagonal entries
91 
92  // values of dvector b and u
93  REAL *b_val = b->val;
94  REAL *u_val = u->val;
95 
96  // local variables
97  INT block = 0;
98  INT point = 0;
99  INT band = 0;
100  INT width = 0;
101  INT size = nc*ngrid;
102  INT nc2 = nc*nc;
103  INT start = 0;
104  INT column = 0;
105  INT start_data = 0;
106  INT start_DATA = 0;
107  INT start_vecb = 0;
108  INT start_vecu = 0;
109 
110  // auxiliary array
111  REAL *b_tmp = NULL;
112 
113  // this should be done once and for all!!
114  b_tmp = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
115 
116  // b_tmp = b_val
117  fasp_array_cp(size,b_val,b_tmp);
118 
119  // It's not necessary to assign the smoothing order since the results doesn't depend on it
120  if (nc == 1) {
121  for (point = 0; point < ngrid; point ++) {
122  for (band = 0; band < nband; band ++) {
123  width = offsets[band];
124  column = point + width;
125  if (width < 0) {
126  if (column >= 0) {
127  b_tmp[point] -= offdiag[band][column]*u_val[column];
128  }
129  }
130  else {// width > 0
131  if (column < ngrid) {
132  b_tmp[point] -= offdiag[band][point]*u_val[column];
133  }
134  }
135  } // end for band
136  } // end for point
137 
138  for (point = 0; point < ngrid; point ++) {
139  // zero-diagonal should be tested previously
140  u_val[point] = b_tmp[point] / diag[point];
141  }
142  } // end if (nc == 1)
143  else if (nc > 1) {
144  for (block = 0; block < ngrid; block ++) {
145  start_DATA = nc2*block;
146  start_vecb = nc*block;
147  for (band = 0; band < nband; band ++) {
148  width = offsets[band];
149  column = block + width;
150  if (width < 0) {
151  if (column >= 0) {
152  start_data = nc2*column;
153  start_vecu = nc*column;
154  blkcontr2_str( start_data, start_vecu, start_vecb,
155  nc, offdiag[band], u_val, b_tmp );
156  }
157  }
158  else {// width > 0
159  if (column < ngrid) {
160  start_vecu = nc*column;
161  blkcontr2_str( start_DATA, start_vecu, start_vecb,
162  nc, offdiag[band], u_val, b_tmp );
163  }
164  }
165  } // end for band
166  } // end for block
167 
168  for (block = 0; block < ngrid; block ++) {
169  start = nc*block;
170  fasp_blas_smat_mxv(diaginv+nc2*block, b_tmp+start, u_val+start, nc);
171  }
172  } // end else if (nc > 1)
173  else {
174  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
175  return;
176  }
177  fasp_mem_free(b_tmp);
178 }
179 
204  dvector *b,
205  dvector *u,
206  const INT order,
207  INT *mark)
208 {
209  INT nc = A->nc; // size of each block (number of components)
210  INT ngrid = A->ngrid; // number of grids
211  REAL *diag = A->diag; // Diagonal entries
212  REAL *diaginv = NULL; // Diagonal inverse(when nc>1),same size and storage scheme as A->diag
213 
214  INT nc2 = nc*nc;
215  INT size = nc2*ngrid;
216  INT block = 0;
217  INT start = 0;
218 
219  if (nc > 1) {
220  // allocate memory
221  diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
222 
223  // diaginv = diag;
224  fasp_array_cp(size,diag,diaginv);
225 
226  // generate diaginv
227  for (block = 0; block < ngrid; block ++) {
228  fasp_blas_smat_inv(diaginv+start, nc);
229  start += nc2;
230  }
231  }
232 
233  fasp_smoother_dstr_gs1(A, b, u, order, mark, diaginv);
234 
235  fasp_mem_free(diaginv);
236 }
237 
264  dvector *b,
265  dvector *u,
266  const INT order,
267  INT *mark,
268  REAL *diaginv)
269 {
270 
271  if (!mark) {
272  if (order == ASCEND) // smooth ascendingly
273  {
274  fasp_smoother_dstr_gs_ascend(A, b, u, diaginv);
275  }
276  else if (order == DESCEND) // smooth descendingly
277  {
278  fasp_smoother_dstr_gs_descend(A, b, u, diaginv);
279  }
280  }
281  else {
282  if (order == USERDEFINED) // smooth according to the order 'mark' defined by user
283  {
284  fasp_smoother_dstr_gs_order(A, b, u, diaginv, mark);
285  }
286  else // smooth according to 'mark', where 'mark' is a CF_marker array
287  {
288  fasp_smoother_dstr_gs_cf(A, b, u, diaginv, mark, order);
289  }
290  }
291 }
292 
309  dvector *b,
310  dvector *u,
311  REAL *diaginv)
312 {
313  // information of A
314  INT ngrid = A->ngrid; // number of grids
315  INT nc = A->nc; // size of each block (number of components)
316  INT nband = A->nband; // number of off-diag band
317  INT *offsets = A->offsets; // offsets of the off-diagals
318  REAL *diag = A->diag; // Diagonal entries
319  REAL **offdiag = A->offdiag; // Off-diagonal entries
320 
321  // values of dvector b and u
322  REAL *b_val = b->val;
323  REAL *u_val = u->val;
324 
325  // local variables
326  INT block = 0;
327  INT point = 0;
328  INT band = 0;
329  INT width = 0;
330  INT nc2 = nc*nc;
331  INT ncb = 0;
332  INT column = 0;
333  INT start_data = 0;
334  INT start_DATA = 0;
335  INT start_vecu = 0;
336  REAL rhs = 0.0;
337 
338  // auxiliary array(nc*1 vector)
339  REAL *vec_tmp = NULL;
340 
341  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
342 
343  if (nc == 1) {
344  for (point = 0; point < ngrid; point ++) {
345  rhs = b_val[point];
346  for (band = 0; band < nband; band ++) {
347  width = offsets[band];
348  column = point + width;
349  if (width < 0) {
350  if (column >= 0) {
351  rhs -= offdiag[band][column]*u_val[column];
352  }
353  }
354  else { // width > 0
355  if (column < ngrid) {
356  rhs -= offdiag[band][point]*u_val[column];
357  }
358  }
359  } // end for band
360 
361  // zero-diagonal should be tested previously
362  u_val[point] = rhs / diag[point];
363 
364  } // end for point
365 
366  } // end if (nc == 1)
367 
368  else if (nc > 1) {
369  for (block = 0; block < ngrid; block ++) {
370  ncb = nc*block;
371  for (point = 0; point < nc; point ++) {
372  vec_tmp[point] = b_val[ncb+point];
373  }
374  start_DATA = nc2*block;
375  for (band = 0; band < nband; band ++) {
376  width = offsets[band];
377  column = block + width;
378  if (width < 0) {
379  if (column >= 0) {
380  start_data = nc2*column;
381  start_vecu = nc*column;
382  blkcontr2_str( start_data, start_vecu, 0, nc,
383  offdiag[band], u_val, vec_tmp );
384  }
385  }
386  else { // width > 0
387  if (column < ngrid) {
388  start_vecu = nc*column;
389  blkcontr2_str( start_DATA, start_vecu, 0, nc,
390  offdiag[band], u_val, vec_tmp );
391  }
392  }
393  } // end for band
394 
395  // subblock smoothing
396  fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
397 
398  } // end for block
399 
400  } // end else if (nc > 1)
401  else {
402  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
403  return;
404  }
405  fasp_mem_free(vec_tmp);
406 }
407 
424  dvector *b,
425  dvector *u,
426  REAL *diaginv)
427 {
428  // information of A
429  INT ngrid = A->ngrid; // number of grids
430  INT nc = A->nc; // size of each block (number of components)
431  INT nband = A->nband ; // number of off-diag band
432  INT *offsets = A->offsets; // offsets of the off-diagals
433  REAL *diag = A->diag; // Diagonal entries
434  REAL **offdiag = A->offdiag; // Off-diagonal entries
435 
436  // values of dvector b and u
437  REAL *b_val = b->val;
438  REAL *u_val = u->val;
439 
440  // local variables
441  INT block = 0;
442  INT point = 0;
443  INT band = 0;
444  INT width = 0;
445  INT nc2 = nc*nc;
446  INT ncb = 0;
447  INT column = 0;
448  INT start_data = 0;
449  INT start_DATA = 0;
450  INT start_vecu = 0;
451  REAL rhs = 0.0;
452 
453  // auxiliary array(nc*1 vector)
454  REAL *vec_tmp = NULL;
455 
456  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
457 
458  if (nc == 1) {
459  for (point = ngrid-1; point >= 0; point --) {
460  rhs = b_val[point];
461  for (band = 0; band < nband; band ++) {
462  width = offsets[band];
463  column = point + width;
464  if (width < 0) {
465  if (column >= 0) {
466  rhs -= offdiag[band][column]*u_val[column];
467  }
468  }
469  else { // width > 0
470  if (column < ngrid) {
471  rhs -= offdiag[band][point]*u_val[column];
472  }
473  }
474  } // end for band
475 
476  // zero-diagonal should be tested previously
477  u_val[point] = rhs / diag[point];
478 
479  } // end for point
480 
481  } // end if (nc == 1)
482 
483  else if (nc > 1) {
484  for (block = ngrid-1; block >= 0; block --) {
485  ncb = nc*block;
486  for (point = 0; point < nc; point ++) {
487  vec_tmp[point] = b_val[ncb+point];
488  }
489  start_DATA = nc2*block;
490  for (band = 0; band < nband; band ++) {
491  width = offsets[band];
492  column = block + width;
493  if (width < 0) {
494  if (column >= 0) {
495  start_data = nc2*column;
496  start_vecu = nc*column;
497  blkcontr2_str( start_data, start_vecu, 0, nc,
498  offdiag[band], u_val, vec_tmp );
499  }
500  }
501  else { // width > 0
502  if (column < ngrid) {
503  start_vecu = nc*column;
504  blkcontr2_str( start_DATA, start_vecu, 0, nc,
505  offdiag[band], u_val, vec_tmp );
506  }
507  }
508  } // end for band
509 
510  // subblock smoothing
511  fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
512 
513  } // end for block
514 
515  } // end else if (nc > 1)
516 
517  else {
518  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
519  return;
520  }
521  fasp_mem_free(vec_tmp);
522 }
523 
541  dvector *b,
542  dvector *u,
543  REAL *diaginv,
544  INT *mark)
545 {
546  // information of A
547  INT ngrid = A->ngrid; // number of grids
548  INT nc = A->nc; // size of each block (number of components)
549  INT nband = A->nband; // number of off-diag band
550  INT *offsets = A->offsets; // offsets of the off-diagals
551  REAL *diag = A->diag; // Diagonal entries
552  REAL **offdiag = A->offdiag; // Off-diagonal entries
553 
554  // values of dvector b and u
555  REAL *b_val = b->val;
556  REAL *u_val = u->val;
557 
558  // local variables
559  INT block = 0;
560  INT point = 0;
561  INT band = 0;
562  INT width = 0;
563  INT nc2 = nc*nc;
564  INT ncb = 0;
565  INT index = 0;
566  INT column = 0;
567  INT start_data = 0;
568  INT start_DATA = 0;
569  INT start_vecu = 0;
570  REAL rhs = 0.0;
571 
572  // auxiliary array(nc*1 vector)
573  REAL *vec_tmp = NULL;
574 
575  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
576 
577  if (nc == 1) {
578  for (index = 0; index < ngrid; index ++) {
579  point = mark[index];
580  rhs = b_val[point];
581  for (band = 0; band < nband; band ++) {
582  width = offsets[band];
583  column = point + width;
584  if (width < 0) {
585  if (column >= 0) {
586  rhs -= offdiag[band][column]*u_val[column];
587  }
588  }
589  else { // width > 0
590  if (column < ngrid) {
591  rhs -= offdiag[band][point]*u_val[column];
592  }
593  }
594  } // end for band
595 
596  // zero-diagonal should be tested previously
597  u_val[point] = rhs / diag[point];
598 
599  } // end for index
600 
601  } // end if (nc == 1)
602 
603  else if (nc > 1) {
604  for (index = 0; index < ngrid; index ++) {
605  block = mark[index];
606  ncb = nc*block;
607  for (point = 0; point < nc; point ++) {
608  vec_tmp[point] = b_val[ncb+point];
609  }
610  start_DATA = nc2*block;
611  for (band = 0; band < nband; band ++) {
612  width = offsets[band];
613  column = block + width;
614  if (width < 0) {
615  if (column >= 0) {
616  start_data = nc2*column;
617  start_vecu = nc*column;
618  blkcontr2_str( start_data, start_vecu, 0, nc,
619  offdiag[band], u_val, vec_tmp );
620  }
621  }
622  else { // width > 0
623  if (column < ngrid) {
624  start_vecu = nc*column;
625  blkcontr2_str( start_DATA, start_vecu, 0, nc,
626  offdiag[band], u_val, vec_tmp );
627  }
628  }
629  } // end for band
630 
631  // subblock smoothing
632  fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
633 
634  } // end for index
635 
636  } // end else if (nc > 1)
637  else {
638  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
639  return;
640  }
641  fasp_mem_free(vec_tmp);
642 }
643 
664  dvector *b,
665  dvector *u,
666  REAL *diaginv,
667  INT *mark,
668  const INT order)
669 {
670  // information of A
671  INT ngrid = A->ngrid; // number of grids
672  INT nc = A->nc; // size of each block (number of components)
673  INT nband = A->nband ; // number of off-diag band
674  INT *offsets = A->offsets; // offsets of the off-diagals
675  REAL *diag = A->diag; // Diagonal entries
676  REAL **offdiag = A->offdiag; // Off-diagonal entries
677 
678  // values of dvector b and u
679  REAL *b_val = b->val;
680  REAL *u_val = u->val;
681 
682  // local variables
683  INT block = 0;
684  INT point = 0;
685  INT band = 0;
686  INT width = 0;
687  INT nc2 = nc*nc;
688  INT ncb = 0;
689  INT column = 0;
690  INT start_data = 0;
691  INT start_DATA = 0;
692  INT start_vecu = 0;
693  INT FIRST = order; // which kind of points to be smoothed firstly?
694  INT SECOND = -order; // which kind of points to be smoothed secondly?
695 
696  REAL rhs = 0.0;
697 
698  // auxiliary array(nc*1 vector)
699  REAL *vec_tmp = NULL;
700 
701  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
702 
703  if (nc == 1) {
704  // deal with the points marked FIRST
705  for (point = 0; point < ngrid; point ++) {
706  if (mark[point] == FIRST) {
707  rhs = b_val[point];
708  for (band = 0; band < nband; band ++) {
709  width = offsets[band];
710  column = point + width;
711  if (width < 0) {
712  if (column >= 0) {
713  rhs -= offdiag[band][column]*u_val[column];
714  }
715  }
716  else { // width > 0
717  if (column < ngrid) {
718  rhs -= offdiag[band][point]*u_val[column];
719  }
720  }
721  } // end for band
722 
723  // zero-diagonal should be tested previously
724  u_val[point] = rhs / diag[point];
725  } // end if (mark[point] == FIRST)
726  } // end for point
727 
728  // deal with the points marked SECOND
729  for (point = 0; point < ngrid; point ++) {
730  if (mark[point] == SECOND) {
731  rhs = b_val[point];
732  for (band = 0; band < nband; band ++) {
733  width = offsets[band];
734  column = point + width;
735  if (width < 0) {
736  if (column >= 0) {
737  rhs -= offdiag[band][column]*u_val[column];
738  }
739  }
740  else { // width > 0
741  if (column < ngrid) {
742  rhs -= offdiag[band][point]*u_val[column];
743  }
744  }
745  } // end for band
746 
747  // zero-diagonal should be tested previously
748  u_val[point] = rhs / diag[point];
749  } // end if (mark[point] == SECOND)
750  } // end for point
751 
752  } // end if (nc == 1)
753 
754  else if (nc > 1) {
755  // deal with the blocks marked FIRST
756  for (block = 0; block < ngrid; block ++) {
757  if (mark[block] == FIRST) {
758  ncb = nc*block;
759  for (point = 0; point < nc; point ++) {
760  vec_tmp[point] = b_val[ncb+point];
761  }
762  start_DATA = nc2*block;
763  for (band = 0; band < nband; band ++) {
764  width = offsets[band];
765  column = block + width;
766  if (width < 0) {
767  if (column >= 0) {
768  start_data = nc2*column;
769  start_vecu = nc*column;
770  blkcontr2_str( start_data, start_vecu, 0, nc,
771  offdiag[band], u_val, vec_tmp );
772  }
773  }
774  else { // width > 0
775  if (column < ngrid) {
776  start_vecu = nc*column;
777  blkcontr2_str( start_DATA, start_vecu, 0, nc,
778  offdiag[band], u_val, vec_tmp );
779  }
780  }
781  } // end for band
782 
783  // subblock smoothing
784  fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
785  } // end if (mark[block] == FIRST)
786 
787  } // end for block
788 
789  // deal with the blocks marked SECOND
790  for (block = 0; block < ngrid; block ++) {
791  if (mark[block] == SECOND) {
792  ncb = nc*block;
793  for (point = 0; point < nc; point ++) {
794  vec_tmp[point] = b_val[ncb+point];
795  }
796  start_DATA = nc2*block;
797  for (band = 0; band < nband; band ++) {
798  width = offsets[band];
799  column = block + width;
800  if (width < 0) {
801  if (column >= 0) {
802  start_data = nc2*column;
803  start_vecu = nc*column;
804  blkcontr2_str( start_data, start_vecu, 0, nc,
805  offdiag[band], u_val, vec_tmp );
806  }
807  }
808  else { // width > 0
809  if (column < ngrid) {
810  start_vecu = nc*column;
811  blkcontr2_str( start_DATA, start_vecu, 0, nc,
812  offdiag[band], u_val, vec_tmp );
813  }
814  }
815  } // end for band
816 
817  // subblock smoothing
818  fasp_blas_smat_mxv(diaginv+start_DATA, vec_tmp, u_val+nc*block, nc);
819  } // end if (mark[block] == SECOND)
820 
821  } // end for block
822 
823  } // end else if (nc > 1)
824  else {
825  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
826  return;
827  }
828  fasp_mem_free(vec_tmp);
829 }
830 
856  dvector *b,
857  dvector *u,
858  const INT order,
859  INT *mark,
860  const REAL weight)
861 {
862  INT nc = A->nc; // size of each block (number of components)
863  INT ngrid = A->ngrid; // number of grids
864  REAL *diag = A->diag; // Diagonal entries
865  REAL *diaginv = NULL; // Diagonal inverse(when nc>1),same size and storage scheme as A->diag
866 
867  INT nc2 = nc*nc;
868  INT size = nc2*ngrid;
869  INT block = 0;
870  INT start = 0;
871 
872  if (nc > 1) {
873  // allocate memory
874  diaginv = (REAL *)fasp_mem_calloc(size,sizeof(REAL));
875 
876  // diaginv = diag;
877  fasp_array_cp(size,diag,diaginv);
878 
879  // generate diaginv
880  for (block = 0; block < ngrid; block ++) {
881  fasp_blas_smat_inv(diaginv+start, nc);
882  start += nc2;
883  }
884  }
885 
886  fasp_smoother_dstr_sor1(A, b, u, order, mark, diaginv, weight);
887 
888  fasp_mem_free(diaginv);
889 }
890 
917  dvector *b,
918  dvector *u,
919  const INT order,
920  INT *mark,
921  REAL *diaginv,
922  const REAL weight)
923 {
924  if (!mark) {
925  if (order == ASCEND) // smooth ascendingly
926  {
927  fasp_smoother_dstr_sor_ascend(A, b, u, diaginv, weight);
928  }
929  else if (order == DESCEND) // smooth descendingly
930  {
931  fasp_smoother_dstr_sor_descend(A, b, u, diaginv, weight);
932  }
933  }
934  else {
935  if (order == USERDEFINED) // smooth according to the order 'mark' defined by user
936  {
937  fasp_smoother_dstr_sor_order(A, b, u, diaginv, mark, weight);
938  }
939  else // smooth according to 'mark', where 'mark' is a CF_marker array
940  {
941  fasp_smoother_dstr_sor_cf(A, b, u, diaginv, mark, order, weight);
942  }
943  }
944 }
945 
963  dvector *b,
964  dvector *u,
965  REAL *diaginv,
966  REAL weight)
967 {
968  // information of A
969  INT ngrid = A->ngrid; // number of grids
970  INT nc = A->nc; // size of each block (number of components)
971  INT nband = A->nband ; // number of off-diag band
972  INT *offsets = A->offsets; // offsets of the off-diagals
973  REAL *diag = A->diag; // Diagonal entries
974  REAL **offdiag = A->offdiag; // Off-diagonal entries
975 
976  // values of dvector b and u
977  REAL *b_val = b->val;
978  REAL *u_val = u->val;
979 
980  // local variables
981  INT block = 0;
982  INT point = 0;
983  INT band = 0;
984  INT width = 0;
985  INT nc2 = nc*nc;
986  INT ncb = 0;
987  INT column = 0;
988  INT start_data = 0;
989  INT start_DATA = 0;
990  INT start_vecu = 0;
991  REAL rhs = 0.0;
992  REAL one_minus_weight = 1.0 - weight;
993 
994  // auxiliary array(nc*1 vector)
995  REAL *vec_tmp = NULL;
996 
997  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
998 
999  if (nc == 1) {
1000  for (point = 0; point < ngrid; point ++) {
1001  rhs = b_val[point];
1002  for (band = 0; band < nband; band ++) {
1003  width = offsets[band];
1004  column = point + width;
1005  if (width < 0) {
1006  if (column >= 0) {
1007  rhs -= offdiag[band][column]*u_val[column];
1008  }
1009  }
1010  else { // width > 0
1011  if (column < ngrid) {
1012  rhs -= offdiag[band][point]*u_val[column];
1013  }
1014  }
1015  } // end for band
1016 
1017  // zero-diagonal should be tested previously
1018  u_val[point] = one_minus_weight*u_val[point] +
1019  weight*(rhs / diag[point]);
1020 
1021  } // end for point
1022 
1023  } // end if (nc == 1)
1024 
1025  else if (nc > 1) {
1026  for (block = 0; block < ngrid; block ++) {
1027  ncb = nc*block;
1028  for (point = 0; point < nc; point ++) {
1029  vec_tmp[point] = b_val[ncb+point];
1030  }
1031  start_DATA = nc2*block;
1032  for (band = 0; band < nband; band ++) {
1033  width = offsets[band];
1034  column = block + width;
1035  if (width < 0) {
1036  if (column >= 0) {
1037  start_data = nc2*column;
1038  start_vecu = nc*column;
1039  blkcontr2_str( start_data, start_vecu, 0, nc,
1040  offdiag[band], u_val, vec_tmp );
1041  }
1042  }
1043  else { // width > 0
1044  if (column < ngrid) {
1045  start_vecu = nc*column;
1046  blkcontr2_str( start_DATA, start_vecu, 0, nc,
1047  offdiag[band], u_val, vec_tmp );
1048  }
1049  }
1050  } // end for band
1051 
1052  // subblock smoothing
1053  aAxpby(weight, one_minus_weight, nc,
1054  diaginv+start_DATA, vec_tmp, u_val+nc*block);
1055 
1056  } // end for block
1057 
1058  } // end else if (nc > 1)
1059  else {
1060  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
1061  return;
1062  }
1063  fasp_mem_free(vec_tmp);
1064 }
1065 
1083  dvector *b,
1084  dvector *u,
1085  REAL *diaginv,
1086  REAL weight)
1087 {
1088  // information of A
1089  INT ngrid = A->ngrid; // number of grids
1090  INT nc = A->nc; // size of each block (number of components)
1091  INT nband = A->nband ; // number of off-diag band
1092  INT *offsets = A->offsets; // offsets of the off-diagals
1093  REAL *diag = A->diag; // Diagonal entries
1094  REAL **offdiag = A->offdiag; // Off-diagonal entries
1095 
1096  // values of dvector b and u
1097  REAL *b_val = b->val;
1098  REAL *u_val = u->val;
1099 
1100  // local variables
1101  INT block = 0;
1102  INT point = 0;
1103  INT band = 0;
1104  INT width = 0;
1105  INT nc2 = nc*nc;
1106  INT ncb = 0;
1107  INT column = 0;
1108  INT start_data = 0;
1109  INT start_DATA = 0;
1110  INT start_vecu = 0;
1111  REAL rhs = 0.0;
1112  REAL one_minus_weight = 1.0 - weight;
1113 
1114  // auxiliary array(nc*1 vector)
1115  REAL *vec_tmp = NULL;
1116 
1117  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1118 
1119  if (nc == 1) {
1120  for (point = ngrid-1; point >= 0; point --) {
1121  rhs = b_val[point];
1122  for (band = 0; band < nband; band ++) {
1123  width = offsets[band];
1124  column = point + width;
1125  if (width < 0) {
1126  if (column >= 0) {
1127  rhs -= offdiag[band][column]*u_val[column];
1128  }
1129  }
1130  else { // width > 0
1131  if (column < ngrid) {
1132  rhs -= offdiag[band][point]*u_val[column];
1133  }
1134  }
1135  } // end for band
1136 
1137  // zero-diagonal should be tested previously
1138  u_val[point] = one_minus_weight*u_val[point] +
1139  weight*(rhs / diag[point]);
1140 
1141  } // end for point
1142 
1143  } // end if (nc == 1)
1144 
1145  else if (nc > 1) {
1146  for (block = ngrid-1; block >= 0; block --) {
1147  ncb = nc*block;
1148  for (point = 0; point < nc; point ++) {
1149  vec_tmp[point] = b_val[ncb+point];
1150  }
1151  start_DATA = nc2*block;
1152  for (band = 0; band < nband; band ++) {
1153  width = offsets[band];
1154  column = block + width;
1155  if (width < 0) {
1156  if (column >= 0) {
1157  start_data = nc2*column;
1158  start_vecu = nc*column;
1159  blkcontr2_str( start_data, start_vecu, 0, nc,
1160  offdiag[band], u_val, vec_tmp );
1161  }
1162  }
1163  else { // width > 0
1164  if (column < ngrid) {
1165  start_vecu = nc*column;
1166  blkcontr2_str( start_DATA, start_vecu, 0, nc,
1167  offdiag[band], u_val, vec_tmp );
1168  }
1169  }
1170  } // end for band
1171 
1172  // subblock smoothing
1173  aAxpby(weight, one_minus_weight, nc,
1174  diaginv+start_DATA, vec_tmp, u_val+nc*block);
1175 
1176  } // end for block
1177 
1178  } // end else if (nc > 1)
1179  else {
1180  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
1181  return;
1182  }
1183  fasp_mem_free(vec_tmp);
1184 }
1185 
1204  dvector *b,
1205  dvector *u,
1206  REAL *diaginv,
1207  INT *mark,
1208  REAL weight)
1209 {
1210  // information of A
1211  INT ngrid = A->ngrid; // number of grids
1212  INT nc = A->nc; // size of each block (number of components)
1213  INT nband = A->nband ; // number of off-diag band
1214  INT *offsets = A->offsets; // offsets of the off-diagals
1215  REAL *diag = A->diag; // Diagonal entries
1216  REAL **offdiag = A->offdiag; // Off-diagonal entries
1217 
1218  // values of dvector b and u
1219  REAL *b_val = b->val;
1220  REAL *u_val = u->val;
1221 
1222  // local variables
1223  INT block = 0;
1224  INT point = 0;
1225  INT band = 0;
1226  INT width = 0;
1227  INT nc2 = nc*nc;
1228  INT ncb = 0;
1229  INT column = 0;
1230  INT index = 0;
1231  INT start_data = 0;
1232  INT start_DATA = 0;
1233  INT start_vecu = 0;
1234  REAL rhs = 0.0;
1235  REAL one_minus_weight = 1.0 - weight;
1236 
1237  // auxiliary array(nc*1 vector)
1238  REAL *vec_tmp = NULL;
1239 
1240  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1241 
1242  if (nc == 1) {
1243  for (index = 0; index < ngrid; index ++) {
1244  point = mark[index];
1245  rhs = b_val[point];
1246  for (band = 0; band < nband; band ++) {
1247  width = offsets[band];
1248  column = point + width;
1249  if (width < 0) {
1250  if (column >= 0) {
1251  rhs -= offdiag[band][column]*u_val[column];
1252  }
1253  }
1254  else { // width > 0
1255  if (column < ngrid) {
1256  rhs -= offdiag[band][point]*u_val[column];
1257  }
1258  }
1259  } // end for band
1260 
1261  // zero-diagonal should be tested previously
1262  u_val[point] = one_minus_weight*u_val[point] +
1263  weight*(rhs / diag[point]);
1264 
1265  } // end for index
1266 
1267  } // end if (nc == 1)
1268 
1269  else if (nc > 1) {
1270  for (index = 0; index < ngrid; index ++) {
1271  block = mark[index];
1272  ncb = nc*block;
1273  for (point = 0; point < nc; point ++) {
1274  vec_tmp[point] = b_val[ncb+point];
1275  }
1276  start_DATA = nc2*block;
1277  for (band = 0; band < nband; band ++) {
1278  width = offsets[band];
1279  column = block + width;
1280  if (width < 0) {
1281  if (column >= 0) {
1282  start_data = nc2*column;
1283  start_vecu = nc*column;
1284  blkcontr2_str( start_data, start_vecu, 0, nc,
1285  offdiag[band], u_val, vec_tmp );
1286  }
1287  }
1288  else { // width > 0
1289  if (column < ngrid) {
1290  start_vecu = nc*column;
1291  blkcontr2_str( start_DATA, start_vecu, 0, nc,
1292  offdiag[band], u_val, vec_tmp );
1293  }
1294  }
1295  } // end for band
1296 
1297  // subblock smoothing
1298  aAxpby(weight, one_minus_weight, nc,
1299  diaginv+start_DATA, vec_tmp, u_val+nc*block);
1300 
1301  } // end for index
1302 
1303  } // end else if (nc > 1)
1304 
1305  else {
1306  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
1307  return;
1308  }
1309 
1310  fasp_mem_free(vec_tmp);
1311 }
1312 
1335  dvector *b,
1336  dvector *u,
1337  REAL *diaginv,
1338  INT *mark,
1339  const INT order,
1340  const REAL weight)
1341 {
1342  // information of A
1343  INT ngrid = A->ngrid; // number of grids
1344  INT nc = A->nc; // size of each block (number of components)
1345  INT nband = A->nband ; // number of off-diag band
1346  INT *offsets = A->offsets; // offsets of the off-diagals
1347  REAL *diag = A->diag; // Diagonal entries
1348  REAL **offdiag = A->offdiag; // Off-diagonal entries
1349 
1350  // values of dvector b and u
1351  REAL *b_val = b->val;
1352  REAL *u_val = u->val;
1353 
1354  // local variables
1355  INT block = 0;
1356  INT point = 0;
1357  INT band = 0;
1358  INT width = 0;
1359  INT nc2 = nc*nc;
1360  INT ncb = 0;
1361  INT column = 0;
1362  INT start_data = 0;
1363  INT start_DATA = 0;
1364  INT start_vecu = 0;
1365  REAL rhs = 0.0;
1366  REAL one_minus_weight = 1.0 - weight;
1367  INT FIRST = order; // which kind of points to be smoothed firstly?
1368  INT SECOND = -order; // which kind of points to be smoothed secondly?
1369 
1370  // auxiliary array(nc*1 vector)
1371  REAL *vec_tmp = NULL;
1372 
1373  vec_tmp = (REAL *)fasp_mem_calloc(nc,sizeof(REAL));
1374 
1375  if (nc == 1) {
1376  // deal with the points marked FIRST
1377  for (point = 0; point < ngrid; point ++) {
1378  if (mark[point] == FIRST) {
1379  rhs = b_val[point];
1380  for (band = 0; band < nband; band ++) {
1381  width = offsets[band];
1382  column = point + width;
1383  if (width < 0) {
1384  if (column >= 0) {
1385  rhs -= offdiag[band][column]*u_val[column];
1386  }
1387  }
1388  else { // width > 0
1389  if (column < ngrid) {
1390  rhs -= offdiag[band][point]*u_val[column];
1391  }
1392  }
1393  } // end for band
1394 
1395  // zero-diagonal should be tested previously
1396  u_val[point] = one_minus_weight*u_val[point] +
1397  weight*(rhs / diag[point]);
1398 
1399  } // end if (mark[point] == FIRST)
1400  } // end for point
1401 
1402  // deal with the points marked SECOND
1403  for (point = 0; point < ngrid; point ++) {
1404  if (mark[point] == SECOND) {
1405  rhs = b_val[point];
1406  for (band = 0; band < nband; band ++) {
1407  width = offsets[band];
1408  column = point + width;
1409  if (width < 0) {
1410  if (column >= 0) {
1411  rhs -= offdiag[band][column]*u_val[column];
1412  }
1413  }
1414  else { // width > 0
1415  if (column < ngrid) {
1416  rhs -= offdiag[band][point]*u_val[column];
1417  }
1418  }
1419  } // end for band
1420 
1421  // zero-diagonal should be tested previously
1422  u_val[point] = rhs / diag[point];
1423  } // end if (mark[point] == SECOND)
1424  } // end for point
1425 
1426  } // end if (nc == 1)
1427 
1428  else if (nc > 1) {
1429  // deal with the blocks marked FIRST
1430  for (block = 0; block < ngrid; block ++) {
1431  if (mark[block] == FIRST) {
1432  ncb = nc*block;
1433  for (point = 0; point < nc; point ++) {
1434  vec_tmp[point] = b_val[ncb+point];
1435  }
1436  start_DATA = nc2*block;
1437  for (band = 0; band < nband; band ++) {
1438  width = offsets[band];
1439  column = block + width;
1440  if (width < 0) {
1441  if (column >= 0) {
1442  start_data = nc2*column;
1443  start_vecu = nc*column;
1444  blkcontr2_str( start_data, start_vecu, 0, nc,
1445  offdiag[band], u_val, vec_tmp );
1446  }
1447  }
1448  else { // width > 0
1449  if (column < ngrid) {
1450  start_vecu = nc*column;
1451  blkcontr2_str( start_DATA, start_vecu, 0, nc,
1452  offdiag[band], u_val, vec_tmp );
1453  }
1454  }
1455  } // end for band
1456 
1457  // subblock smoothing
1458  aAxpby(weight, one_minus_weight, nc,
1459  diaginv+start_DATA, vec_tmp, u_val+nc*block);
1460  } // end if (mark[block] == FIRST)
1461 
1462  } // end for block
1463 
1464  // deal with the blocks marked SECOND
1465  for (block = 0; block < ngrid; block ++) {
1466  if (mark[block] == SECOND) {
1467  ncb = nc*block;
1468  for (point = 0; point < nc; point ++) {
1469  vec_tmp[point] = b_val[ncb+point];
1470  }
1471  start_DATA = nc2*block;
1472  for (band = 0; band < nband; band ++) {
1473  width = offsets[band];
1474  column = block + width;
1475  if (width < 0) {
1476  if (column >= 0) {
1477  start_data = nc2*column;
1478  start_vecu = nc*column;
1479  blkcontr2_str( start_data, start_vecu, 0, nc,
1480  offdiag[band], u_val, vec_tmp );
1481  }
1482  }
1483  else { // width > 0
1484  if (column < ngrid) {
1485  start_vecu = nc*column;
1486  blkcontr2_str( start_DATA, start_vecu, 0, nc,
1487  offdiag[band], u_val, vec_tmp );
1488  }
1489  }
1490  } // end for band
1491 
1492  // subblock smoothing
1493  aAxpby(weight, one_minus_weight, nc,
1494  diaginv+start_DATA, vec_tmp, u_val+nc*block);
1495  } // end if (mark[block] == SECOND)
1496 
1497  } // end for block
1498 
1499  } // end else if (nc > 1)
1500  else {
1501  printf("### ERROR: nc is illegal! %s : %d\n", __FILE__, __LINE__);
1502  return;
1503  }
1504  fasp_mem_free(vec_tmp);
1505 }
1506 
1522  ivector *neigh,
1523  dvector *diaginv,
1524  ivector *pivot)
1525 {
1526  // information about A
1527  const INT nc = A->nc;
1528  const INT ngrid = A->ngrid;
1529  const INT nband = A->nband;
1530 
1531  INT *offsets = A->offsets;
1532  REAL *diag = A->diag;
1533  REAL **offdiag = A->offdiag;
1534 
1535  // information about neighbors
1536  INT nneigh;
1537  if (!neigh) {
1538  nneigh = 0;
1539  }
1540  else {
1541  nneigh= neigh->row/ngrid;
1542  }
1543 
1544  // local variable
1545  INT i, j, k, l, m, n, nbd, p;
1546  INT count;
1547  INT block_size;
1548  INT mem_inv = 0;
1549  INT mem_pivot = 0;
1550 
1551  // allocation
1552  REAL *temp = (REAL *)fasp_mem_calloc(((nneigh+1)*nc)*((nneigh+1)*nc)*ngrid, sizeof(REAL));
1553  INT *tmp = (INT *)fasp_mem_calloc(((nneigh+1)*nc)*ngrid, sizeof(INT));
1554 
1555  // main loop
1556  for (i=0; i<ngrid; ++i) {
1557  // count number of neighbors of node i
1558  count = 1;
1559  for (l=0; l<nneigh; ++l) {
1560  if (neigh->val[i*nneigh+l] >= 0) count++ ;
1561  }
1562 
1563  // prepare the inverse of diagonal block i
1564  block_size = count*nc;
1565 
1566  diaginv[i].row = block_size*block_size;
1567  diaginv[i].val = temp + mem_inv;
1568  mem_inv += diaginv[i].row;
1569 
1570  pivot[i].row = block_size;
1571  pivot[i].val = tmp + mem_pivot;
1572  mem_pivot += pivot[i].row;
1573 
1574  // put the diagonal block corresponding to node i
1575  for (j=0; j<nc; ++j) {
1576  for (k=0; k<nc; ++k) {
1577  diaginv[i].val[j*block_size+k] = diag[i*nc*nc + j*nc + k];
1578  }
1579  }
1580 
1581  // put the blocks corresponding to the neighbor of node i
1582  count = 1;
1583  for (l=0; l<nneigh; ++l) {
1584  p = neigh->val[i*nneigh+l];
1585  if (p >= 0){
1586  // put the diagonal block corresponding to this neighbor
1587  for (j=0; j<nc; ++j) {
1588  for (k=0; k<nc; ++k) {
1589  m = count*nc + j; n = count*nc+k;
1590  diaginv[i].val[m*block_size+n] = diag[p*nc*nc+j*nc+k];
1591  }
1592  }
1593 
1594  for (nbd=0; nbd<nband; nbd++) {
1595  // put the block corresponding to (i, p)
1596  if ( offsets[nbd] == (p-i) ) {
1597  for (j=0; j<nc; ++j) {
1598  for(k=0; k<nc; ++k) {
1599  m = j; n = count*nc + k;
1600  diaginv[i].val[m*block_size+n] = offdiag[nbd][(p-MAX(p-i, 0))*nc*nc+j*nc+k];
1601  }
1602  }
1603  }
1604 
1605  // put the block corresponding to (p, i)
1606  if ( offsets[nbd] == (i-p) ) {
1607  for (j=0; j<nc; ++j) {
1608  for(k=0; k<nc; ++k) {
1609  m = count*nc + j; n = k;
1610  diaginv[i].val[m*block_size+n] = offdiag[nbd][(i-MAX(i-p, 0))*nc*nc+j*nc+k];
1611  }
1612  }
1613  }
1614  }
1615  count++;
1616  } //end if
1617  } // end for (l=0; l<nneigh; ++l)
1618 
1619  //fasp_blas_smat_inv(diaginv[i].val, block_size);
1620  fasp_smat_lu_decomp(diaginv[i].val, pivot[i].val, block_size);
1621 
1622  } // end of main loop
1623 }
1624 
1644  dvector *b,
1645  dvector *u,
1646  dvector *diaginv,
1647  ivector *pivot,
1648  ivector *neigh,
1649  ivector *order)
1650 {
1651  // information about A
1652  const INT ngrid = A->ngrid;
1653  const INT nc = A->nc;
1654 
1655  // information about neighbors
1656  INT nneigh;
1657  if (!neigh) {
1658  nneigh = 0;
1659  }
1660  else {
1661  nneigh= neigh->row/ngrid;
1662  }
1663 
1664  // local variable
1665  INT i, j, k, l, p, ti;
1666 
1667  // work space
1668  REAL *temp = (REAL *)fasp_mem_calloc(b->row + (nneigh+1)*nc + (nneigh+1)*nc, sizeof(REAL));
1669  dvector r, e, ri;
1670  r.row = b->row; r.val = temp;
1671  e.row = (nneigh+1)*nc; e.val = temp + b->row;
1672  ri.row = (nneigh+1)*nc; ri.val = temp + b->row + (nneigh+1)*nc;
1673 
1674  // initial residual
1675  fasp_dvec_cp(b,&r);fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1676 
1677  // main loop
1678  if (!order) {
1679  for (i=0; i<ngrid; ++i) {
1680  //-----------------------------------------------------
1681  // right hand side for A_ii e_i = r_i
1682  // rhs corresponding to node i
1683  for (j=0; j<nc; ++j) {
1684  ri.val[j] = r.val[i*nc + j];
1685  }
1686  // rhs corrsponding to the neighbors of node i
1687  k = 1;
1688  for (l=0; l<nneigh; ++l) {
1689  p=neigh->val[nneigh*i+l];
1690  if ( p>=0 ) {
1691  for (j=0; j<nc; ++j) {
1692  ri.val[k*nc+j] = r.val[p*nc+j];
1693  }
1694 
1695  ++k;
1696  } // end if
1697  }
1698 
1699  ri.row = k*nc;
1700  //----------------------------------------------------
1701  //----------------------------------------------------
1702  // solve local problem
1703  e.row = k*nc;
1704  //fasp_blas_smat_mxv(diaginv[ti].val, ri.val, k*nc, e.val);
1705  fasp_smat_lu_solve(diaginv[i].val, ri.val, pivot[i].val, e.val, k*nc);
1706  //----------------------------------------------------
1707  //----------------------------------------------------
1708  // update solution
1709  // solution corresponding to node i
1710  for (j=0; j<nc; ++j) {
1711  u->val[i*nc + j] += e.val[j];
1712  }
1713  // solution corresponding to the neighbor of node i
1714  k = 1;
1715  for (l=0; l<nneigh; ++l) {
1716  p=neigh->val[nneigh*i+l];
1717  if ( p>=0 ) {
1718  for (j=0; j<nc; ++j) {
1719  u->val[p*nc+j] += e.val[k*nc+j];
1720  }
1721 
1722  ++k;
1723  } // end if
1724  }
1725  //----------------------------------------------------
1726  //----------------------------------------------------
1727  // update residule
1728  fasp_dvec_cp(b,&r); fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1729  }
1730  }
1731  else {
1732  for (i=0; i<ngrid; ++i) {
1733  ti = order->val[i];
1734  //-----------------------------------------------------
1735  // right hand side for A_ii e_i = r_i
1736  // rhs corresponding to node i
1737  for (j=0; j<nc; ++j) {
1738  ri.val[j] = r.val[ti*nc + j];
1739  }
1740  // rhs corrsponding to the neighbors of node i
1741  k = 1;
1742  for (l=0; l<nneigh; ++l) {
1743  p=neigh->val[nneigh*ti+l];
1744  if ( p>=0 ) {
1745  for (j=0; j<nc; ++j) {
1746  ri.val[k*nc+j] = r.val[p*nc+j];
1747  }
1748 
1749  ++k;
1750  } // end if
1751  }
1752 
1753  ri.row = k*nc;
1754  //----------------------------------------------------
1755  //----------------------------------------------------
1756  // solve local problem
1757  e.row = k*nc;
1758  //fasp_blas_smat_mxv(diaginv[ti].val, ri.val, k*nc, e.val);
1759  fasp_smat_lu_solve(diaginv[ti].val, ri.val, pivot[ti].val, e.val, k*nc);
1760  //----------------------------------------------------
1761  //----------------------------------------------------
1762  // update solution
1763  // solution corresponding to node i
1764  for (j=0; j<nc; ++j) {
1765  u->val[ti*nc + j] += e.val[j];
1766  }
1767  // solution corresponding to the neighbor of node i
1768  k = 1;
1769  for (l=0; l<nneigh; ++l) {
1770  p=neigh->val[nneigh*ti+l];
1771  if ( p>=0 ) {
1772  for (j=0; j<nc; ++j) {
1773  u->val[p*nc+j] += e.val[k*nc+j];
1774  }
1775 
1776  ++k;
1777  } // end if
1778  }
1779  //----------------------------------------------------
1780  //----------------------------------------------------
1781  // update residule
1782  fasp_dvec_cp(b,&r); fasp_blas_dstr_aAxpy(-1.0,A,u->val,r.val);
1783  }
1784  }// end of main loop
1785 }
1786 
1787 
1788 /*---------------------------------*/
1789 /*-- Private Functions --*/
1790 /*---------------------------------*/
1791 
1811 static void blkcontr2_str (INT start_data, INT start_vecx, INT start_vecy, INT nc,
1812  REAL *data, REAL *x, REAL *y)
1813 {
1814  INT i,j,k,m;
1815  if (start_vecy == 0) {
1816  for (i = 0; i < nc; i ++) {
1817  k = start_data + i*nc;
1818  for (j = 0; j < nc; j ++) {
1819  y[i] -= data[k+j]*x[start_vecx+j];
1820  }
1821  }
1822  }
1823  else {
1824  for (i = 0; i < nc; i ++) {
1825  k = start_data + i*nc;
1826  m = start_vecy + i;
1827  for (j = 0; j < nc; j ++) {
1828  y[m] -= data[k+j]*x[start_vecx+j];
1829  }
1830  }
1831  }
1832 }
1833 
1849 static void aAxpby (REAL alpha, REAL beta, INT size, REAL *A, REAL *x, REAL *y)
1850 {
1851  INT i,j;
1852  REAL tmp = 0.0;
1853  if (alpha == 0) {
1854  for (i = 0; i < size; i ++) {
1855  y[i] *= beta;
1856  }
1857  return;
1858  }
1859  tmp = beta / alpha;
1860  // y:=(beta/alpha)y
1861  for (i = 0; i < size; i ++) {
1862  y[i] *= tmp;
1863  }
1864  // y:=y+Ax
1865  for (i = 0; i < size; i ++) {
1866  for (j = 0; j < size; j ++) {
1867  y[i] += A[i*size+j]*x[j];
1868  }
1869  }
1870  // y:=alpha*y
1871  for (i = 0; i < size; i ++) {
1872  y[i] *= alpha;
1873  }
1874 }
1875 
1876 /*---------------------------------*/
1877 /*-- End of File --*/
1878 /*---------------------------------*/
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.
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
void fasp_blas_dstr_aAxpy(const REAL alpha, dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_str.c:47
void fasp_generate_diaginv_block(dSTRmat *A, ivector *neigh, dvector *diaginv, ivector *pivot)
Generate inverse of diagonal block for block smoothers.
INT ngrid
number of grids
Definition: fasp.h:322
INT fasp_blas_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: smat.c:554
INT nc
size of each block (number of components)
Definition: fasp.h:319
void fasp_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: smoother_str.c:263
REAL ** offdiag
off-diagonal entries (dimension is nband * [(ngrid-|offsets|) * nc^2])
Definition: fasp.h:334
INT * offsets
offsets of the off-diagonals (length is nband)
Definition: fasp.h:331
Structure matrix of REAL type.
Definition: fasp.h:304
void fasp_dvec_cp(dvector *x, dvector *y)
Copy dvector x to dvector y.
Definition: vec.c:345
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void fasp_smoother_dstr_gs(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark)
Gauss-Seidel method as the smoother.
Definition: smoother_str.c:203
Vector with n entries of REAL type.
Definition: fasp.h:342
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: smoother_str.c:962
#define INT
Definition: fasp.h:64
#define DESCEND
Definition: fasp_const.h:232
INT nband
number of off-diag bands
Definition: fasp.h:328
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.
void fasp_smoother_dstr_schwarz(dSTRmat *A, dvector *b, dvector *u, dvector *diaginv, ivector *pivot, ivector *neigh, ivector *order)
Schwarz method as the smoother.
#define ASCEND
Definition: fasp_const.h:231
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
Main header file for FASP.
SHORT fasp_smat_lu_decomp(REAL *A, INT pivot[], const INT n)
LU decomposition of A usind Doolittle's method.
Definition: lu.c:46
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: smoother_str.c:423
INT * val
actual vector entries
Definition: fasp.h:362
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: smoother_str.c:663
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: smoother_str.c:916
INT row
number of rows
Definition: fasp.h:345
REAL * diag
diagonal entries (length is ngrid*(nc^2))
Definition: fasp.h:325
#define USERDEFINED
Type of ordering for smoothers.
Definition: fasp_const.h:228
INT count
INT row
number of rows
Definition: fasp.h:359
SHORT fasp_smat_lu_solve(REAL *A, REAL b[], INT pivot[], REAL x[], const INT n)
Solving Ax=b using LU decomposition.
Definition: lu.c:117
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: smoother_str.c:308
void fasp_smoother_dstr_jacobi(dSTRmat *A, dvector *b, dvector *u)
Jacobi method as the smoother.
Definition: smoother_str.c:31
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
void fasp_smoother_dstr_sor(dSTRmat *A, dvector *b, dvector *u, const INT order, INT *mark, const REAL weight)
SOR method as the smoother.
Definition: smoother_str.c:855
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.
void fasp_smoother_dstr_jacobi1(dSTRmat *A, dvector *b, dvector *u, REAL *diaginv)
Jacobi method as the smoother with diag_inv given.
Definition: smoother_str.c:79
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: smoother_str.c:540
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
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72