/[escript]/trunk/paso/src/AMG.c
ViewVC logotype

Diff of /trunk/paso/src/AMG.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3439 by plaub, Mon Jan 10 02:06:07 2011 UTC revision 3440 by gross, Fri Jan 14 00:04:53 2011 UTC
# Line 93  Paso_Preconditioner_LocalAMG* Paso_Preco Line 93  Paso_Preconditioner_LocalAMG* Paso_Preco
93    Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;    Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;
94    const dim_t n=A_p->numRows;    const dim_t n=A_p->numRows;
95    const dim_t n_block=A_p->row_block_size;    const dim_t n_block=A_p->row_block_size;
96    index_t* split_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree=NULL;    index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree_S=NULL;
97    dim_t n_F=0, n_C=0, i;    dim_t n_F=0, n_C=0, i;
98    double time0=0;    double time0=0;
99    const double theta = options->coarsening_threshold;    const double theta = options->coarsening_threshold;
# Line 113  Paso_Preconditioner_LocalAMG* Paso_Preco Line 113  Paso_Preconditioner_LocalAMG* Paso_Preco
113                        - 'SIZE' = min_coarse_matrix_size exceeded                        - 'SIZE' = min_coarse_matrix_size exceeded
114                        - 'LEVEL' = level_max exceeded                        - 'LEVEL' = level_max exceeded
115            */            */
116            printf("Paso_Preconditioner: Stopping condition: ");            printf("Paso_Preconditioner AMG: termination of coarsening by ");
117    
118            if (A_p->pattern->len >= options->min_coarse_sparsity * n * n)            if (A_p->pattern->len >= options->min_coarse_sparsity * n * n)
119                printf("SPAR ");                printf("SPAR ");
# Line 135  Paso_Preconditioner_LocalAMG* Paso_Preco Line 135  Paso_Preconditioner_LocalAMG* Paso_Preco
135    }    }
136       /* Start Coarsening : */       /* Start Coarsening : */
137    
138       split_marker=TMPMEMALLOC(n,index_t);       F_marker=TMPMEMALLOC(n,index_t);
139       counter=TMPMEMALLOC(n,index_t);       counter=TMPMEMALLOC(n,index_t);
140       degree=TMPMEMALLOC(n, dim_t);       degree_S=TMPMEMALLOC(n, dim_t);
141       S=TMPMEMALLOC(A_p->pattern->len, index_t);       S=TMPMEMALLOC(A_p->pattern->len, index_t);
142       if ( !( Esys_checkPtr(split_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree) || Esys_checkPtr(S) ) ) {       if ( !( Esys_checkPtr(F_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree_S) || Esys_checkPtr(S) ) ) {
143       /*       /*
144            set splitting of unknows:            set splitting of unknows:
145                
146           */           */
147       time0=Esys_timer();       time0=Esys_timer();
148       if (n_block>1) {       if (n_block>1) {
149             Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree, S, theta,tau);             Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree_S, S, theta,tau);
150       } else {       } else {
151             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree, S, theta,tau);             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree_S, S, theta,tau);
152       }       }
153       Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree, S, split_marker, options->usePanel);       Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);
154    
155             /* in BoomerAMG interpolation is used FF connectiovity is required :*/
156             if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
157                                 Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);  
158       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
159            
160       if (Esys_noError() ) {       if (Esys_noError() ) {
161          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
162          for (i = 0; i < n; ++i) split_marker[i]= (split_marker[i] == PASO_AMG_IN_F);          for (i = 0; i < n; ++i) F_marker[i]=(F_marker[i] ==  PASO_AMG_IN_F);
163            
164          /*          /*
165             count number of unkowns to be eliminated:             count number of unkowns to be eliminated:
166          */          */
167          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);          n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);
168          n_C=n-n_F;          n_C=n-n_F;
169          if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);          if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
170            
# Line 210  Paso_Preconditioner_LocalAMG* Paso_Preco Line 214  Paso_Preconditioner_LocalAMG* Paso_Preco
214                 {                 {
215                    #pragma omp for schedule(static)                    #pragma omp for schedule(static)
216                    for (i = 0; i < n; ++i) {                    for (i = 0; i < n; ++i) {
217                   if  (split_marker[i]) rows_in_F[counter[i]]=i;                   if  (F_marker[i]) rows_in_F[counter[i]]=i;
218                    }                    }
219                 }                 }
220                 /*  create mask of C nodes with value >-1 gives new id */                 /*  create mask of C nodes with value >-1 gives new id */
221                 i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);                 i=Paso_Util_cumsum_maskedFalse(n,counter, F_marker);
222    
223                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
224                 for (i = 0; i < n; ++i) {                 for (i = 0; i < n; ++i) {
225                    if  (split_marker[i]) {                    if  (F_marker[i]) {
226                   mask_C[i]=-1;                   mask_C[i]=-1;
227                    } else {                    } else {
228                   mask_C[i]=counter[i];;                   mask_C[i]=counter[i];;
229                    }                    }
230                 }                 }
231                 /*                 /*
232                    get Restriction :                      get Prolongation :    
233                 */                                   */                  
234                 time0=Esys_timer();                 time0=Esys_timer();
235                 out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);                 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);
236                 if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);                 if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);
237              }              }
238              /*                    /*      
239                 construct Prolongation operator as transposed of restriction operator:                 construct Restriction operator as transposed of Prolongation operator:
240              */              */
241              if ( Esys_noError()) {              if ( Esys_noError()) {
242                 time0=Esys_timer();                 time0=Esys_timer();
# Line 294  Paso_Preconditioner_LocalAMG* Paso_Preco Line 298  Paso_Preconditioner_LocalAMG* Paso_Preco
298       }       }
299    }    }
300    TMPMEMFREE(counter);    TMPMEMFREE(counter);
301    TMPMEMFREE(split_marker);    TMPMEMFREE(F_marker);
302    TMPMEMFREE(degree);    TMPMEMFREE(degree_S);
303    TMPMEMFREE(S);    TMPMEMFREE(S);
304    
305    if (Esys_noError()) {    if (Esys_noError()) {
# Line 370  in the sense that |A_{ij}| >= theta * ma Line 374  in the sense that |A_{ij}| >= theta * ma
374  */  */
375    
376  void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,  void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,
377                        dim_t *degree, index_t *S,                        dim_t *degree_S, index_t *S,
378                        const double theta, const double tau)                        const double theta, const double tau)
379  {  {
380     const dim_t n=A->numRows;     const dim_t n=A->numRows;
# Line 408  void Paso_Preconditioner_AMG_setStrongCo Line 412  void Paso_Preconditioner_AMG_setStrongCo
412            kdeg++;            kdeg++;
413             }             }
414          }          }
415       }           }
416       degree[i]=kdeg;       degree_S[i]=kdeg;
417        }        }
418    
419  }  }
# Line 421  void Paso_Preconditioner_AMG_setStrongCo Line 425  void Paso_Preconditioner_AMG_setStrongCo
425  in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F  in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
426  */  */
427  void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,  void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,
428                              dim_t *degree, index_t *S,                              dim_t *degree_S, index_t *S,
429                              const double theta, const double tau)                              const double theta, const double tau)
430    
431  {  {
# Line 447  void Paso_Preconditioner_AMG_setStrongCo Line 451  void Paso_Preconditioner_AMG_setStrongCo
451          max_offdiagonal = 0.;          max_offdiagonal = 0.;
452          sum_row=0;          sum_row=0;
453          main_row=0;          main_row=0;
454    
455          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
456             j=A->pattern->index[iptr];             j=A->pattern->index[iptr];
457             fnorm=0;             fnorm=0;
# Line 474  void Paso_Preconditioner_AMG_setStrongCo Line 479  void Paso_Preconditioner_AMG_setStrongCo
479            }            }
480             }             }
481          }          }
482          degree[i]=kdeg;          degree_S[i]=kdeg;
483       }             }      
484       TMPMEMFREE(rtmp);       TMPMEMFREE(rtmp);
485        } /* end of parallel region */        } /* end of parallel region */
# Line 482  void Paso_Preconditioner_AMG_setStrongCo Line 487  void Paso_Preconditioner_AMG_setStrongCo
487  }    }  
488    
489  /* the runge stueben coarsening algorithm: */  /* the runge stueben coarsening algorithm: */
490  void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset,  void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset_S,
491                          const dim_t* degree, const index_t* S,                          const dim_t* degree_S, const index_t* S,
492                          index_t*split_marker, const bool_t usePanel)                          index_t*split_marker, const bool_t usePanel)
493  {  {
494        
495     index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;     index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;
496     dim_t i,k, p, q, *degreeT=NULL, len_panel, len_panel_new;     dim_t i,k, p, q, *degree_ST=NULL, len_panel, len_panel_new;
497     register index_t j, itmp;     register index_t j, itmp;
498        
499     if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */     if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
500        
501     lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);     lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);
502     degreeT=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degreeT);     degree_ST=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degree_ST);
503     ST=TMPMEMALLOC(offset[n], index_t);  Esys_checkPtr(ST);     ST=TMPMEMALLOC(offset_S[n], index_t);  Esys_checkPtr(ST);
504     if (usePanel) {     if (usePanel) {
505        notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);        notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);
506        panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);        panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);
# Line 508  void Paso_Preconditioner_AMG_RungeStuebe Line 513  void Paso_Preconditioner_AMG_RungeStuebe
513        /* those unknows which are not influenced go into F, the rest is available for F or C */        /* those unknows which are not influenced go into F, the rest is available for F or C */
514        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
515        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
516       degreeT[i]=0;       degree_ST[i]=0;
517       if (degree[i]>0) {       if (degree_S[i]>0) {
518          lambda[i]=0;          lambda[i]=0;
519          split_marker[i]=PASO_AMG_UNDECIDED;          split_marker[i]=PASO_AMG_UNDECIDED;
520       } else {       } else {
# Line 519  void Paso_Preconditioner_AMG_RungeStuebe Line 524  void Paso_Preconditioner_AMG_RungeStuebe
524        }        }
525        /* create transpose :*/        /* create transpose :*/
526        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
527          for (p=0; p<degree[i]; ++p) {          for (p=0; p<degree_S[i]; ++p) {
528             j=S[offset[i]+p];             j=S[offset_S[i]+p];
529             ST[offset[j]+degreeT[j]]=i;             ST[offset_S[j]+degree_ST[j]]=i;
530             degreeT[j]++;             degree_ST[j]++;
531          }          }
532        }        }
533        /* lambda[i] = |undecided k in ST[i]| + 2 * |F-unknown in ST[i]| */        /* lambda[i] = |undecided k in ST[i]| + 2 * |F-unknown in ST[i]| */
# Line 530  void Paso_Preconditioner_AMG_RungeStuebe Line 535  void Paso_Preconditioner_AMG_RungeStuebe
535        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
536       if (split_marker[i]==PASO_AMG_UNDECIDED) {       if (split_marker[i]==PASO_AMG_UNDECIDED) {
537          itmp=lambda[i];          itmp=lambda[i];
538          for (p=0; p<degreeT[i]; ++p) {          for (p=0; p<degree_ST[i]; ++p) {
539             j=ST[offset[i]+p];             j=ST[offset_S[i]+p];
540             if (split_marker[j]==PASO_AMG_UNDECIDED) {             if (split_marker[j]==PASO_AMG_UNDECIDED) {
541            itmp++;            itmp++;
542             } else {  /* at this point there are no C points */             } else {  /* at this point there are no C points */
# Line 557  void Paso_Preconditioner_AMG_RungeStuebe Line 562  void Paso_Preconditioner_AMG_RungeStuebe
562             lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */             lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
563                        
564             /* all undecided unknown strongly coupled to i are moved to F */             /* all undecided unknown strongly coupled to i are moved to F */
565             for (p=0; p<degreeT[i]; ++p) {             for (p=0; p<degree_ST[i]; ++p) {
566            j=ST[offset[i]+p];            j=ST[offset_S[i]+p];
567                        
568            if (split_marker[j]==PASO_AMG_UNDECIDED) {            if (split_marker[j]==PASO_AMG_UNDECIDED) {
569                            
570               split_marker[j]=PASO_AMG_IN_F;               split_marker[j]=PASO_AMG_IN_F;
571               lambda[j]=-1;               lambda[j]=-1;
572                            
573               for (q=0; q<degreeT[j]; ++q) {               for (q=0; q<degree_ST[j]; ++q) {
574              k=ST[offset[j]+q];              k=ST[offset_S[j]+q];
575              if (split_marker[k]==PASO_AMG_UNDECIDED) {              if (split_marker[k]==PASO_AMG_UNDECIDED) {
576                 lambda[k]++;                 lambda[k]++;
577                 if (notInPanel[k]) {                 if (notInPanel[k]) {
# Line 582  void Paso_Preconditioner_AMG_RungeStuebe Line 587  void Paso_Preconditioner_AMG_RungeStuebe
587                            
588            }            }
589             }             }
590             for (p=0; p<degree[i]; ++p) {             for (p=0; p<degree_S[i]; ++p) {
591            j=S[offset[i]+p];            j=S[offset_S[i]+p];
592            if (split_marker[j]==PASO_AMG_UNDECIDED) {            if (split_marker[j]==PASO_AMG_UNDECIDED) {
593               lambda[j]--;               lambda[j]--;
594               if (notInPanel[j]) {               if (notInPanel[j]) {
# Line 622  void Paso_Preconditioner_AMG_RungeStuebe Line 627  void Paso_Preconditioner_AMG_RungeStuebe
627          lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */          lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
628                    
629          /* all undecided unknown strongly coupled to i are moved to F */          /* all undecided unknown strongly coupled to i are moved to F */
630          for (p=0; p<degreeT[i]; ++p) {          for (p=0; p<degree_ST[i]; ++p) {
631             j=ST[offset[i]+p];             j=ST[offset_S[i]+p];
632             if (split_marker[j]==PASO_AMG_UNDECIDED) {             if (split_marker[j]==PASO_AMG_UNDECIDED) {
633                    
634            split_marker[j]=PASO_AMG_IN_F;            split_marker[j]=PASO_AMG_IN_F;
635            lambda[j]=-1;            lambda[j]=-1;
636                    
637            for (q=0; q<degreeT[j]; ++q) {            for (q=0; q<degree_ST[j]; ++q) {
638               k=ST[offset[j]+q];               k=ST[offset_S[j]+q];
639               if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;               if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;
640            }            }
641    
642             }             }
643          }          }
644          for (p=0; p<degree[i]; ++p) {          for (p=0; p<degree_S[i]; ++p) {
645             j=S[offset[i]+p];             j=S[offset_S[i]+p];
646             if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;             if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;
647          }          }
648                    
649       }       }
650       i=Paso_Util_arg_max(n,lambda);       i=Paso_Util_arg_max(n,lambda);
651        }        }
652              
653     }     }
654     TMPMEMFREE(lambda);     TMPMEMFREE(lambda);
655     TMPMEMFREE(ST);     TMPMEMFREE(ST);
656     TMPMEMFREE(degreeT);     TMPMEMFREE(degree_ST);
657     TMPMEMFREE(panel);     TMPMEMFREE(panel);
658     TMPMEMFREE(notInPanel);     TMPMEMFREE(notInPanel);
659  }  }
660    /* ensures that two F nodes are connected via a C node :*/
661    void Paso_Preconditioner_AMG_enforceFFConnectivity(const dim_t n, const index_t* offset_S,
662                            const dim_t* degree_S, const index_t* S,
663                            index_t*split_marker)
664    {
665          dim_t i, p, q;
666    
667          /* now we make sure that two (strongly) connected F nodes are (strongly) connected via a C node. */
668          for (i=0;i<n;++i) {
669                if ( (split_marker[i]==PASO_AMG_IN_F) && (degree_S[i]>0) ) {
670               for (p=0; p<degree_S[i]; ++p) {
671                      register index_t j=S[offset_S[i]+p];
672                      if ( (split_marker[j]==PASO_AMG_IN_F)  && (degree_S[j]>0) )  {
673                          /* i and j are now two F nodes which are strongly connected */
674                          /* is there a C node they share ? */
675                          register index_t sharing=-1;
676                      for (q=0; q<degree_S[i]; ++q) {
677                             index_t k=S[offset_S[i]+q];
678                             if (split_marker[k]==PASO_AMG_IN_C) {
679                                register index_t* where_k=(index_t*)bsearch(&k, &(S[offset_S[j]]), degree_S[j], sizeof(index_t), Paso_comparIndex);
680                                if (where_k != NULL) {
681                                   sharing=k;
682                                   break;
683                                }
684                             }
685                          }
686                          if (sharing<0) {
687                               if (i<j) {
688                                  split_marker[j]=PASO_AMG_IN_C;
689                               } else {
690                                  split_marker[i]=PASO_AMG_IN_C;
691                                  break;  /* no point to look any further as i is now a C node */
692                               }
693                          }
694                      }
695                   }
696                }
697           }
698    }

Legend:
Removed from v.3439  
changed lines
  Added in v.3440

  ViewVC Help
Powered by ViewVC 1.1.26