Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreAMGInterp.c
Go to the documentation of this file.
1
21#include <math.h>
22#include <time.h>
23
24#ifdef _OPENMP
25#include <omp.h>
26#endif
27
28#include "fasp.h"
29#include "fasp_functs.h"
30
31/*---------------------------------*/
32/*-- Declare Private Functions --*/
33/*---------------------------------*/
34
35static void interp_DIR(dCSRmat*, ivector*, dCSRmat*, AMG_param*);
36static void interp_RDC(dCSRmat*, ivector*, dCSRmat*, AMG_param*);
37static void interp_STD(dCSRmat*, ivector*, dCSRmat*, iCSRmat*, AMG_param*);
38static void interp_EXT(dCSRmat*, ivector*, dCSRmat*, iCSRmat*, AMG_param*);
39static void amg_interp_trunc(dCSRmat*, AMG_param*);
40
41/*---------------------------------*/
42/*-- Public Functions --*/
43/*---------------------------------*/
44
64void fasp_amg_interp(
65 dCSRmat* A, ivector* vertices, dCSRmat* P, iCSRmat* S, AMG_param* param)
66{
67 const INT coarsening_type = param->coarsening_type;
68 INT interp_type = param->interpolation_type;
69
70 // make sure standard interpolation is used for aggressive coarsening
71 if (coarsening_type == COARSE_AC) interp_type = INTERP_STD;
72
73#if DEBUG_MODE > 0
74 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
75#endif
76
77 switch (interp_type) {
78
79 case INTERP_DIR: // Direct interpolation
80 interp_DIR(A, vertices, P, param);
81 break;
82
83 case INTERP_STD: // Standard interpolation
84 interp_STD(A, vertices, P, S, param);
85 break;
86
87 case INTERP_EXT: // Extended interpolation
88 interp_EXT(A, vertices, P, S, param);
89 break;
90
91 case INTERP_ENG: // Energy-min interpolation defined in PreAMGInterpEM.c
92 fasp_amg_interp_em(A, vertices, P, param);
93 break;
94
95 case INTERP_RDC: // Reduction-based interpolation
96 interp_RDC(A, vertices, P, param);
97 break;
98
99 default:
100 fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
101 }
102
103#if DEBUG_MODE > 0
104 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
105#endif
106}
107
108/*---------------------------------*/
109/*-- Private Functions --*/
110/*---------------------------------*/
111
127static void amg_interp_trunc(dCSRmat* P, AMG_param* param)
128{
129 const INT row = P->row;
130 const INT nnzold = P->nnz;
131 const INT prtlvl = param->print_level;
132 const REAL eps_tr = param->truncation_threshold;
133
134 // local variables
135 INT num_nonzero = 0; // number of non zeros after truncation
136 REAL Min_neg, Max_pos; // min negative and max positive entries
137 REAL Fac_neg, Fac_pos; // factors for negative and positive entries
138 REAL Sum_neg, TSum_neg; // sum and truncated sum of negative entries
139 REAL Sum_pos, TSum_pos; // sum and truncated sum of positive entries
140
141 INT index1 = 0, index2 = 0, begin, end;
142 INT i, j;
143
144#if DEBUG_MODE > 0
145 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
146#endif
147
148 for (i = 0; i < row; ++i) {
149
150 begin = P->IA[i];
151 end = P->IA[i + 1];
152
153 P->IA[i] = num_nonzero;
154 Min_neg = Max_pos = 0;
155 Sum_neg = Sum_pos = 0;
156 TSum_neg = TSum_pos = 0;
157
158 // 1. Summations of positive and negative entries
159 for (j = begin; j < end; ++j) {
160
161 if (P->val[j] > 0) {
162 Sum_pos += P->val[j];
163 Max_pos = MAX(Max_pos, P->val[j]);
164 }
165
166 else {
167 Sum_neg += P->val[j];
168 Min_neg = MIN(Min_neg, P->val[j]);
169 }
170 }
171
172 // Truncate according to max and min values!!!
173 Max_pos *= eps_tr;
174 Min_neg *= eps_tr;
175
176 // 2. Set JA of truncated P
177 for (j = begin; j < end; ++j) {
178
179 if (P->val[j] >= Max_pos) {
180 num_nonzero++;
181 P->JA[index1++] = P->JA[j];
182 TSum_pos += P->val[j];
183 }
184
185 else if (P->val[j] <= Min_neg) {
186 num_nonzero++;
187 P->JA[index1++] = P->JA[j];
188 TSum_neg += P->val[j];
189 }
190 }
191
192 // 3. Compute factors and set values of truncated P
193 if (TSum_pos > SMALLREAL) {
194 Fac_pos = Sum_pos / TSum_pos; // factor for positive entries
195 } else {
196 Fac_pos = 1.0;
197 }
198
199 if (TSum_neg < -SMALLREAL) {
200 Fac_neg = Sum_neg / TSum_neg; // factor for negative entries
201 } else {
202 Fac_neg = 1.0;
203 }
204
205 for (j = begin; j < end; ++j) {
206
207 if (P->val[j] >= Max_pos)
208 P->val[index2++] = P->val[j] * Fac_pos;
209
210 else if (P->val[j] <= Min_neg)
211 P->val[index2++] = P->val[j] * Fac_neg;
212 }
213 }
214
215 // resize the truncated prolongation P
216 P->nnz = P->IA[row] = num_nonzero;
217 P->JA = (INT*)fasp_mem_realloc(P->JA, num_nonzero * sizeof(INT));
218 P->val = (REAL*)fasp_mem_realloc(P->val, num_nonzero * sizeof(REAL));
219
220 if (prtlvl >= PRINT_MOST) {
221 printf("NNZ in prolongator: before truncation = %10d, after = %10d\n", nnzold,
222 num_nonzero);
223 }
224
225#if DEBUG_MODE > 0
226 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
227#endif
228}
229
240static void interp_RDC(dCSRmat* A, ivector* vertices, dCSRmat* P, AMG_param* param)
241{
242 INT row = A->row;
243 INT* vec = vertices->val;
244
245 // local variables
246 INT i, j, k, l, index, idiag;
247 const REAL alpha = 2.0 - 1.0 / param->theta;
248
249 // first pass: P->val
250 for (index = i = 0; i < row; ++i) {
251 if (vec[i] == CGPT) { // identity for coarse points
252 P->val[index++] = 1.0;
253 } else { // interpolation for fine points
254 for (j = A->IA[i]; j < A->IA[i + 1]; ++j) {
255 if (A->JA[j] == i) {
256 idiag = j; // found diagonal entry
257 break;
258 }
259 }
260 REAL Dii = alpha * A->val[idiag];
261 // fill in entries for P
262 for (j = A->IA[i]; j < A->IA[i + 1]; ++j) {
263 if (vec[A->JA[j]] == CGPT) {
264 P->val[index++] = -A->val[j] / Dii;
265 }
266 }
267 }
268 }
269
270 // second pass: P->JA renumber column index for coarse space
271 INT* cindex = (INT*)fasp_mem_calloc(row, sizeof(INT));
272 // index of coarse points from fine-space to coarse-space
273 for (index = i = 0; i < row; ++i) {
274 if (vec[i] == CGPT) cindex[i] = index++;
275 }
276 for (i = 0; i < P->nnz; ++i) P->JA[i] = cindex[P->JA[i]];
277
278 // clean up
279 fasp_mem_free(cindex);
280}
281
302static void interp_DIR(dCSRmat* A, ivector* vertices, dCSRmat* P, AMG_param* param)
303{
304 INT row = A->row;
305 INT* vec = vertices->val;
306
307 // local variables
308 SHORT IS_STRONG; // is the variable strong coupled to i?
309 INT num_pcouple; // number of positive strong couplings
310 INT begin_row, end_row;
311 INT i, j, k, l, index = 0, idiag;
312
313 // a_minus and a_plus for Neighbors and Prolongation support
314 REAL amN, amP, apN, apP;
315 REAL alpha, beta, aii = 0.0;
316
317 // indices of C-nodes
318 INT* cindex = (INT*)fasp_mem_calloc(row, sizeof(INT));
319
320 SHORT use_openmp = FALSE;
321
322#ifdef _OPENMP
323 INT myid, mybegin, myend, stride_i, nthreads;
324 // row = MIN(P->IA[P->row], row);
325 if (MIN(P->IA[P->row], row) > OPENMP_HOLDS) {
326 use_openmp = TRUE;
327 nthreads = fasp_get_num_threads();
328 }
329#endif
330
331 // Step 1. Fill in values for interpolation operator P
332 if (use_openmp) {
333#ifdef _OPENMP
334 stride_i = row / nthreads;
335#pragma omp parallel private(myid, mybegin, myend, i, begin_row, end_row, idiag, aii, \
336 amN, amP, apN, apP, num_pcouple, j, k, alpha, beta, l) \
337 num_threads(nthreads)
338 {
339 myid = omp_get_thread_num();
340 mybegin = myid * stride_i;
341 if (myid < nthreads - 1)
342 myend = mybegin + stride_i;
343 else
344 myend = row;
345 aii = 0.0;
346
347 for (i = mybegin; i < myend; ++i) {
348 begin_row = A->IA[i];
349 end_row = A->IA[i + 1] - 1;
350 for (idiag = begin_row; idiag <= end_row; idiag++) {
351 if (A->JA[idiag] == i) {
352 aii = A->val[idiag];
353 break;
354 }
355 }
356 if (vec[i] == 0) { // if node i is on fine grid
357 amN = 0, amP = 0, apN = 0, apP = 0, num_pcouple = 0;
358 for (j = begin_row; j <= end_row; ++j) {
359 if (j == idiag) continue;
360 for (k = P->IA[i]; k < P->IA[i + 1]; ++k) {
361 if (P->JA[k] == A->JA[j]) break;
362 }
363 if (A->val[j] > 0) {
364 apN += A->val[j];
365 if (k < P->IA[i + 1]) {
366 apP += A->val[j];
367 num_pcouple++;
368 }
369 } else {
370 amN += A->val[j];
371 if (k < P->IA[i + 1]) {
372 amP += A->val[j];
373 }
374 }
375 } // j
376
377 // avoid division by zero for amP and apP
378 amP = (amP < -SMALLREAL) ? amP : -SMALLREAL;
379 apP = (apP > SMALLREAL) ? apP : SMALLREAL;
380 alpha = amN / amP;
381 if (num_pcouple > 0) {
382 beta = apN / apP;
383 } else {
384 beta = 0;
385 aii += apN;
386 }
387 for (j = P->IA[i]; j < P->IA[i + 1]; ++j) {
388 k = P->JA[j];
389 for (l = A->IA[i]; l < A->IA[i + 1]; l++) {
390 if (A->JA[l] == k) break;
391 }
392 if (A->val[l] > 0) {
393 P->val[j] = -beta * A->val[l] / aii;
394 } else {
395 P->val[j] = -alpha * A->val[l] / aii;
396 }
397 }
398 } else if (vec[i] == 2) // if node i is a special fine node
399 {
400
401 } else { // if node i is on coarse grid
402 P->val[P->IA[i]] = 1;
403 }
404 }
405 }
406#endif
407 }
408
409 else {
410
411 for (i = 0; i < row; ++i) {
412
413 begin_row = A->IA[i];
414 end_row = A->IA[i + 1];
415
416 // find diagonal entry first!!!
417 for (idiag = begin_row; idiag < end_row; idiag++) {
418 if (A->JA[idiag] == i) {
419 aii = A->val[idiag];
420 break;
421 }
422 }
423
424 if (vec[i] == FGPT) { // fine grid nodes
425
426 amN = amP = apN = apP = 0.0;
427
428 num_pcouple = 0;
429
430 for (j = begin_row; j < end_row; ++j) {
431
432 if (j == idiag) continue; // skip diagonal
433
434 // check a point strong-coupled to i or not
435 IS_STRONG = FALSE;
436 for (k = P->IA[i]; k < P->IA[i + 1]; ++k) {
437 if (P->JA[k] == A->JA[j]) {
438 IS_STRONG = TRUE;
439 break;
440 }
441 }
442
443 if (A->val[j] > 0) {
444 apN += A->val[j]; // sum up positive entries
445 if (IS_STRONG) {
446 apP += A->val[j];
447 num_pcouple++;
448 }
449 } else {
450 amN += A->val[j]; // sum up negative entries
451 if (IS_STRONG) {
452 amP += A->val[j];
453 }
454 }
455 } // end for j
456
457 // set weight factors
458 // avoid division by zero for amP and apP
459 amP = (amP < -SMALLREAL) ? amP : -SMALLREAL;
460 apP = (apP > SMALLREAL) ? apP : SMALLREAL;
461 alpha = amN / amP;
462 if (num_pcouple > 0) {
463 beta = apN / apP;
464 } else {
465 beta = 0.0;
466 aii += apN;
467 }
468
469 // keep aii inside the loop to avoid floating pt error! --Chensong
470 for (j = P->IA[i]; j < P->IA[i + 1]; ++j) {
471 k = P->JA[j];
472 for (l = A->IA[i]; l < A->IA[i + 1]; l++) {
473 if (A->JA[l] == k) break;
474 }
475 if (A->val[l] > 0) {
476 P->val[j] = -beta * A->val[l] / aii;
477 } else {
478 P->val[j] = -alpha * A->val[l] / aii;
479 }
480 }
481
482 } // end if vec
483
484 else if (vec[i] == CGPT) { // coarse grid nodes
485 P->val[P->IA[i]] = 1.0;
486 }
487 }
488 }
489
490 // Step 2. Generate coarse level indices and set values of P.JA
491 for (index = i = 0; i < row; ++i) {
492 if (vec[i] == CGPT) cindex[i] = index++;
493 }
494 P->col = index;
495
496 if (use_openmp) {
497#ifdef _OPENMP
498 stride_i = P->IA[P->row] / nthreads;
499#pragma omp parallel private(myid, mybegin, myend, i, j) num_threads(nthreads)
500 {
501 myid = omp_get_thread_num();
502 mybegin = myid * stride_i;
503 if (myid < nthreads - 1)
504 myend = mybegin + stride_i;
505 else
506 myend = P->IA[P->row];
507 for (i = mybegin; i < myend; ++i) {
508 j = P->JA[i];
509 P->JA[i] = cindex[j];
510 }
511 }
512#endif
513 } else {
514 for (i = 0; i < P->nnz; ++i) {
515 j = P->JA[i];
516 P->JA[i] = cindex[j];
517 }
518 }
519
520 // clean up
521 fasp_mem_free(cindex);
522 cindex = NULL;
523
524 // Step 3. Truncate the prolongation operator to reduce cost
525 amg_interp_trunc(P, param);
526}
527
547static void interp_STD(dCSRmat* A, ivector* vertices, dCSRmat* P, iCSRmat* S, AMG_param* param)
548{
549 const INT row = A->row;
550 INT* vec = vertices->val;
551
552 // local variables
553 INT i, j, k, l, m, index;
554 REAL alpha = 1.0, factor, alN, alP;
555 REAL akk, akl, aik, aki;
556
557 // indices for coarse neighbor node for every node
558 INT* cindex = (INT*)fasp_mem_calloc(row, sizeof(INT));
559
560 // indices from column number to index in nonzeros in i-th row
561 INT* rindi = (INT*)fasp_mem_calloc(2 * row, sizeof(INT));
562
563 // indices from column number to index in nonzeros in k-th row
564 INT* rindk = (INT*)fasp_mem_calloc(2 * row, sizeof(INT));
565
566 // sums of strongly connected C neighbors
567 REAL* csum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
568
569#if RS_C1
570 // sums of all neighbors except ISPT
571 REAL* psum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
572#endif
573
574 // sums of all neighbors
575 REAL* nsum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
576
577 // diagonal entries
578 REAL* diag = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
579
580 // coefficients hat a_ij for relevant CGPT of the i-th node
581 REAL* Ahat = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
582
583 // Step 0. Prepare diagonal, Cs-sum, and N-sum
584 fasp_iarray_set(row, cindex, -1);
585 fasp_darray_set(row, csum, 0.0);
586 fasp_darray_set(row, nsum, 0.0);
587
588 for (i = 0; i < row; i++) {
589
590 // set flags for strong-connected C nodes
591 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
592 k = S->JA[j];
593 if (vec[k] == CGPT) cindex[k] = i;
594 }
595
596 for (j = A->IA[i]; j < A->IA[i + 1]; j++) {
597 k = A->JA[j];
598
599 if (cindex[k] == i) csum[i] += A->val[j]; // strong C-couplings
600
601 if (k == i) diag[i] = A->val[j];
602#if RS_C1
603 else {
604 nsum[i] += A->val[j];
605 if (vec[k] != ISPT) {
606 psum[i] += A->val[j];
607 }
608 }
609#else
610 else
611 nsum[i] += A->val[j];
612#endif
613 }
614 }
615
616 // Step 1. Fill in values for interpolation operator P
617 for (i = 0; i < row; i++) {
618
619 if (vec[i] == FGPT) {
620#if RS_C1
621 alN = psum[i];
622#else
623 alN = nsum[i];
624#endif
625 alP = csum[i];
626
627 // form the reverse indices for i-th row
628 for (j = A->IA[i]; j < A->IA[i + 1]; j++) rindi[A->JA[j]] = j;
629
630 // clean up Ahat for relevant nodes only
631 for (j = P->IA[i]; j < P->IA[i + 1]; j++) Ahat[P->JA[j]] = 0.0;
632
633 // set values of Ahat
634 Ahat[i] = diag[i];
635
636 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
637
638 k = S->JA[j];
639 aik = A->val[rindi[k]];
640
641 if (vec[k] == CGPT)
642 Ahat[k] += aik;
643
644 else if (vec[k] == FGPT) {
645
646 akk = diag[k];
647
648 // form the reverse indices for k-th row
649 for (m = A->IA[k]; m < A->IA[k + 1]; m++) rindk[A->JA[m]] = m;
650
651 factor = aik / akk;
652
653 // visit the strong-connected C neighbors of k, compute
654 // Ahat in the i-th row, set aki if found
655 aki = 0.0;
656#if 0 // modified by Xiaoqiang Yue 12/25/2013
657 for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
658 l = S->JA[m];
659 akl = A->val[rindk[l]];
660 if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
661 else if ( l == i ) {
662 aki = akl; Ahat[l] -= factor * aki;
663 }
664 } // end for m
665#else
666 for (m = A->IA[k]; m < A->IA[k + 1]; m++) {
667 if (A->JA[m] == i) {
668 aki = A->val[m];
669 Ahat[i] -= factor * aki;
670 }
671 } // end for m
672#endif
673 for (m = S->IA[k]; m < S->IA[k + 1]; m++) {
674 l = S->JA[m];
675 akl = A->val[rindk[l]];
676 if (vec[l] == CGPT) Ahat[l] -= factor * akl;
677 } // end for m
678
679 // compute Cs-sum and N-sum for Ahat
680 alN -= factor * (nsum[k] - aki + akk);
681 alP -= factor * csum[k];
682
683 } // end if vec[k]
684
685 } // end for j
686
687 // Originally: alpha = alN/alP, do this only if P is not empty!
688 if (P->IA[i] < P->IA[i + 1]) alpha = alN / alP;
689
690 // How about positive entries? --Chensong
691 for (j = P->IA[i]; j < P->IA[i + 1]; j++) {
692 k = P->JA[j];
693 P->val[j] = -alpha * Ahat[k] / Ahat[i];
694 }
695
696 }
697
698 else if (vec[i] == CGPT) {
699 P->val[P->IA[i]] = 1.0;
700 }
701
702 } // end for i
703
704 // Step 2. Generate coarse level indices and set values of P.JA
705 for (index = i = 0; i < row; ++i) {
706 if (vec[i] == CGPT) cindex[i] = index++;
707 }
708 P->col = index;
709
710#ifdef _OPENMP
711#pragma omp parallel for private(i, j) if (P->IA[P->row] > OPENMP_HOLDS)
712#endif
713 for (i = 0; i < P->IA[P->row]; ++i) {
714 j = P->JA[i];
715 P->JA[i] = cindex[j];
716 }
717
718 // clean up
719 fasp_mem_free(cindex);
720 cindex = NULL;
721 fasp_mem_free(rindi);
722 rindi = NULL;
723 fasp_mem_free(rindk);
724 rindk = NULL;
725 fasp_mem_free(nsum);
726 nsum = NULL;
727 fasp_mem_free(csum);
728 csum = NULL;
729 fasp_mem_free(diag);
730 diag = NULL;
731 fasp_mem_free(Ahat);
732 Ahat = NULL;
733
734#if RS_C1
735 fasp_mem_free(psum);
736 psum = NULL;
737#endif
738
739 // Step 3. Truncate the prolongation operator to reduce cost
740 amg_interp_trunc(P, param);
741}
742
760static void interp_EXT(dCSRmat* A, ivector* vertices, dCSRmat* P, iCSRmat* S, AMG_param* param)
761{
762 const INT row = A->row;
763 INT* vec = vertices->val;
764
765 // local variables
766 INT i, j, k, l, m, index;
767 REAL alpha = 1.0, factor, alN, alP;
768 REAL akk, akl, aik, aki;
769
770 // indices for coarse neighbor node for every node
771 INT* cindex = (INT*)fasp_mem_calloc(row, sizeof(INT));
772
773 // indices from column number to index in nonzeros in i-th row
774 INT* rindi = (INT*)fasp_mem_calloc(2 * row, sizeof(INT));
775
776 // indices from column number to index in nonzeros in k-th row
777 INT* rindk = (INT*)fasp_mem_calloc(2 * row, sizeof(INT));
778
779 // sums of strongly connected C neighbors
780 REAL* csum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
781
782#if RS_C1
783 // sums of all neighbors except ISPT
784 REAL* psum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
785#endif
786 // sums of all neighbors
787 REAL* nsum = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
788
789 // diagonal entries
790 REAL* diag = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
791
792 // coefficients hat a_ij for relevant CGPT of the i-th node
793 REAL* Ahat = (REAL*)fasp_mem_calloc(row, sizeof(REAL));
794
795 // Step 0. Prepare diagonal, Cs-sum, and N-sum
796 fasp_iarray_set(row, cindex, -1);
797 fasp_darray_set(row, csum, 0.0);
798 fasp_darray_set(row, nsum, 0.0);
799
800 for (i = 0; i < row; i++) {
801
802 // set flags for strong-connected C nodes
803 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
804 k = S->JA[j];
805 if (vec[k] == CGPT) cindex[k] = i;
806 }
807
808 for (j = A->IA[i]; j < A->IA[i + 1]; j++) {
809 k = A->JA[j];
810
811 if (cindex[k] == i) csum[i] += A->val[j]; // strong C-couplings
812
813 if (k == i) diag[i] = A->val[j];
814#if RS_C1
815 else {
816 nsum[i] += A->val[j];
817 if (vec[k] != ISPT) {
818 psum[i] += A->val[j];
819 }
820 }
821#else
822 else
823 nsum[i] += A->val[j];
824#endif
825 }
826 }
827
828 // Step 1. Fill in values for interpolation operator P
829 for (i = 0; i < row; i++) {
830
831 if (vec[i] == FGPT) {
832#if RS_C1
833 alN = psum[i];
834#else
835 alN = nsum[i];
836#endif
837 alP = csum[i];
838
839 // form the reverse indices for i-th row
840 for (j = A->IA[i]; j < A->IA[i + 1]; j++) rindi[A->JA[j]] = j;
841
842 // clean up Ahat for relevant nodes only
843 for (j = P->IA[i]; j < P->IA[i + 1]; j++) Ahat[P->JA[j]] = 0.0;
844
845 // set values of Ahat
846 Ahat[i] = diag[i];
847
848 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
849
850 k = S->JA[j];
851 aik = A->val[rindi[k]];
852
853 if (vec[k] == CGPT)
854 Ahat[k] += aik;
855
856 else if (vec[k] == FGPT) {
857
858 akk = diag[k];
859
860 // form the reverse indices for k-th row
861 for (m = A->IA[k]; m < A->IA[k + 1]; m++) rindk[A->JA[m]] = m;
862
863 factor = aik / akk;
864
865 // visit the strong-connected C neighbors of k, compute
866 // Ahat in the i-th row, set aki if found
867 aki = 0.0;
868#if 0 // modified by Xiaoqiang Yue 12/25/2013
869 for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
870 l = S->JA[m];
871 akl = A->val[rindk[l]];
872 if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
873 else if ( l == i ) {
874 aki = akl; Ahat[l] -= factor * aki;
875 }
876 } // end for m
877#else
878 for (m = A->IA[k]; m < A->IA[k + 1]; m++) {
879 if (A->JA[m] == i) {
880 aki = A->val[m];
881 Ahat[i] -= factor * aki;
882 }
883 } // end for m
884#endif
885 for (m = S->IA[k]; m < S->IA[k + 1]; m++) {
886 l = S->JA[m];
887 akl = A->val[rindk[l]];
888 if (vec[l] == CGPT) Ahat[l] -= factor * akl;
889 } // end for m
890
891 // compute Cs-sum and N-sum for Ahat
892 alN -= factor * (nsum[k] - aki + akk);
893 alP -= factor * csum[k];
894
895 } // end if vec[k]
896
897 } // end for j
898
899 // Originally: alpha = alN/alP, do this only if P is not empty!
900 if (P->IA[i] < P->IA[i + 1]) alpha = alN / alP;
901
902 // How about positive entries? --Chensong
903 for (j = P->IA[i]; j < P->IA[i + 1]; j++) {
904 k = P->JA[j];
905 P->val[j] = -alpha * Ahat[k] / Ahat[i];
906 }
907
908 }
909
910 else if (vec[i] == CGPT) {
911 P->val[P->IA[i]] = 1.0;
912 }
913
914 } // end for i
915
916 // Step 2. Generate coarse level indices and set values of P.JA
917 for (index = i = 0; i < row; ++i) {
918 if (vec[i] == CGPT) cindex[i] = index++;
919 }
920 P->col = index;
921
922#ifdef _OPENMP
923#pragma omp parallel for private(i, j) if (P->IA[P->row] > OPENMP_HOLDS)
924#endif
925 for (i = 0; i < P->IA[P->row]; ++i) {
926 j = P->JA[i];
927 P->JA[i] = cindex[j];
928 }
929
930 // clean up
931 fasp_mem_free(cindex);
932 cindex = NULL;
933 fasp_mem_free(rindi);
934 rindi = NULL;
935 fasp_mem_free(rindk);
936 rindk = NULL;
937 fasp_mem_free(nsum);
938 nsum = NULL;
939 fasp_mem_free(csum);
940 csum = NULL;
941 fasp_mem_free(diag);
942 diag = NULL;
943 fasp_mem_free(Ahat);
944 Ahat = NULL;
945
946#if RS_C1
947 fasp_mem_free(psum);
948 psum = NULL;
949#endif
950
951 // Step 3. Truncate the prolongation operator to reduce cost
952 amg_interp_trunc(P, param);
953}
954
955/*---------------------------------*/
956/*-- End of File --*/
957/*---------------------------------*/
fasp_darray_set
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:41
fasp_iarray_set
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:98
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_realloc
void * fasp_mem_realloc(void *oldmem, const LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: AuxMemory.c:113
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_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_amg_interp_em
void fasp_amg_interp_em(dCSRmat *A, ivector *vertices, dCSRmat *P, AMG_param *param)
Energy-min interpolation.
Definition: PreAMGInterpEM.c:63
fasp_amg_interp
void fasp_amg_interp(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Generate interpolation operator P.
Definition: PreAMGInterp.c:64
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
ISPT
#define ISPT
Definition: fasp_const.h:235
ERROR_AMG_INTERP_TYPE
#define ERROR_AMG_INTERP_TYPE
Definition: fasp_const.h:35
PRINT_MOST
#define PRINT_MOST
Definition: fasp_const.h:77
FGPT
#define FGPT
Definition: fasp_const.h:233
INTERP_EXT
#define INTERP_EXT
Definition: fasp_const.h:225
CGPT
#define CGPT
Definition: fasp_const.h:234
INTERP_RDC
#define INTERP_RDC
Definition: fasp_const.h:226
COARSE_AC
#define COARSE_AC
Definition: fasp_const.h:216
INTERP_STD
#define INTERP_STD
Definition: fasp_const.h:223
INTERP_DIR
#define INTERP_DIR
Definition of interpolation types.
Definition: fasp_const.h:222
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
TRUE
#define TRUE
Definition of logic type.
Definition: fasp_const.h:61
FALSE
#define FALSE
Definition: fasp_const.h:62
SMALLREAL
#define SMALLREAL
Definition: fasp_const.h:256
INTERP_ENG
#define INTERP_ENG
Definition: fasp_const.h:224
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
AMG_param::interpolation_type
SHORT interpolation_type
interpolation type
Definition: fasp.h:524
AMG_param::print_level
SHORT print_level
print level for AMG
Definition: fasp.h:461
AMG_param::coarsening_type
SHORT coarsening_type
coarsening type
Definition: fasp.h:515
AMG_param::theta
REAL theta
theta for reduction-based amg
Definition: fasp.h:593
AMG_param::truncation_threshold
REAL truncation_threshold
truncation threshold
Definition: fasp.h:533
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
iCSRmat
Sparse matrix of INT type in CSR format.
Definition: fasp.h:190
iCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:202
iCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:205
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374