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

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

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

revision 3641 by lgao, Thu Apr 7 01:12:59 2011 UTC revision 3642 by caltinay, Thu Oct 27 03:41:51 2011 UTC
# Line 60  index_t Paso_Preconditioner_LocalAMG_get Line 60  index_t Paso_Preconditioner_LocalAMG_get
60        return Paso_Preconditioner_LocalAMG_getMaxLevel(in->AMG_C);        return Paso_Preconditioner_LocalAMG_getMaxLevel(in->AMG_C);
61     }     }
62  }  }
63    
64  double Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(const Paso_Preconditioner_LocalAMG * in) {  double Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(const Paso_Preconditioner_LocalAMG * in) {
65        if (in->AMG_C == NULL) {        if (in->AMG_C == NULL) {
66       if (in->A_C == NULL) {       if (in->A_C == NULL) {
# Line 83  dim_t Paso_Preconditioner_LocalAMG_getNu Line 84  dim_t Paso_Preconditioner_LocalAMG_getNu
84       return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);       return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);
85     }     }
86  }  }
87    
88  /*****************************************************************  /*****************************************************************
89    
90     constructs AMG     Constructs AMG
91        
92  ******************************************************************/  ******************************************************************/
93  Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {  Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
# Line 122  Paso_Preconditioner_LocalAMG* Paso_Preco Line 124  Paso_Preconditioner_LocalAMG* Paso_Preco
124      printf("Paso_Preconditioner: AMG: termination of coarsening by ");      printf("Paso_Preconditioner: AMG: termination of coarsening by ");
125            
126      if (sparsity >= options->min_coarse_sparsity)      if (sparsity >= options->min_coarse_sparsity)
127         printf("SPAR ");         printf("SPAR");
128            
129      if (total_n <= options->min_coarse_matrix_size)      if (total_n <= options->min_coarse_matrix_size)
130         printf("SIZE ");         printf("SIZE");
131            
132      if (level > options->level_max)      if (level > options->level_max)
133         printf("LEVEL ");         printf("LEVEL");
134            
135      printf("\n");      printf("\n");
136            
137      printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",      printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. Sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
138             level,  options->level_max, sparsity, options->min_coarse_sparsity, total_n, options->min_coarse_matrix_size);               level,  options->level_max, sparsity, options->min_coarse_sparsity, total_n, options->min_coarse_matrix_size);  
139                        
140       }       }
# Line 147  Paso_Preconditioner_LocalAMG* Paso_Preco Line 149  Paso_Preconditioner_LocalAMG* Paso_Preco
149       S=TMPMEMALLOC(A_p->pattern->len, index_t);       S=TMPMEMALLOC(A_p->pattern->len, index_t);
150       if ( !( Esys_checkPtr(F_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree_S) || Esys_checkPtr(S) ) ) {       if ( !( Esys_checkPtr(F_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree_S) || Esys_checkPtr(S) ) ) {
151       /*       /*
152            set splitting of unknows:            set splitting of unknowns:
         
153           */           */
154       time0=Esys_timer();       time0=Esys_timer();
155       if (n_block>1) {       if (n_block>1) {
# Line 159  Paso_Preconditioner_LocalAMG* Paso_Preco Line 160  Paso_Preconditioner_LocalAMG* Paso_Preco
160    
161       Paso_Preconditioner_LocalAMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);       Paso_Preconditioner_LocalAMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);
162    
163           /* in BoomerAMG interpolation is used FF connectiovity is required :*/           /* in BoomerAMG if interpolation is used FF connectivity is required: */
164           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
165                               Paso_Preconditioner_LocalAMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);                                 Paso_Preconditioner_LocalAMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);  
166       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);
# Line 169  Paso_Preconditioner_LocalAMG* Paso_Preco Line 170  Paso_Preconditioner_LocalAMG* Paso_Preco
170          for (i = 0; i < n; ++i) F_marker[i]=(F_marker[i] ==  PASO_AMG_IN_F);          for (i = 0; i < n; ++i) F_marker[i]=(F_marker[i] ==  PASO_AMG_IN_F);
171            
172          /*          /*
173             count number of unkowns to be eliminated:             count number of unknowns to be eliminated:
174          */          */
175          n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);          n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);
176          n_C=n-n_F;          n_C=n-n_F;
177          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);
178            
179          if ( n_F == 0 ) {  /*  is a nasty case. a direct solver should be used, return NULL */          if ( n_F == 0 ) {  /* This is a nasty case. A direct solver should be used, return NULL */
180             out = NULL;             out = NULL;
181          } else {          } else {
182             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);
# Line 204  Paso_Preconditioner_LocalAMG* Paso_Preco Line 205  Paso_Preconditioner_LocalAMG* Paso_Preco
205            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
206                
207            if (n_C != 0) {            if (n_C != 0) {
208                 /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */                 /* if nothing has been removed we have a diagonal
209                                * dominant matrix and we just run a few steps of
210                                * the smoother */
211        
212              /* allocate helpers :*/              /* allocate helpers :*/
213              out->x_C=MEMALLOC(n_block*n_C,double);              out->x_C=MEMALLOC(n_block*n_C,double);
# Line 216  Paso_Preconditioner_LocalAMG* Paso_Preco Line 219  Paso_Preconditioner_LocalAMG* Paso_Preco
219              Esys_checkPtr(out->b_C);              Esys_checkPtr(out->b_C);
220                            
221              if ( Esys_noError() ) {              if ( Esys_noError() ) {
222                 /* creates index for F:*/                 /* creates index for F */
223                 #pragma omp parallel private(i)                 #pragma omp parallel private(i)
224                 {                 {
225                    #pragma omp for schedule(static)                    #pragma omp for schedule(static)
# Line 266  Paso_Preconditioner_LocalAMG* Paso_Preco Line 269  Paso_Preconditioner_LocalAMG* Paso_Preco
269    
270                            
271              /*              /*
272                 constructe courser level:                 construct coarser level:
273                                
274              */              */
275              if ( Esys_noError()) {              if ( Esys_noError()) {
# Line 276  Paso_Preconditioner_LocalAMG* Paso_Preco Line 279  Paso_Preconditioner_LocalAMG* Paso_Preco
279                 if ( out->AMG_C == NULL ) {                 if ( out->AMG_C == NULL ) {
280                    out->reordering = options->reordering;                    out->reordering = options->reordering;
281                    out->refinements = options->coarse_matrix_refinements;                    out->refinements = options->coarse_matrix_refinements;
282                    /* no coarse level matrix has been constructed. use direct solver */                    /* no coarse level matrix has been constructed. Use direct solver */
283                    #ifdef MKL                    #ifdef MKL
284                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
285                      Paso_SparseMatrix_free(A_C);                      Paso_SparseMatrix_free(A_C);
# Line 331  void Paso_Preconditioner_LocalAMG_solve( Line 334  void Paso_Preconditioner_LocalAMG_solve(
334       time0=Esys_timer();       time0=Esys_timer();
335       Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);       Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
336       time0=Esys_timer()-time0;       time0=Esys_timer()-time0;
337       if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);       if (SHOW_TIMING) printf("timing: level %d: Presmoothing: %e\n",amg->level, time0);
338       /* end of presmoothing */       /* end of presmoothing */
339            
340       if (amg->n_F < amg->n) { /* is there work on the coarse level? */       if (amg->n_F < amg->n) { /* is there work on the coarse level? */
# Line 411  void Paso_Preconditioner_LocalAMG_setStr Line 414  void Paso_Preconditioner_LocalAMG_setStr
414       {       {
415          const double threshold = theta*max_offdiagonal;          const double threshold = theta*max_offdiagonal;
416          register dim_t kdeg=0;          register dim_t kdeg=0;
417          if (tau*main_row < sum_row) { /* no diagonal domainance */          if (tau*main_row < sum_row) { /* no diagonal dominance */
418             #pragma ivdep             #pragma ivdep
419             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) {
420                register index_t j=A->pattern->index[iptr];                register index_t j=A->pattern->index[iptr];
# Line 481  void Paso_Preconditioner_LocalAMG_setStr Line 484  void Paso_Preconditioner_LocalAMG_setStr
484              {              {
485             const double threshold = theta*max_offdiagonal;             const double threshold = theta*max_offdiagonal;
486             register dim_t kdeg=0;             register dim_t kdeg=0;
487             if (tau*main_row < sum_row) { /* no diagonal domainance */             if (tau*main_row < sum_row) { /* no diagonal dominance */
488                #pragma ivdep                #pragma ivdep
489                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) {
490               register index_t j=A->pattern->index[iptr];               register index_t j=A->pattern->index[iptr];
# Line 499  void Paso_Preconditioner_LocalAMG_setStr Line 502  void Paso_Preconditioner_LocalAMG_setStr
502    
503  }    }  
504    
505  /* the runge stueben coarsening algorithm: */  /* the Runge Stueben coarsening algorithm: */
506  void Paso_Preconditioner_LocalAMG_RungeStuebenSearch(const dim_t n, const index_t* offset_S,  void Paso_Preconditioner_LocalAMG_RungeStuebenSearch(const dim_t n, const index_t* offset_S,
507                          const dim_t* degree_S, const index_t* S,                          const dim_t* degree_S, const index_t* S,
508                          index_t*split_marker, const bool_t usePanel)                          index_t*split_marker, const bool_t usePanel)
# Line 524  void Paso_Preconditioner_LocalAMG_RungeS Line 527  void Paso_Preconditioner_LocalAMG_RungeS
527        
528        
529     if (Esys_noError() ) {     if (Esys_noError() ) {
530        /* initialize  split_marker and split_marker :*/        /* initialize split_marker: */
531        /* those unknows which are not influenced go into F, the rest is available for F or C */        /* Those unknowns which are not influenced go into F, the rest is available for F or C */
532        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
533        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
534       degree_ST[i]=0;       degree_ST[i]=0;
# Line 576  void Paso_Preconditioner_LocalAMG_RungeS Line 579  void Paso_Preconditioner_LocalAMG_RungeS
579                 split_marker[i]=PASO_AMG_IN_C;                 split_marker[i]=PASO_AMG_IN_C;
580                 lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */                 lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
581    
582                 /* all undecided unknown strongly coupled to i are moved to F */                 /* all undecided unknowns strongly coupled to i are moved to F */
583                 for (p=0; p<degree_ST[i]; ++p) {                 for (p=0; p<degree_ST[i]; ++p) {
584                    j=ST[offset_S[i]+p];                    j=ST[offset_S[i]+p];
585    
# Line 622  void Paso_Preconditioner_LocalAMG_RungeS Line 625  void Paso_Preconditioner_LocalAMG_RungeS
625             split_marker[i]=PASO_AMG_IN_C;             split_marker[i]=PASO_AMG_IN_C;
626             lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */             lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
627                        
628             /* all undecided unknown strongly coupled to i are moved to F */             /* all undecided unknowns strongly coupled to i are moved to F */
629             for (p=0; p<degree_ST[i]; ++p) {             for (p=0; p<degree_ST[i]; ++p) {
630            j=ST[offset_S[i]+p];            j=ST[offset_S[i]+p];
631                        
# Line 687  void Paso_Preconditioner_LocalAMG_RungeS Line 690  void Paso_Preconditioner_LocalAMG_RungeS
690          split_marker[i]=PASO_AMG_IN_C;          split_marker[i]=PASO_AMG_IN_C;
691          lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */          lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
692                    
693          /* all undecided unknown strongly coupled to i are moved to F */          /* all undecided unknowns strongly coupled to i are moved to F */
694          for (p=0; p<degree_ST[i]; ++p) {          for (p=0; p<degree_ST[i]; ++p) {
695             j=ST[offset_S[i]+p];             j=ST[offset_S[i]+p];
696             if (split_marker[j]==PASO_AMG_UNDECIDED) {             if (split_marker[j]==PASO_AMG_UNDECIDED) {
# Line 720  void Paso_Preconditioner_LocalAMG_RungeS Line 723  void Paso_Preconditioner_LocalAMG_RungeS
723       if (!SMALL_PANEL) TMPMEMFREE(notInPanel);       if (!SMALL_PANEL) TMPMEMFREE(notInPanel);
724     }     }
725  }  }
726    
727  /* ensures that two F nodes are connected via a C node :*/  /* ensures that two F nodes are connected via a C node :*/
728  void Paso_Preconditioner_LocalAMG_enforceFFConnectivity(const dim_t n, const index_t* offset_S,  void Paso_Preconditioner_LocalAMG_enforceFFConnectivity(const dim_t n, const index_t* offset_S,
729                          const dim_t* degree_S, const index_t* S,                          const dim_t* degree_S, const index_t* S,
# Line 759  void Paso_Preconditioner_LocalAMG_enforc Line 763  void Paso_Preconditioner_LocalAMG_enforc
763              }              }
764         }         }
765  }  }
766    

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

  ViewVC Help
Powered by ViewVC 1.1.26