29#include "fasp_functs.h"
76 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
79 switch ( interp_type ) {
82 interp_DIR(A, vertices, P, param);
break;
85 interp_STD(A, vertices, P, S, param);
break;
88 interp_EXT(A, vertices, P, S, param);
break;
99 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
122static void amg_interp_trunc (
dCSRmat *P,
126 const INT nnzold = P->
nnz;
132 REAL Min_neg, Max_pos;
133 REAL Fac_neg, Fac_pos;
134 REAL Sum_neg, TSum_neg;
135 REAL Sum_pos, TSum_pos;
137 INT index1 = 0, index2 = 0, begin, end;
141 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
144 for ( i = 0; i < row; ++i ) {
146 begin = P->
IA[i]; end = P->
IA[i+1];
148 P->
IA[i] = num_nonzero;
149 Min_neg = Max_pos = 0;
150 Sum_neg = Sum_pos = 0;
151 TSum_neg = TSum_pos = 0;
154 for ( j = begin; j < end; ++j ) {
156 if ( P->
val[j] > 0 ) {
157 Sum_pos += P->
val[j];
158 Max_pos =
MAX(Max_pos, P->
val[j]);
162 Sum_neg += P->
val[j];
163 Min_neg =
MIN(Min_neg, P->
val[j]);
169 Max_pos *= eps_tr; Min_neg *= eps_tr;
172 for ( j = begin; j < end; ++j ) {
174 if ( P->
val[j] >= Max_pos ) {
176 P->
JA[index1++] = P->
JA[j];
177 TSum_pos += P->
val[j];
180 else if ( P->
val[j] <= Min_neg ) {
182 P->
JA[index1++] = P->
JA[j];
183 TSum_neg += P->
val[j];
190 Fac_pos = Sum_pos / TSum_pos;
197 Fac_neg = Sum_neg / TSum_neg;
203 for ( j = begin; j < end; ++j ) {
205 if ( P->
val[j] >= Max_pos )
206 P->
val[index2++] = P->
val[j] * Fac_pos;
208 else if ( P->
val[j] <= Min_neg )
209 P->
val[index2++] = P->
val[j] * Fac_neg;
215 P->
nnz = P->
IA[row] = num_nonzero;
220 printf(
"NNZ in prolongator: before truncation = %10d, after = %10d\n",
221 nnzold, num_nonzero);
225 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
250static void interp_DIR (
dCSRmat *A,
261 INT begin_row, end_row;
262 INT i, j, k, l, index = 0, idiag;
265 REAL amN, amP, apN, apP;
266 REAL alpha, beta, aii = 0.0;
274 INT myid, mybegin, myend, stride_i, nthreads;
278 nthreads = fasp_get_num_threads();
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) \
290 myid = omp_get_thread_num();
291 mybegin = myid*stride_i;
292 if (myid < nthreads-1) myend = mybegin+stride_i;
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) {
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;
334 for (j=P->
IA[i];j<P->IA[i+1];++j) {
336 for (l=A->
IA[i];l<A->IA[i+1];l++) {
337 if (A->
JA[l]==k)
break;
340 P->
val[j]=-beta*A->
val[l]/aii;
343 P->
val[j]=-alpha*A->
val[l]/aii;
361 for ( i = 0; i < row; ++i ) {
363 begin_row = A->
IA[i]; end_row = A->
IA[i+1];
366 for ( idiag = begin_row; idiag < end_row; idiag++ ) {
367 if ( A->
JA[idiag] == i ) {
368 aii = A->
val[idiag];
break;
372 if ( vec[i] ==
FGPT ) {
374 amN = amP = apN = apP = 0.0;
378 for ( j = begin_row; j < end_row; ++j ) {
380 if ( j == idiag )
continue;
384 for ( k = P->
IA[i]; k < P->IA[i+1]; ++k ) {
385 if ( P->
JA[k] == A->
JA[j] ) { IS_STRONG =
TRUE;
break; }
388 if ( A->
val[j] > 0 ) {
390 if ( IS_STRONG ) { apP += A->
val[j]; num_pcouple++; }
394 if ( IS_STRONG ) { amP += A->
val[j]; }
400 if ( num_pcouple > 0 ) {
404 beta = 0.0; aii += apN;
408 for ( j = P->
IA[i]; j < P->IA[i+1]; ++j ) {
410 for ( l = A->
IA[i]; l < A->IA[i+1]; l++ ) {
411 if ( A->
JA[l] == k )
break;
413 if ( A->
val[l] > 0 ) {
414 P->
val[j] = -beta * A->
val[l] / aii;
417 P->
val[j] = -alpha * A->
val[l] / aii;
423 else if ( vec[i] ==
CGPT ) {
424 P->
val[P->
IA[i]] = 1.0;
430 for ( index = i = 0; i < row; ++i ) {
431 if ( vec[i] ==
CGPT ) cindex[i] = index++;
437 stride_i = P->
IA[P->
row]/nthreads;
438#pragma omp parallel private(myid,mybegin,myend,i,j) num_threads(nthreads)
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 ) {
446 P->
JA[i] = cindex[j];
452 for ( i = 0; i < P->
nnz; ++i ) {
454 P->
JA[i] = cindex[j];
462 amg_interp_trunc(P, param);
484static void interp_STD (
dCSRmat *A,
494 INT i, j, k, l, m, index;
495 REAL alpha = 1.0, factor, alN, alP;
496 REAL akk, akl, aik, aki;
529 for ( i = 0; i < row; i++ ) {
532 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
534 if ( vec[k] ==
CGPT ) cindex[k] = i;
537 for ( j = A->
IA[i]; j < A->IA[i+1]; j++ ) {
540 if ( cindex[k] == i ) csum[i] += A->
val[j];
542 if ( k == i ) diag[i] = A->
val[j];
545 nsum[i] += A->
val[j];
546 if ( vec[k] !=
ISPT ) {
547 psum[i] += A->
val[j];
551 else nsum[i] += A->
val[j];
558 for ( i = 0; i < row; i++ ) {
560 if ( vec[i] ==
FGPT ) {
569 for ( j = A->
IA[i]; j < A->IA[i+1]; j++ ) rindi[A->
JA[j]] = j;
572 for ( j = P->
IA[i]; j < P->IA[i+1]; j++ ) Ahat[P->
JA[j]] = 0.0;
577 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
579 k = S->
JA[j]; aik = A->
val[rindi[k]];
581 if ( vec[k] ==
CGPT ) Ahat[k] += aik;
583 else if ( vec[k] ==
FGPT ) {
588 for ( m = A->
IA[k]; m < A->IA[k+1]; m++ ) rindk[A->
JA[m]] = m;
596 for ( m = S->
IA[k]; m < S->IA[k+1]; m++ ) {
598 akl = A->
val[rindk[l]];
599 if ( vec[l] ==
CGPT ) Ahat[l] -= factor * akl;
601 aki = akl; Ahat[l] -= factor * aki;
605 for ( m = A->
IA[k]; m < A->IA[k+1]; m++ ) {
606 if ( A->
JA[m] == i ) {
608 Ahat[i] -= factor * aki;
612 for ( m = S->
IA[k]; m < S->IA[k+1]; m++ ) {
614 akl = A->
val[rindk[l]];
615 if ( vec[l] ==
CGPT ) Ahat[l] -= factor * akl;
619 alN -= factor * (nsum[k]-aki+akk);
620 alP -= factor * csum[k];
627 if ( P->
IA[i] < P->
IA[i+1] ) alpha = alN/alP;
630 for ( j = P->
IA[i]; j < P->IA[i+1]; j++ ) {
632 P->
val[j] = -alpha*Ahat[k]/Ahat[i];
637 else if ( vec[i] ==
CGPT ) {
638 P->
val[P->
IA[i]] = 1.0;
644 for ( index = i = 0; i < row; ++i ) {
645 if ( vec[i] ==
CGPT ) cindex[i] = index++;
650#pragma omp parallel for private(i,j) if(P->IA[P->row]>OPENMP_HOLDS)
652 for ( i = 0; i < P->
IA[P->
row]; ++i ) {
654 P->
JA[i] = cindex[j];
671 amg_interp_trunc(P, param);
691static void interp_EXT (
dCSRmat *A,
701 INT i, j, k, l, m, index;
702 REAL alpha = 1.0, factor, alN, alP;
703 REAL akk, akl, aik, aki;
735 for ( i = 0; i < row; i++ ) {
738 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
740 if ( vec[k] ==
CGPT ) cindex[k] = i;
743 for ( j = A->
IA[i]; j < A->IA[i+1]; j++ ) {
746 if ( cindex[k] == i ) csum[i] += A->
val[j];
748 if ( k == i ) diag[i] = A->
val[j];
751 nsum[i] += A->
val[j];
752 if ( vec[k] !=
ISPT ) {
753 psum[i] += A->
val[j];
757 else nsum[i] += A->
val[j];
764 for ( i = 0; i < row; i++ ) {
766 if ( vec[i] ==
FGPT ) {
775 for ( j = A->
IA[i]; j < A->IA[i+1]; j++ ) rindi[A->
JA[j]] = j;
778 for ( j = P->
IA[i]; j < P->IA[i+1]; j++ ) Ahat[P->
JA[j]] = 0.0;
783 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
785 k = S->
JA[j]; aik = A->
val[rindi[k]];
787 if ( vec[k] ==
CGPT ) Ahat[k] += aik;
789 else if ( vec[k] ==
FGPT ) {
794 for ( m = A->
IA[k]; m < A->IA[k+1]; m++ ) rindk[A->
JA[m]] = m;
802 for ( m = S->
IA[k]; m < S->IA[k+1]; m++ ) {
804 akl = A->
val[rindk[l]];
805 if ( vec[l] ==
CGPT ) Ahat[l] -= factor * akl;
807 aki = akl; Ahat[l] -= factor * aki;
811 for ( m = A->
IA[k]; m < A->IA[k+1]; m++ ) {
812 if ( A->
JA[m] == i ) {
814 Ahat[i] -= factor * aki;
818 for ( m = S->
IA[k]; m < S->IA[k+1]; m++ ) {
820 akl = A->
val[rindk[l]];
821 if ( vec[l] ==
CGPT ) Ahat[l] -= factor * akl;
825 alN -= factor * (nsum[k]-aki+akk);
826 alP -= factor * csum[k];
833 if ( P->
IA[i] < P->
IA[i+1] ) alpha = alN/alP;
836 for ( j = P->
IA[i]; j < P->IA[i+1]; j++ ) {
838 P->
val[j] = -alpha*Ahat[k]/Ahat[i];
843 else if ( vec[i] ==
CGPT ) {
844 P->
val[P->
IA[i]] = 1.0;
850 for ( index = i = 0; i < row; ++i ) {
851 if ( vec[i] ==
CGPT ) cindex[i] = index++;
856#pragma omp parallel for private(i,j) if(P->IA[P->row]>OPENMP_HOLDS)
858 for ( i = 0; i < P->
IA[P->
row]; ++i ) {
860 P->
JA[i] = cindex[j];
877 amg_interp_trunc(P, param);
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
void * fasp_mem_realloc(void *oldmem, const LONGLONG tsize)
Reallocate, initiate, and check memory.
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
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.
Main header file for the FASP project.
#define SHORT
FASP integer and floating point numbers.
#define MAX(a, b)
Definition of max, min, abs.
#define ERROR_AMG_INTERP_TYPE
#define INTERP_DIR
Definition of interpolation types.
#define TRUE
Definition of logic type.
Parameters for AMG methods.
SHORT interpolation_type
interpolation type
SHORT print_level
print level for AMG
SHORT coarsening_type
coarsening type
REAL truncation_threshold
truncation threshold
Sparse matrix of REAL type in CSR format.
INT col
column of matrix A, n
REAL * val
nonzero entries of A
INT row
row number of matrix A, m
INT * IA
integer array of row pointers, the size is m+1
INT nnz
number of nonzero entries
INT * JA
integer array of column indexes, the size is nnz
Sparse matrix of INT type in CSR format.
INT * IA
integer array of row pointers, the size is m+1
INT * JA
integer array of column indexes, the size is nnz
Vector with n entries of INT type.
INT * val
actual vector entries