# Diff of /trunk/paso/src/AMG_Prolongation.c

revision 3641 by gross, Fri Jan 14 00:04:53 2011 UTC revision 3642 by caltinay, Thu Oct 27 03:41:51 2011 UTC
# Line 29  Line 29
29
30  /**************************************************************  /**************************************************************
31
32      Methods nessecary for AMG preconditioner      Methods necessary for AMG preconditioner
33
34      construct n x n_C the prolongation matrix P from A_p.      construct n x n_C the prolongation matrix P from A_p.
35
36      the columns in A_p to be considered are marked by counter_C[n] where      The columns in A_p to be considered are marked by counter_C[n] where
37      an unknown i is to be considered in P is marked by 0<= counter_C[i] < n_C      an unknown i to be considered in P is marked by 0<= counter_C[i] < n_C
38      and counter_C[i]  gives the new column number in P. S defines the strong connections.      and counter_C[i]  gives the new column number in P.
39        S defines the strong connections.
40
41      The pattern of P is formed as follows:      The pattern of P is formed as follows:
42
43      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
44      If row i is not C, then P[i,j] <> 0 if counter_C[k]==j (k in C) and (i,k) is strong connection.        If row i is not C, then P[i,j] <> 0 if counter_C[k]==j (k in C) and (i,k) is a strong connection.
45
46      two settings for P are implemented (see below)      Two settings for P are implemented (see below)
47
48  */  */
49
# Line 143  Paso_SparseMatrix* Paso_Preconditioner_A Line 144  Paso_SparseMatrix* Paso_Preconditioner_A
144      Direct Prolongation:      Direct Prolongation:
145      -------------------      -------------------
146
147      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
148      If row i is not C, then P[i,j] = - a[i] * A[i,k]/A[i,i] with j=counter_C[k]>=0 and k in S      If row i is not C, then P[i,j] = - a[i] * A[i,k]/A[i,i] with j=counter_C[k]>=0 and k in S
149
150     and    a[i]=     and    a[i]=
# Line 189  void Paso_Preconditioner_AMG_setDirectPr Line 190  void Paso_Preconditioner_AMG_setDirectPr
190             }             }
191
192             if (counter_C[j]>=0) {             if (counter_C[j]>=0) {
193            /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */            /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
194            start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);            start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
195            where_p=(index_t*)bsearch(&(counter_C[j]), start_p,            where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
196                          P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],                          P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
197                          sizeof(index_t),                          sizeof(index_t),
198                          Paso_comparIndex);                          Paso_comparIndex);
199            if (! (where_p == NULL) ) { /* yes i stronly connect with j */            if (! (where_p == NULL) ) { /* yes i strongly connected with j */
200              offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);              offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
201              P_p->val[offset]=A_ij; /* will be modified later */              P_p->val[offset]=A_ij; /* will be modified later */
202              if (A_ij< 0)  {              if (A_ij< 0)  {
# Line 287  void Paso_Preconditioner_AMG_setDirectPr Line 288  void Paso_Preconditioner_AMG_setDirectPr
288            }            }
289
290            if (counter_C[j]>=0) {            if (counter_C[j]>=0) {
291               /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */               /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
292               start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);               start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
293               where_p=(index_t*)bsearch(&(counter_C[j]), start_p,               where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
294                           P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],                           P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
295                           sizeof(index_t),                           sizeof(index_t),
296                           Paso_comparIndex);                           Paso_comparIndex);
297               if (! (where_p == NULL) ) { /* yes i stronly connect with j */               if (! (where_p == NULL) ) { /* yes i strongly connected with j */
298                    offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                    offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
299                    for (ib =0; ib<row_block; ++ib) {                    for (ib =0; ib<row_block; ++ib) {
300                   A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];                   A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
# Line 354  void Paso_Preconditioner_AMG_setDirectPr Line 355  void Paso_Preconditioner_AMG_setDirectPr
355      Classic Prolongation:      Classic Prolongation:
356      -------------------      -------------------
357
358      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise      If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
359      If row i is not C, then P[i,j] = - 1/a[i] * ( A[i,k] + sum_{l} A[i,l]*A+[l,k]/B[i,k])      If row i is not C, then P[i,j] = - 1/a[i] * ( A[i,k] + sum_{l} A[i,l]*A+[l,k]/B[i,k])
360               where the summation over l is considering columns which are strongly connected               where the summation over l is considering columns which are strongly connected
361               to i (l in S[i]) and not in C (counter_C[l]<0) and               to i (l in S[i]) and not in C (counter_C[l]<0) and
# Line 392  void Paso_Preconditioner_AMG_setClassicP Line 393  void Paso_Preconditioner_AMG_setClassicP
393             const index_t *start_s = &(S[offset_S[i]]);             const index_t *start_s = &(S[offset_S[i]]);
394             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
395                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
396                /* this loop sums up the weak connections in a and creates a list of the strong connected columns                /* this loop sums up the weak connections in A and creates a list of the strong connected columns
397                                                                        which are not in C (=no interpolation nodes) */                   which are not in C (=no interpolation nodes) */
398                const double A_ii = A_p->val[ptr_main_A[i]];                const double A_ii = A_p->val[ptr_main_A[i]];
399                double a=A_ii;                double a=A_ii;
400
# Line 403  void Paso_Preconditioner_AMG_setClassicP Line 404  void Paso_Preconditioner_AMG_setClassicP
404                   if ( (i!=j) && (degree_S[j]>0) ) {                   if ( (i!=j) && (degree_S[j]>0) ) {
405                      /* is (i,j) a strong connection ?*/                      /* is (i,j) a strong connection ?*/
406                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
407                  if (where_s == NULL) { /* weak connections are accummulated */                  if (where_s == NULL) { /* weak connections are accumulated */
408                          a+=A_ij;                            a+=A_ij;
409                      } else {   /* yes i stronly connect with j */                      } else {   /* yes i strongly connected with j */
410                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */
411                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);
412                                 if (where_p == NULL)  {                                 if (where_p == NULL)  {
413                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setBoomerProlongation: interpolation point is missing.");                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");
414                                 } else {                                 } else {
415                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
416                                  P_p->val[offset]+=A_ij;                                  P_p->val[offset]+=A_ij;
# Line 456  void Paso_Preconditioner_AMG_setClassicP Line 457  void Paso_Preconditioner_AMG_setClassicP
457                     }                     }
458                }                }
459            }            }
460          }  /* endo of row i loop */          }  /* end of row i loop */
461          TMPMEMFREE(D_s);          TMPMEMFREE(D_s);
462          TMPMEMFREE(D_s_offset);          TMPMEMFREE(D_s_offset);
463       }    /* end of parallel region */       }    /* end of parallel region */
# Line 491  void Paso_Preconditioner_AMG_setClassicP Line 492  void Paso_Preconditioner_AMG_setClassicP
492             const index_t *start_s = &(S[offset_S[i]]);             const index_t *start_s = &(S[offset_S[i]]);
493             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
494                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
495                /* this loop sums up the weak connections in a and creates a list of the strong connected columns                /* this loop sums up the weak connections in A and creates a list of the strong connected columns
496                                                                        which are not in C (=no interpolation nodes) */                   which are not in C (=no interpolation nodes) */
497                const double *A_ii = &(A_p->val[ptr_main_A[i]*A_block]);                const double *A_ii = &(A_p->val[ptr_main_A[i]*A_block]);
498                for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];                for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];
499
# Line 504  void Paso_Preconditioner_AMG_setClassicP Line 505  void Paso_Preconditioner_AMG_setClassicP
505                   if ( (i!=j) && (degree_S[j]>0) ) {                   if ( (i!=j) && (degree_S[j]>0) ) {
506                      /* is (i,j) a strong connection ?*/                      /* is (i,j) a strong connection ?*/
507                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
508                  if (where_s == NULL) { /* weak connections are accummulated */                  if (where_s == NULL) { /* weak connections are accumulated */
509                          for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];                          for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
510                      } else {   /* yes i stronly connect with j */                      } else {   /* yes i strongly connected with j */
511                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */
512                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);
513                                 if (where_p == NULL)  {                                 if (where_p == NULL)  {
514                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setBoomerProlongation: interpolation point is missing.");                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");
515                                 } else {                                 } else {
516                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
517                                      for (ib=0; ib<row_block; ib++) P_p->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];                                      for (ib=0; ib<row_block; ib++) P_p->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
# Line 565  void Paso_Preconditioner_AMG_setClassicP Line 566  void Paso_Preconditioner_AMG_setClassicP
566                     }                     }
567                }                }
568            }            }
569          }  /* endo of row i loop */          }  /* end of row i loop */
570          TMPMEMFREE(D_s);          TMPMEMFREE(D_s);
571          TMPMEMFREE(D_s_offset);          TMPMEMFREE(D_s_offset);
572          TMPMEMFREE(a);          TMPMEMFREE(a);

Legend:
 Removed from v.3641 changed lines Added in v.3642