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