Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
interpolation.c
Go to the documentation of this file.
1 
11 #include <math.h>
12 #include <time.h>
13 
14 #ifdef _OPENMP
15 #include <omp.h>
16 #endif
17 
18 #include "fasp.h"
19 #include "fasp_functs.h"
20 
21 static void interp_DIR (dCSRmat *, ivector *, dCSRmat *, AMG_param *);
22 static void interp_DIR1(dCSRmat *, ivector *, dCSRmat *, AMG_param *, INT *);
23 static void interp_STD (dCSRmat *, ivector *, dCSRmat *, iCSRmat *, AMG_param *);
24 
25 /*---------------------------------*/
26 /*-- Public Functions --*/
27 /*---------------------------------*/
28 
49  ivector *vertices,
50  dCSRmat *P,
51  iCSRmat *S,
52  AMG_param *param)
53 {
54  const INT coarsening_type = param->coarsening_type;
55  INT interp_type = param->interpolation_type;
56 
57 #if DEBUG_MODE > 0
58  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
59 #endif
60 
61  // make sure standard interpolation is used for aggressive coarsening
62  if ( coarsening_type == COARSE_AC ) interp_type = INTERP_STD;
63 
64  switch ( interp_type ) {
65 
66  case INTERP_DIR: // Direct interpolation
67  interp_DIR(A, vertices, P, param); break;
68 
69  case INTERP_STD: // Standard interpolation
70  interp_STD(A, vertices, P, S, param); break;
71 
72  case INTERP_ENG: // Energy-min interpolation
73  fasp_amg_interp_em(A, vertices, P, param); break;
74 
75  default:
76  fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
77 
78  }
79 
80 #if DEBUG_MODE > 0
81  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
82 #endif
83 }
84 
106  ivector *vertices,
107  dCSRmat *P,
108  AMG_param *param,
109  iCSRmat *S,
110  INT *icor_ysk)
111 {
112  const INT coarsening_type = param->coarsening_type;
113  INT interp_type = param->interpolation_type;
114 
115 #if DEBUG_MODE > 0
116  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
117 #endif
118 
119  // make sure standard interpolation is used for aggressive coarsening
120  if ( coarsening_type == COARSE_AC ) interp_type = INTERP_STD;
121 
122  switch ( interp_type ) {
123 
124  case INTERP_DIR: // Direct interpolation
125  interp_DIR1(A, vertices, P, param, icor_ysk); break;
126  // interp_DIR(A, vertices, P, param); break;
127 
128  case INTERP_STD: // Standard interpolation
129  interp_STD(A, vertices, P, S, param); break;
130 
131  case INTERP_ENG: // Energy-min interpolation
132  fasp_amg_interp_em(A, vertices, P, param); break;
133 
134  default:
135  fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
136 
137  }
138 
139 #if DEBUG_MODE > 0
140  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
141 #endif
142 }
143 
160  AMG_param *param)
161 {
162  const INT row = P->row;
163  const INT nnzold = P->nnz;
164  const INT prtlvl = param->print_level;
165  const REAL eps_tr = param->truncation_threshold;
166 
167  // local variables
168  INT num_nonzero = 0; // number of non zeros after truncation
169  REAL Min_neg, Max_pos; // min negative and max positive entries
170  REAL Fac_neg, Fac_pos; // factors for negative and positive entries
171  REAL Sum_neg, TSum_neg; // sum and truncated sum of negative entries
172  REAL Sum_pos, TSum_pos; // sum and truncated sum of positive entries
173 
174  INT index1 = 0, index2 = 0, begin, end;
175  INT i, j;
176 
177 #if DEBUG_MODE > 0
178  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
179 #endif
180 
181  for ( i = 0; i < row; ++i ) {
182 
183  begin = P->IA[i]; end = P->IA[i+1];
184 
185  P->IA[i] = num_nonzero;
186  Min_neg = Max_pos = 0;
187  Sum_neg = Sum_pos = 0;
188  TSum_neg = TSum_pos = 0;
189 
190  // 1. Summations of positive and negative entries
191  for ( j = begin; j < end; ++j ) {
192 
193  if ( P->val[j] > 0 ) {
194  Sum_pos += P->val[j];
195  Max_pos = MAX(Max_pos, P->val[j]);
196  }
197 
198  else if ( P->val[j] < 0 ) {
199  Sum_neg += P->val[j];
200  Min_neg = MIN(Min_neg, P->val[j]);
201  }
202 
203  }
204 
205  Max_pos *= eps_tr; Min_neg *= eps_tr;
206 
207  // 2. Set JA of truncated P
208  for ( j = begin; j < end; ++j ) {
209 
210  if ( P->val[j] >= Max_pos ) {
211  num_nonzero++;
212  P->JA[index1++] = P->JA[j];
213  TSum_pos += P->val[j];
214  }
215 
216  else if ( P->val[j] <= Min_neg ) {
217  num_nonzero++;
218  P->JA[index1++] = P->JA[j];
219  TSum_neg += P->val[j];
220  }
221 
222  }
223 
224  // 3. Compute factors and set values of truncated P
225  if ( TSum_pos > SMALLREAL ) {
226  Fac_pos = Sum_pos / TSum_pos; // factor for positive entries
227  }
228  else {
229  Fac_pos = 1.0;
230  }
231 
232  if ( TSum_neg < -SMALLREAL ) {
233  Fac_neg = Sum_neg / TSum_neg; // factor for negative entries
234  }
235  else {
236  Fac_neg = 1.0;
237  }
238 
239  for ( j = begin; j < end; ++j ) {
240 
241  if ( P->val[j] >= Max_pos )
242  P->val[index2++] = P->val[j] * Fac_pos;
243 
244  else if ( P->val[j] <= Min_neg )
245  P->val[index2++] = P->val[j] * Fac_neg;
246  }
247 
248  }
249 
250  // resize the truncated prolongation P
251  P->nnz = P->IA[row] = num_nonzero;
252  P->JA = (INT *)fasp_mem_realloc(P->JA, num_nonzero*sizeof(INT));
253  P->val = (REAL *)fasp_mem_realloc(P->val, num_nonzero*sizeof(REAL));
254 
255  if ( prtlvl >= PRINT_MOST ) {
256  printf("Truncate prolongation, nnz before: %10d, after: %10d\n",
257  nnzold, num_nonzero);
258  }
259 
260 #if DEBUG_MODE > 0
261  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
262 #endif
263 
264 }
265 
266 /*---------------------------------*/
267 /*-- Private Functions --*/
268 /*---------------------------------*/
269 
288 static void interp_DIR (dCSRmat *A,
289  ivector *vertices,
290  dCSRmat *P,
291  AMG_param *param )
292 {
293  INT row = A->row;
294  INT *vec = vertices->val;
295 
296  // local variables
297  SHORT IS_STRONG; // is the variable strong coupled to i?
298  INT num_pcouple; // number of positive strong couplings
299  INT begin_row, end_row;
300  INT i, j, k, l, index = 0, idiag;
301 
302  // a_minus and a_plus for Neighbors and Prolongation support
303  REAL amN, amP, apN, apP;
304  REAL alpha, beta, aii = 0;
305 
306  // indices of C-nodes
307  INT * cindex = (INT *)fasp_mem_calloc(row, sizeof(INT));
308 
309  INT use_openmp = FALSE;
310 
311 #ifdef _OPENMP
312  INT myid, mybegin, myend, stride_i, nthreads;
313  row = MIN(P->IA[P->row], row);
314  if ( row > OPENMP_HOLDS ) {
315  use_openmp = TRUE;
316  nthreads = FASP_GET_NUM_THREADS();
317  }
318 #endif
319 
320  // Step 1. Fill in values for interpolation operator P
321  if (use_openmp) {
322 #ifdef _OPENMP
323  stride_i = row/nthreads;
324 #pragma omp parallel private(myid,mybegin,myend,i,begin_row,end_row,idiag,aii, \
325  amN,amP,apN,apP,num_pcouple,j,k,alpha,beta,l) \
326  num_threads(nthreads)
327  {
328  myid = omp_get_thread_num();
329  mybegin = myid*stride_i;
330  if(myid < nthreads-1) myend = mybegin+stride_i;
331  else myend = row;
332  for (i=mybegin; i<myend; ++i){
333  begin_row=A->IA[i]; end_row=A->IA[i+1]-1;
334  for(idiag=begin_row;idiag<=end_row;idiag++){
335  if (A->JA[idiag]==i) {
336  aii=A->val[idiag];
337  break;
338  }
339  }
340  if(vec[i]==0){ // if node i is on fine grid
341  amN=0, amP=0, apN=0, apP=0, num_pcouple=0;
342  for(j=begin_row;j<=end_row;++j){
343  if(j==idiag) continue;
344  for(k=P->IA[i];k<P->IA[i+1];++k) {
345  if(P->JA[k]==A->JA[j]) break;
346  }
347  if(A->val[j]>0) {
348  apN+=A->val[j];
349  if(k<P->IA[i+1]) {
350  apP+=A->val[j];
351  num_pcouple++;
352  }
353  }
354  else {
355  amN+=A->val[j];
356  if(k<P->IA[i+1]) {
357  amP+=A->val[j];
358  }
359  }
360  } // j
361 
362  alpha=amN/amP;
363  if(num_pcouple>0) {
364  beta=apN/apP;
365  }
366  else {
367  beta=0;
368  aii+=apN;
369  }
370  for(j=P->IA[i];j<P->IA[i+1];++j){
371  k=P->JA[j];
372  for(l=A->IA[i];l<A->IA[i+1];l++){
373  if(A->JA[l]==k) break;
374  }
375  if(A->val[l]>0){
376  P->val[j]=-beta*A->val[l]/aii;
377  }
378  else {
379  P->val[j]=-alpha*A->val[l]/aii;
380  }
381  }
382  }
383  else if(vec[i]==2) // if node i is a special fine node
384  {
385 
386  }
387  else {// if node i is on coarse grid
388  P->val[P->IA[i]]=1;
389  }
390  }
391  }
392 #endif
393  }
394 
395  else {
396 
397  for ( i = 0; i < row; ++i ) {
398 
399  begin_row = A->IA[i]; end_row = A->IA[i+1];
400 
401  // find diagonal entry first!!!
402  for ( idiag = begin_row; idiag < end_row; idiag++ ) {
403  if ( A->JA[idiag] == i ) {
404  aii = A->val[idiag]; break;
405  }
406  }
407 
408  if ( vec[i] == FGPT ) { // fine grid nodes
409 
410  amN = amP = apN = apP = 0.0;
411 
412  num_pcouple = 0;
413 
414  for ( j = begin_row; j < end_row; ++j ) {
415 
416  if ( j == idiag ) continue; // skip diagonal
417 
418  // check a point strong-coupled to i or not
419  IS_STRONG = FALSE;
420  for ( k = P->IA[i]; k < P->IA[i+1]; ++k ) {
421  if ( P->JA[k] == A->JA[j] ) { IS_STRONG = TRUE; break; }
422  }
423 
424  if ( A->val[j] > 0 ) {
425  apN += A->val[j]; // sum up positive entries
426  if ( IS_STRONG ) { apP += A->val[j]; num_pcouple++; }
427  }
428  else {
429  amN += A->val[j]; // sum up negative entries
430  if ( IS_STRONG ) { amP += A->val[j]; }
431  }
432  } // end for j
433 
434  // set weight factors
435  alpha = amN / amP;
436  if ( num_pcouple > 0 ) {
437  beta = apN / apP;
438  }
439  else {
440  beta = 0.0; aii += apN;
441  }
442 
443  // keep aii inside the loop to avoid floating pt error! --Chensong
444  for ( j = P->IA[i]; j < P->IA[i+1]; ++j ) {
445  k = P->JA[j];
446  for ( l = A->IA[i]; l < A->IA[i+1]; l++ ) {
447  if ( A->JA[l] == k ) break;
448  }
449  if ( A->val[l] > 0 ) {
450  P->val[j] = -beta * A->val[l] / aii;
451  }
452  else {
453  P->val[j] = -alpha * A->val[l] / aii;
454  }
455  }
456 
457  } // end if vec
458 
459  else if ( vec[i] == CGPT ) { // coarse grid nodes
460  P->val[P->IA[i]] = 1.0;
461  }
462  }
463  }
464 
465  // Step 2. Generate coarse level indices and set values of P.JA
466  for ( index = i = 0; i < row; ++i ) {
467  if ( vec[i] == CGPT ) cindex[i] = index++;
468  }
469  P->col = index;
470 
471  if (use_openmp) {
472 #ifdef _OPENMP
473  stride_i = P->IA[P->row]/nthreads;
474 #pragma omp parallel private(myid,mybegin,myend,i,j) num_threads(nthreads)
475  {
476  myid = omp_get_thread_num();
477  mybegin = myid*stride_i;
478  if ( myid < nthreads-1 ) myend = mybegin+stride_i;
479  else myend = P->IA[P->row];
480  for ( i = mybegin; i < myend; ++i ) {
481  j = P->JA[i];
482  P->JA[i] = cindex[j];
483  }
484  }
485 #endif
486  }
487  else {
488  for ( i = 0; i < P->nnz; ++i ) {
489  j = P->JA[i];
490  P->JA[i] = cindex[j];
491  }
492  }
493 
494  // clean up
495  fasp_mem_free(cindex);
496 
497  // Step 3. Truncate the prolongation operator to reduce cost
498  fasp_amg_interp_trunc(P, param);
499 }
500 
520 static void interp_STD (dCSRmat *A,
521  ivector *vertices,
522  dCSRmat *P,
523  iCSRmat *S,
524  AMG_param *param)
525 {
526  const INT row = A->row;
527  INT *vec = vertices->val;
528 
529  // local variables
530  INT i, j, k, l, m, index;
531  REAL alpha, factor, alN, alP;
532  REAL akk, akl, aik, aki;
533 
534  // indices for coarse neighbor node for every node
535  INT * cindex = (INT *)fasp_mem_calloc(row, sizeof(INT));
536 
537  // indices from column number to index in nonzeros in i-th row
538  INT * rindi = (INT *)fasp_mem_calloc(2*row, sizeof(INT));
539 
540  // indices from column number to index in nonzeros in k-th row
541  INT * rindk = (INT *)fasp_mem_calloc(2*row, sizeof(INT));
542 
543  // sums of strongly connected C neighbors
544  REAL * csum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
545 
546 #if RS_C1
547  // sums of all neighbors except ISPT
548  REAL * psum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
549 #endif
550  // sums of all neighbors
551  REAL * nsum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
552 
553  // diagonal entries
554  REAL * diag = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
555 
556  // coefficients hat a_ij for relevant CGPT of the i-th node
557  REAL * Ahat = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
558 
559  // Step 0. Prepare diagonal, Cs-sum, and N-sum
560  fasp_iarray_set(row, cindex, -1);
561  fasp_array_set(row, csum, 0.0);
562  fasp_array_set(row, nsum, 0.0);
563 
564  for ( i = 0; i < row; i++ ) {
565 
566  // set flags for strong-connected C nodes
567  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
568  k = S->JA[j];
569  if ( vec[k] == CGPT ) cindex[k] = i;
570  }
571 
572  for ( j = A->IA[i]; j < A->IA[i+1]; j++ ) {
573  k = A->JA[j];
574 
575  if ( cindex[k] == i ) csum[i] += A->val[j]; // strong C-couplings
576 
577  if ( k == i ) diag[i] = A->val[j];
578 #if RS_C1
579  else {
580  nsum[i] += A->val[j];
581  if ( vec[k] != ISPT ) {
582  psum[i] += A->val[j];
583  }
584  }
585 #else
586  else nsum[i] += A->val[j];
587 #endif
588  }
589 
590  }
591 
592  // Step 1. Fill in values for interpolation operator P
593  for ( i = 0; i < row; i++ ) {
594 
595  if ( vec[i] == FGPT ) {
596 #if RS_C1
597  alN = psum[i];
598 #else
599  alN = nsum[i];
600 #endif
601  alP = csum[i];
602 
603  // form the reverse indices for i-th row
604  for ( j = A->IA[i]; j < A->IA[i+1]; j++ ) rindi[A->JA[j]] = j;
605 
606  // clean up Ahat for relevant nodes only
607  for ( j = P->IA[i]; j < P->IA[i+1]; j++ ) Ahat[P->JA[j]] = 0.0;
608 
609  // set values of Ahat
610  Ahat[i] = diag[i];
611 
612  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
613 
614  k = S->JA[j]; aik = A->val[rindi[k]];
615 
616  if ( vec[k] == CGPT ) Ahat[k] += aik;
617 
618  else if ( vec[k] == FGPT ) {
619 
620  akk = diag[k];
621 
622  // form the reverse indices for k-th row
623  for ( m = A->IA[k]; m < A->IA[k+1]; m++ ) rindk[A->JA[m]] = m;
624 
625  factor = aik / akk;
626 
627  // visit the strong-connected C neighbors of k, compute
628  // Ahat in the i-th row, set aki if found
629  aki = 0.0;
630 #if 0 // modified by Xiaoqiang Yue 12/25/2013
631  for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
632  l = S->JA[m];
633  akl = A->val[rindk[l]];
634  if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
635  else if ( l == i ) {
636  aki = akl; Ahat[l] -= factor * aki;
637  }
638  } // end for m
639 #else
640  for ( m = A->IA[k]; m < A->IA[k+1]; m++ ) {
641  if ( A->JA[m] == i ) {
642  aki = A->val[m];
643  Ahat[i] -= factor * aki;
644  }
645  } // end for m
646 #endif
647  for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
648  l = S->JA[m];
649  akl = A->val[rindk[l]];
650  if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
651  } // end for m
652 
653  // compute Cs-sum and N-sum for Ahat
654  alN -= factor * (nsum[k]-aki+akk);
655  alP -= factor * csum[k];
656 
657  } // end if vec[k]
658 
659  } // end for j
660 
661  // How about positive entries? --Chensong
662  alpha = alN/alP;
663  for ( j = P->IA[i]; j < P->IA[i+1]; j++ ) {
664  k = P->JA[j];
665  P->val[j] = -alpha*Ahat[k]/Ahat[i];
666  }
667 
668  }
669 
670  else if ( vec[i] == CGPT ) {
671  P->val[P->IA[i]] = 1.0;
672  }
673 
674  } // end for i
675 
676  // Step 2. Generate coarse level indices and set values of P.JA
677  for ( index = i = 0; i < row; ++i ) {
678  if ( vec[i] == CGPT ) cindex[i] = index++;
679  }
680  P->col = index;
681 
682 #ifdef _OPENMP
683 #pragma omp parallel for private(i,j) if(P->IA[P->row]>OPENMP_HOLDS)
684 #endif
685  for ( i = 0; i < P->IA[P->row]; ++i ) {
686  j = P->JA[i];
687  P->JA[i] = cindex[j];
688  }
689 
690  // clean up
691  fasp_mem_free(cindex);
692  fasp_mem_free(rindi);
693  fasp_mem_free(rindk);
694  fasp_mem_free(nsum);
695 
696 #if RS_C1
697  fasp_mem_free(psum);
698 #endif
699 
700  fasp_mem_free(csum);
701  fasp_mem_free(diag);
702  fasp_mem_free(Ahat);
703 
704  // Step 3. Truncate the prolongation operator to reduce cost
705  fasp_amg_interp_trunc(P, param);
706 }
707 
725 static void interp_EXT (dCSRmat *A,
726  ivector *vertices,
727  dCSRmat *P,
728  iCSRmat *S,
729  AMG_param *param)
730 {
731  const INT row = A->row;
732  INT *vec = vertices->val;
733 
734  // local variables
735  INT i, j, k, l, m, index;
736  REAL alpha, factor, alN, alP;
737  REAL akk, akl, aik, aki;
738 
739  // indices for coarse neighbor node for every node
740  INT * cindex = (INT *)fasp_mem_calloc(row, sizeof(INT));
741 
742  // indices from column number to index in nonzeros in i-th row
743  INT * rindi = (INT *)fasp_mem_calloc(2*row, sizeof(INT));
744 
745  // indices from column number to index in nonzeros in k-th row
746  INT * rindk = (INT *)fasp_mem_calloc(2*row, sizeof(INT));
747 
748  // sums of strongly connected C neighbors
749  REAL * csum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
750 
751 #if RS_C1
752  // sums of all neighbors except ISPT
753  REAL * psum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
754 #endif
755  // sums of all neighbors
756  REAL * nsum = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
757 
758  // diagonal entries
759  REAL * diag = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
760 
761  // coefficients hat a_ij for relevant CGPT of the i-th node
762  REAL * Ahat = (REAL *)fasp_mem_calloc(row, sizeof(REAL));
763 
764  // Step 0. Prepare diagonal, Cs-sum, and N-sum
765  fasp_iarray_set(row, cindex, -1);
766  fasp_array_set(row, csum, 0.0);
767  fasp_array_set(row, nsum, 0.0);
768 
769  for ( i = 0; i < row; i++ ) {
770 
771  // set flags for strong-connected C nodes
772  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
773  k = S->JA[j];
774  if ( vec[k] == CGPT ) cindex[k] = i;
775  }
776 
777  for ( j = A->IA[i]; j < A->IA[i+1]; j++ ) {
778  k = A->JA[j];
779 
780  if ( cindex[k] == i ) csum[i] += A->val[j]; // strong C-couplings
781 
782  if ( k == i ) diag[i] = A->val[j];
783 #if RS_C1
784  else {
785  nsum[i] += A->val[j];
786  if ( vec[k] != ISPT ) {
787  psum[i] += A->val[j];
788  }
789  }
790 #else
791  else nsum[i] += A->val[j];
792 #endif
793  }
794 
795  }
796 
797  // Step 1. Fill in values for interpolation operator P
798  for ( i = 0; i < row; i++ ) {
799 
800  if ( vec[i] == FGPT ) {
801 #if RS_C1
802  alN = psum[i];
803 #else
804  alN = nsum[i];
805 #endif
806  alP = csum[i];
807 
808  // form the reverse indices for i-th row
809  for ( j = A->IA[i]; j < A->IA[i+1]; j++ ) rindi[A->JA[j]] = j;
810 
811  // clean up Ahat for relevant nodes only
812  for ( j = P->IA[i]; j < P->IA[i+1]; j++ ) Ahat[P->JA[j]] = 0.0;
813 
814  // set values of Ahat
815  Ahat[i] = diag[i];
816 
817  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
818 
819  k = S->JA[j]; aik = A->val[rindi[k]];
820 
821  if ( vec[k] == CGPT ) Ahat[k] += aik;
822 
823  else if ( vec[k] == FGPT ) {
824 
825  akk = diag[k];
826 
827  // form the reverse indices for k-th row
828  for ( m = A->IA[k]; m < A->IA[k+1]; m++ ) rindk[A->JA[m]] = m;
829 
830  factor = aik / akk;
831 
832  // visit the strong-connected C neighbors of k, compute
833  // Ahat in the i-th row, set aki if found
834  aki = 0.0;
835 #if 0 // modified by Xiaoqiang Yue 12/25/2013
836  for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
837  l = S->JA[m];
838  akl = A->val[rindk[l]];
839  if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
840  else if ( l == i ) {
841  aki = akl; Ahat[l] -= factor * aki;
842  }
843  } // end for m
844 #else
845  for ( m = A->IA[k]; m < A->IA[k+1]; m++ ) {
846  if ( A->JA[m] == i ) {
847  aki = A->val[m];
848  Ahat[i] -= factor * aki;
849  }
850  } // end for m
851 #endif
852  for ( m = S->IA[k]; m < S->IA[k+1]; m++ ) {
853  l = S->JA[m];
854  akl = A->val[rindk[l]];
855  if ( vec[l] == CGPT ) Ahat[l] -= factor * akl;
856  } // end for m
857 
858  // compute Cs-sum and N-sum for Ahat
859  alN -= factor * (nsum[k]-aki+akk);
860  alP -= factor * csum[k];
861 
862  } // end if vec[k]
863 
864  } // end for j
865 
866  // How about positive entries? --Chensong
867  alpha = alN/alP;
868  for ( j = P->IA[i]; j < P->IA[i+1]; j++ ) {
869  k = P->JA[j];
870  P->val[j] = -alpha*Ahat[k]/Ahat[i];
871  }
872 
873  }
874 
875  else if ( vec[i] == CGPT ) {
876  P->val[P->IA[i]] = 1.0;
877  }
878 
879  } // end for i
880 
881  // Step 2. Generate coarse level indices and set values of P.JA
882  for ( index = i = 0; i < row; ++i ) {
883  if ( vec[i] == CGPT ) cindex[i] = index++;
884  }
885  P->col = index;
886 
887 #ifdef _OPENMP
888 #pragma omp parallel for private(i,j) if(P->IA[P->row]>OPENMP_HOLDS)
889 #endif
890  for ( i = 0; i < P->IA[P->row]; ++i ) {
891  j = P->JA[i];
892  P->JA[i] = cindex[j];
893  }
894 
895  // clean up
896  fasp_mem_free(cindex);
897  fasp_mem_free(rindi);
898  fasp_mem_free(rindk);
899  fasp_mem_free(nsum);
900 
901 #if RS_C1
902  fasp_mem_free(psum);
903 #endif
904 
905  fasp_mem_free(csum);
906  fasp_mem_free(diag);
907  fasp_mem_free(Ahat);
908 
909  // Step 3. Truncate the prolongation operator to reduce cost
910  fasp_amg_interp_trunc(P, param);
911 }
924 static void get_nwidth (dCSRmat *A,
925  INT *nbl_ptr,
926  INT *nbr_ptr )
927 {
928 #ifdef _OPENMP
929  INT *IA = A->IA;
930  INT *JA = A->JA;
931  INT myid, mybegin, myend, nthreads;
932  INT max_l, max_r;
933  INT i, end_row_A, j;
934 
935  nthreads = FASP_GET_NUM_THREADS();
936 
937  INT *max_left_right = (INT *)fasp_mem_calloc(2*nthreads, sizeof(INT));
938 
939 #pragma omp parallel for private(myid,mybegin,myend,max_l,max_r,i,end_row_A,j)
940  for (myid = 0; myid < nthreads; myid ++)
941  {
942  FASP_GET_START_END(myid, nthreads, A->row, &mybegin, &myend);
943  max_l = 0;
944  max_r = 0;
945  for (i = mybegin; i < myend; i ++) {
946  end_row_A = IA[i+1];
947  for (j = IA[i]; j < end_row_A; j ++) {
948  max_l = MAX(i-JA[j], max_l);
949  max_r = MAX(JA[j]-i, max_r);
950  }
951  }
952  max_left_right[myid*2] = max_l;
953  max_left_right[myid*2+1] = max_r;
954  }
955  max_l = max_left_right[0];
956  max_r = max_left_right[1];
957  for (i = 1; i < nthreads; i ++) {
958  max_l = MAX(max_l, max_left_right[i*2]);
959  max_r = MAX(max_r, max_left_right[i*2+1]);
960  }
961  fasp_mem_free(max_left_right);
962  *nbl_ptr = max_l;
963  *nbr_ptr = max_r;
964 #endif
965 }
966 
978 static void mod_cindex (INT nrows,
979  INT *cindex)
980 {
981 #ifdef _OPENMP
982  INT myid, mybegin, myend, i, nthreads;
983 
984  nthreads = FASP_GET_NUM_THREADS();
985 
986 #pragma omp parallel for private(myid,mybegin,myend,i)
987  for (myid = 0; myid < nthreads; myid ++)
988  {
989  FASP_GET_START_END(myid, nthreads, nrows, &mybegin, &myend);
990  if (myid == 0)
991  {
992  mybegin ++;
993  }
994  for (i = mybegin; i < myend; i ++)
995  {
996  if (cindex[i] < cindex[i-1])
997  {
998  cindex[i] = cindex[i-1];
999  }
1000  }
1001  }
1002 #endif
1003 }
1004 
1023 static void get_cindex (INT nrows,
1024  INT ncols,
1025  INT *cindex,
1026  INT nbl_ysk,
1027  INT nbr_ysk,
1028  INT *CF_marker,
1029  INT *icor_ysk)
1030 {
1031 #ifdef _OPENMP
1032  INT myid, FiveMyid, mybegin, myend, min_A, max_A, i, first_f_node, min_P, max_P, myend_minus_one;
1033  INT lengthAA = 0, lengthPP = 0;
1034  INT nthreads;
1035 
1036  nthreads = FASP_GET_NUM_THREADS();
1037 
1038 #pragma omp parallel for private(myid,FiveMyid,mybegin,myend,min_A,max_A,i, \
1039  first_f_node,min_P,max_P,myend_minus_one) \
1040  reduction(+: lengthAA,lengthPP)
1041  for (myid = 0; myid < nthreads; myid ++)
1042  {
1043  FiveMyid = myid * 5;
1044  FASP_GET_START_END(myid, nthreads, ncols, &mybegin, &myend);
1045  icor_ysk[FiveMyid] = mybegin;
1046  if (mybegin == myend) {
1047  lengthAA = 0;
1048  lengthPP = 0;
1049  icor_ysk[FiveMyid+1] = 0;
1050  icor_ysk[FiveMyid+3] = 0;
1051  }
1052  else {
1053  first_f_node = fasp_BinarySearch(cindex, mybegin, nrows);
1054  for (i = first_f_node; i > -1; i --) {
1055  if (cindex[i] != mybegin) {
1056  break;
1057  }
1058  }
1059  min_A = i + 1;
1060  min_A = MAX(0, min_A-2*nbl_ysk);
1061  myend_minus_one = myend - 1;
1062  first_f_node = fasp_BinarySearch(cindex, myend_minus_one, nrows);
1063  for (i = first_f_node; i > -1; i --) {
1064  if (cindex[i] != myend_minus_one) {
1065  max_A = i;
1066  break;
1067  }
1068  }
1069  max_A = MIN(nrows, max_A+2*nbr_ysk+1);
1070  lengthAA = max_A - min_A + 2;
1071  icor_ysk[FiveMyid+1] = lengthAA;
1072  icor_ysk[FiveMyid+2] = min_A;
1073  for (i = min_A; i >= 0; i --) {
1074  if (CF_marker[i] == 0) {
1075  first_f_node = i;
1076  break;
1077  }
1078  }
1079  if (i == -1) {
1080  min_P = 0;
1081  }
1082  else {
1083  first_f_node -= nbl_ysk;
1084  if (first_f_node <= 0) {
1085  min_P = 0;
1086  }
1087  else {
1088  for (i = first_f_node; i >= 0; i --) {
1089  if (CF_marker[i] == 1) {
1090  min_P = cindex[i];
1091  break;
1092  }
1093  }
1094  if (i == -1) {
1095  min_P = 0;
1096  }
1097  }
1098  }
1099  for (i = max_A-1; i < nrows; i ++) {
1100  if (CF_marker[i] == 0) {
1101  first_f_node = i;
1102  break;
1103  }
1104  }
1105  if (i == nrows) {
1106  max_P = ncols;
1107  }
1108  else {
1109  first_f_node += nbr_ysk;
1110  if (first_f_node >= nrows) {
1111  max_P = ncols;
1112  }
1113  else {
1114  for (i = first_f_node; i < nrows; i ++) {
1115  if (CF_marker[i] == 1) {
1116  max_P = cindex[i] + 1;
1117  break;
1118  }
1119  }
1120  if (i == nrows) {
1121  max_P = ncols;
1122  }
1123  }
1124  }
1125  lengthPP = max_P - min_P + 2;
1126  icor_ysk[FiveMyid+3] = lengthPP;
1127  icor_ysk[FiveMyid+4] = min_P;
1128  }
1129  }
1130  icor_ysk[5*nthreads] = lengthAA;
1131  icor_ysk[5*nthreads+1] = lengthPP;
1132 #endif
1133 }
1134 
1153 static void interp_DIR1 (dCSRmat *A,
1154  ivector *vertices,
1155  dCSRmat *Ptr,
1156  AMG_param *param,
1157  INT *icor_ysk)
1158 {
1159  REAL eps_tr = param->truncation_threshold;
1160  REAL amN, amP, apN, apP;
1161  REAL alpha, beta, aii=0;
1162  INT *vec = vertices->val;
1163  INT num_pcouple, idiag;
1164 
1165  INT i,j,k,l,index=0;
1166  INT begin_row, end_row;
1167  INT myid;
1168  INT mybegin;
1169  INT myend;
1170  INT *indexs = NULL;
1171  INT shift;
1172  INT nbl_ysk, nbr_ysk;
1173 
1174  /* Generate interpolation P */
1175  dCSRmat P=fasp_dcsr_create(Ptr->row,Ptr->col,Ptr->nnz);
1176 
1177  INT nthreads = 1, use_openmp = FALSE;
1178 
1179 #ifdef _OPENMP
1180  INT row = MIN(P.IA[P.row], A->row);
1181  if ( row > OPENMP_HOLDS ) {
1182  use_openmp = TRUE;
1183  nthreads = FASP_GET_NUM_THREADS();
1184  }
1185 #endif
1186 
1187  /* step 1: Find first the structure IA of P */
1188  fasp_iarray_cp(P.row+1, Ptr->IA, P.IA);
1189 
1190  /* step 2: Find the structure JA of P */
1191  fasp_iarray_cp(P.nnz, Ptr->JA, P.JA);
1192 
1193  /* step 3: Fill the data of P */
1194  if (use_openmp) {
1195 #ifdef _OPENMP
1196 #pragma omp parallel for private(myid,mybegin,myend,i,begin_row,end_row,idiag,aii, \
1197  amN,amP,apN,apP,num_pcouple,j,k,alpha,beta,l)
1198 #endif
1199  for (myid = 0; myid < nthreads; myid++ )
1200  {
1201  FASP_GET_START_END(myid, nthreads, A->row, &mybegin, &myend);
1202  for (i=mybegin; i<myend; ++i)
1203  {
1204  begin_row=A->IA[i]; end_row=A->IA[i+1]-1;
1205 
1206  for(idiag=begin_row;idiag<=end_row;idiag++) {
1207  if (A->JA[idiag]==i) {
1208  aii=A->val[idiag];
1209  break;
1210  }
1211  }
1212 
1213  if(vec[i]==0) // if node i is on fine grid
1214  {
1215  amN=0, amP=0, apN=0, apP=0, num_pcouple=0;
1216 
1217  for(j=begin_row;j<=end_row;++j)
1218  {
1219  if(j==idiag) continue;
1220 
1221  for(k=Ptr->IA[i];k<Ptr->IA[i+1];++k) {
1222  if(Ptr->JA[k]==A->JA[j]) break;
1223  }
1224 
1225  if(A->val[j]>0) {
1226  apN+=A->val[j];
1227  if(k<Ptr->IA[i+1]) {
1228  apP+=A->val[j];
1229  num_pcouple++;
1230  }
1231  }
1232  else
1233  {
1234  amN+=A->val[j];
1235  if(k<Ptr->IA[i+1]) {
1236  amP+=A->val[j];
1237  }
1238  }
1239  } // j
1240 
1241  alpha=amN/amP;
1242  if(num_pcouple>0) {
1243  beta=apN/apP;
1244  }
1245  else {
1246  beta=0;
1247  aii+=apN;
1248  }
1249 
1250  for(j=P.IA[i];j<P.IA[i+1];++j)
1251  {
1252  k=P.JA[j];
1253  for(l=A->IA[i];l<A->IA[i+1];l++)
1254  {
1255  if(A->JA[l]==k) break;
1256  }
1257 
1258  if(A->val[l]>0)
1259  {
1260  P.val[j]=-beta*A->val[l]/aii;
1261  }
1262  else
1263  {
1264  P.val[j]=-alpha*A->val[l]/aii;
1265  }
1266  }
1267  }
1268  else if(vec[i]==2)
1269  {
1270  // if node i is a special fine node
1271  }
1272  else
1273  {
1274  // if node i is on coarse grid
1275  P.val[P.IA[i]]=1;
1276  }
1277  }
1278  }
1279  }
1280  else {
1281  for(i=0;i<A->row;++i) {
1282  begin_row=A->IA[i]; end_row=A->IA[i+1]-1;
1283 
1284  for(idiag=begin_row;idiag<=end_row;idiag++) {
1285  if (A->JA[idiag]==i) {
1286  aii=A->val[idiag];
1287  break;
1288  }
1289  }
1290 
1291  if(vec[i]==0) // if node i is on fine grid
1292  {
1293  amN=0, amP=0, apN=0, apP=0, num_pcouple=0;
1294 
1295  for(j=begin_row;j<=end_row;++j)
1296  {
1297  if(j==idiag) continue;
1298 
1299  for(k=Ptr->IA[i];k<Ptr->IA[i+1];++k) {
1300  if(Ptr->JA[k]==A->JA[j]) break;
1301  }
1302 
1303  if(A->val[j]>0) {
1304  apN+=A->val[j];
1305  if(k<Ptr->IA[i+1]) {
1306  apP+=A->val[j];
1307  num_pcouple++;
1308  }
1309  }
1310  else
1311  {
1312  amN+=A->val[j];
1313  if(k<Ptr->IA[i+1]) {
1314  amP+=A->val[j];
1315  }
1316  }
1317  } // j
1318 
1319  alpha=amN/amP;
1320  if(num_pcouple>0) {
1321  beta=apN/apP;
1322  }
1323  else {
1324  beta=0;
1325  aii+=apN;
1326  }
1327 
1328  for(j=P.IA[i];j<P.IA[i+1];++j)
1329  {
1330  k=P.JA[j];
1331  for(l=A->IA[i];l<A->IA[i+1];l++)
1332  {
1333  if(A->JA[l]==k) break;
1334  }
1335 
1336  if(A->val[l]>0)
1337  {
1338  P.val[j]=-beta*A->val[l]/aii;
1339  }
1340  else
1341  {
1342  P.val[j]=-alpha*A->val[l]/aii;
1343  }
1344  }
1345  }
1346  else if(vec[i]==2) // if node i is a special fine node
1347  {
1348  }
1349  else // if node i is on coarse grid
1350  {
1351  P.val[P.IA[i]]=1;
1352  }
1353  }
1354  }
1355  fasp_mem_free(Ptr->IA);
1356  fasp_mem_free(Ptr->JA);
1357  fasp_mem_free(Ptr->val);
1358 
1359  INT *cindex;
1360 
1361  cindex=(INT*)fasp_mem_calloc(A->row, sizeof(INT));
1362 
1363 #if CHMEM_MODE
1364  total_alloc_mem += (A->row)*sizeof(INT);
1365 #endif
1366 
1367  // The following is one of OPTIMAL parts ...0802...
1368  // Generate cindex in parallel
1369  if (use_openmp) {
1370 #ifdef _OPENMP
1371 #pragma omp master
1372  {
1373  indexs = (INT *)fasp_mem_calloc(nthreads, sizeof(INT));
1374  }
1375 #pragma omp parallel for private(myid, mybegin, myend, index, i)
1376 #endif
1377  for (myid = 0; myid < nthreads; myid ++)
1378  {
1379  FASP_GET_START_END(myid, nthreads, A->row, &mybegin, &myend);
1380  index = 0;
1381  for (i=mybegin;i<myend;++i) {
1382  if(vec[i]==1)
1383  {
1384  cindex[i]=index;
1385  index++;
1386  }
1387  }
1388  indexs[myid] = index;
1389  }
1390  for (i = 1; i < nthreads; i ++) {
1391  indexs[i] += indexs[i-1];
1392  }
1393 #ifdef _OPENMP
1394 #pragma omp parallel for private(myid, mybegin, myend, shift, i)
1395 #endif
1396  for (myid = 0; myid < nthreads; myid ++)
1397  {
1398  FASP_GET_START_END(myid, nthreads, A->row, &mybegin, &myend);
1399  shift = 0;
1400  if (myid > 0) {
1401  shift = indexs[myid-1];
1402  }
1403  for (i=mybegin;i<myend;++i) {
1404  if(vec[i]==1)
1405  {
1406  cindex[i] += shift;
1407  }
1408  }
1409  }
1410  fasp_mem_free(indexs);
1411  }
1412  else {
1413  index=0;
1414  for(i=0;i<A->row;++i) {
1415  if(vec[i]==1)
1416  {
1417  cindex[i]=index;
1418  index++;
1419  }
1420  }
1421  }
1422 
1423  if (use_openmp) {
1424 #ifdef _OPENMP
1425 #pragma omp parallel for private(myid, mybegin,myend,i,j)
1426 #endif
1427  for (myid = 0; myid < nthreads; myid++ )
1428  {
1429  FASP_GET_START_END(myid, nthreads, P.IA[P.row], &mybegin, &myend);
1430  for (i=mybegin; i<myend; ++i)
1431  {
1432  j=P.JA[i];
1433  P.JA[i]=cindex[j];
1434  }
1435  }
1436  }
1437  else {
1438  for(i=0;i<P.IA[P.row];++i)
1439  {
1440  j=P.JA[i];
1441  P.JA[i]=cindex[j];
1442  }
1443  }
1444  // if (Ptr->col > OPENMP_HOLDS) {
1445  // The following is another OPTIMAL part ...0802...
1447  get_nwidth(A, &nbl_ysk, &nbr_ysk);
1448  mod_cindex(A->row, cindex);
1449  get_cindex(A->row, Ptr->col, cindex, nbl_ysk, nbr_ysk, vec, icor_ysk);
1451  // }
1452  fasp_mem_free(cindex);
1453 
1454  /* Truncation of interpolation */
1455  REAL Min_neg, Max_pos;
1456  REAL Sum_neg, Sum_pos;
1457  REAL TSum_neg, TSum_pos;
1458  INT num_neg, num_pos;
1459  INT num_nonzero=0;
1460 
1461 
1462  Ptr->val=(REAL*)fasp_mem_calloc(P.IA[Ptr->row],sizeof(REAL));
1463 
1464 #if CHMEM_MODE
1465  total_alloc_mem += (P.IA[Ptr->row])*sizeof(REAL);
1466 #endif
1467  Ptr->JA=(INT*)fasp_mem_calloc(P.IA[Ptr->row],sizeof(INT));
1468 #if CHMEM_MODE
1469  total_alloc_mem += (P.IA[Ptr->row])*sizeof(INT);
1470 #endif
1471  Ptr->IA=(INT*)fasp_mem_calloc(Ptr->row+1, sizeof(INT));
1472 #if CHMEM_MODE
1473  total_alloc_mem += (Ptr->row+1)*sizeof(INT);
1474 #endif
1475 
1476  INT index1=0, index2=0;
1477  for(i=0;i<P.row;++i)
1478  {
1479  Min_neg=0;
1480  Max_pos=0;
1481  Sum_neg=0;
1482  Sum_pos=0;
1483  TSum_neg=0;
1484  TSum_pos=0;
1485  num_neg=0;
1486  num_pos=0;
1487 
1488  Ptr->IA[i]-=num_nonzero;
1489 
1490  for(j=P.IA[i];j<P.IA[i+1];++j)
1491  {
1492  if(P.val[j]<0)
1493  {
1494  Sum_neg+=P.val[j];
1495  if(P.val[j]<Min_neg)
1496  {
1497  Min_neg=P.val[j];
1498  }
1499  }
1500 
1501  if(P.val[j]>0)
1502  {
1503  Sum_pos+=P.val[j];
1504  if(P.val[j]>Max_pos)
1505  {
1506  Max_pos=P.val[j];
1507  }
1508  }
1509  }
1510 
1511  for(j=P.IA[i];j<P.IA[i+1];++j)
1512  {
1513  if(P.val[j]<0)
1514  {
1515  if(P.val[j]>Min_neg*eps_tr)
1516  {
1517  num_neg++;
1518  }
1519  else
1520  {
1521  num_nonzero--;
1522  }
1523  }
1524 
1525  if(P.val[j]>0)
1526  {
1527  if(P.val[j]<Max_pos*eps_tr)
1528  {
1529  num_pos++;
1530  }
1531  else
1532  {
1533  num_nonzero--;
1534  }
1535  }
1536  }
1537 
1538  // step 2: Find the structure JA and fill the data A of Ptr
1539  for(j=P.IA[i];j<P.IA[i+1];++j)
1540  {
1541  if(P.val[j]<0)
1542  {
1543  if(!(P.val[j]>Min_neg*eps_tr))
1544  {
1545  Ptr->JA[index1]=P.JA[j];
1546  TSum_neg+=P.val[j];
1547  index1++;
1548  }
1549  }
1550 
1551  if(P.val[j]>0)
1552  {
1553  if(!(P.val[j]<Max_pos*eps_tr))
1554  {
1555  Ptr->JA[index1]=P.JA[j];
1556  TSum_pos+=P.val[j];
1557  index1++;
1558  }
1559  }
1560  }
1561 
1562  // step 3: Fill the data A of Ptr
1563  for(j=P.IA[i];j<P.IA[i+1];++j)
1564  {
1565  if(P.val[j]<0)
1566  {
1567  if(!(P.val[j]>Min_neg*eps_tr))
1568  {
1569  Ptr->val[index2]=P.val[j]/TSum_neg*Sum_neg;
1570  index2++;
1571  }
1572  }
1573 
1574  if(P.val[j]>0)
1575  {
1576  if(!(P.val[j]<Max_pos*eps_tr))
1577  {
1578  Ptr->val[index2]=P.val[j]/TSum_pos*Sum_pos;
1579  index2++;
1580  }
1581  }
1582  }
1583  }
1584  Ptr->IA[P.row]-=num_nonzero;
1585  Ptr->nnz=Ptr->IA[Ptr->row];
1586 
1587  Ptr->JA=(INT*)fasp_mem_realloc(Ptr->JA, Ptr->IA[Ptr->row]*sizeof(INT));
1588  Ptr->val=(REAL*)fasp_mem_realloc(Ptr->val, Ptr->IA[Ptr->row]*sizeof(REAL));
1589 
1590  fasp_mem_free(P.IA);
1591  fasp_mem_free(P.JA);
1592  fasp_mem_free(P.val);
1593 }
1594 
1595 /*---------------------------------*/
1596 /*-- End of File --*/
1597 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: sparse_csr.c:34
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
Parameters for AMG solver.
Definition: fasp.h:583
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
SHORT coarsening_type
coarsening type
Definition: fasp.h:643
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
#define COARSE_AC
Definition: fasp_const.h:200
void fasp_amg_interp_em(dCSRmat *A, ivector *vertices, dCSRmat *P, AMG_param *param)
Energy-min interpolation.
#define ERROR_AMG_INTERP_TYPE
Definition: fasp_const.h:42
REAL * val
nonzero entries of A
Definition: fasp.h:166
#define INTERP_DIR
Definition of interpolation types.
Definition: fasp_const.h:206
#define PRINT_MOST
Definition: fasp_const.h:83
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
#define INTERP_STD
Definition: fasp_const.h:207
#define INT
Definition: fasp.h:64
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:193
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_amg_interp1(dCSRmat *A, ivector *vertices, dCSRmat *P, AMG_param *param, iCSRmat *S, INT *icor_ysk)
Generate interpolation operator P.
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
void fasp_iarray_cp(const INT n, INT *x, INT *y)
Copy an array to the other y=x.
Definition: array.c:185
INT fasp_BinarySearch(INT *list, const INT value, const INT nlist)
Binary Search.
Definition: ordering.c:30
#define OPENMP_HOLDS
Definition: fasp_const.h:248
INT nnz
number of nonzero entries
Definition: fasp.h:157
REAL truncation_threshold
truncation threshold
Definition: fasp.h:658
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:190
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
void fasp_amg_interp(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Generate interpolation operator P.
Definition: interpolation.c:48
INT * val
actual vector entries
Definition: fasp.h:362
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define CGPT
Definition: fasp_const.h:216
void fasp_amg_interp_trunc(dCSRmat *P, AMG_param *param)
Truncation step for prolongation operators.
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MIN(a, b)
Definition: fasp.h:73
SHORT interpolation_type
interpolation type
Definition: fasp.h:649
SHORT print_level
print level for AMG
Definition: fasp.h:589
unsigned INT total_alloc_mem
Definition: memory.c:34
#define ISPT
Definition: fasp_const.h:217
#define INTERP_ENG
Definition: fasp_const.h:208
#define FALSE
Definition: fasp_const.h:68
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: array.c:107
Sparse matrix of INT type in CSR format.
Definition: fasp.h:178
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define SMALLREAL
Definition: fasp_const.h:238
#define FGPT
Definition: fasp_const.h:215