/[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 3489 by caltinay, Wed Mar 30 00:46:04 2011 UTC revision 3642 by caltinay, Thu Oct 27 03:41:51 2011 UTC
# Line 18  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /* Author: artak@uq.edu.au, l.gross@uq.edu.au                                */  /* Author: artak@uq.edu.au, l.gross@uq.edu.au                 */
22    
23  /**************************************************************/  /**************************************************************/
24    
# Line 107  Paso_Preconditioner_AMG* Paso_Preconditi Line 107  Paso_Preconditioner_AMG* Paso_Preconditi
107    
108        
109    /*    /*
110        is the input matrix A suitable for coarsening        is the input matrix A suitable for coarsening?
111                
112    */    */
113    if ( (sparsity >= options->min_coarse_sparsity) ||    if ( (sparsity >= options->min_coarse_sparsity) ||
# Line 124  Paso_Preconditioner_AMG* Paso_Preconditi Line 124  Paso_Preconditioner_AMG* Paso_Preconditi
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    
# Line 196  return NULL; Line 196  return NULL;
196       Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);       Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);
197  */  */
198            
199           /* in BoomerAMG interpolation is used FF connectiovity is required :*/           /* in BoomerAMG if interpolation is used FF connectivity is required */
200  /*MPI:  /*MPI:
201           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
202                               Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);                                 Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);  
# Line 216  return NULL; Line 216  return NULL;
216          n_C=n-n_F;          n_C=n-n_F;
217          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);
218            
219          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 */
220             out = NULL;             out = NULL;
221          } else {          } else {
222             out=MEMALLOC(1,Paso_Preconditioner_AMG);             out=MEMALLOC(1,Paso_Preconditioner_AMG);
# Line 245  return NULL; Line 245  return NULL;
245            out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);            out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
246                
247            if (n_C != 0) {            if (n_C != 0) {
248                 /* 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 dominant matrix and we just run a few steps of the smoother */
249        
250              /* allocate helpers :*/              /* allocate helpers :*/
251              out->x_C=MEMALLOC(n_block*n_C,double);              out->x_C=MEMALLOC(n_block*n_C,double);
# Line 265  return NULL; Line 265  return NULL;
265                   if  (F_marker[i]) rows_in_F[counter[i]]=i;                   if  (F_marker[i]) rows_in_F[counter[i]]=i;
266                    }                    }
267                 }                 }
268                 /*  create mask of C nodes with value >-1 gives new id */                 /*  create mask of C nodes with value >-1, gives new id */
269                 i=Paso_Util_cumsum_maskedFalse(n,counter, F_marker);                 i=Paso_Util_cumsum_maskedFalse(n,counter, F_marker);
270    
271                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
# Line 279  return NULL; Line 279  return NULL;
279                 /*                 /*
280                    get Prolongation :                        get Prolongation :    
281                 */                                   */                  
                time0=Esys_timer();  
282  /*MPI:  /*MPI:
283                   time0=Esys_timer();
284                 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);                 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);
 */  
285                 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);
286    */
287              }              }
288              /*                    /*      
289                 construct Restriction operator as transposed of Prolongation operator:                 construct Restriction operator as transposed of Prolongation operator:
290              */              */
291    /*MPI:
292              if ( Esys_noError()) {              if ( Esys_noError()) {
293                 time0=Esys_timer();                 time0=Esys_timer();
 /*MPI:  
294                 out->R=Paso_SystemMatrix_getTranspose(out->P);                 out->R=Paso_SystemMatrix_getTranspose(out->P);
 */  
295                 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);                 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
296              }                    }      
297    */
298              /*              /*
299              construct coarse level matrix:              construct coarse level matrix:
300              */              */
301    /*MPI:
302              if ( Esys_noError()) {              if ( Esys_noError()) {
303                 time0=Esys_timer();                 time0=Esys_timer();
 /*MPI:  
304                 Atemp=Paso_SystemMatrix_MatrixMatrix(A_p,out->P);                 Atemp=Paso_SystemMatrix_MatrixMatrix(A_p,out->P);
305                 A_C=Paso_SystemMatrix_MatrixMatrix(out->R,Atemp);                 A_C=Paso_SystemMatrix_MatrixMatrix(out->R,Atemp);
306                 Paso_Preconditioner_AMG_setStrongConnections                 Paso_Preconditioner_AMG_setStrongConnections
307                 Paso_SystemMatrix_free(Atemp);                 Paso_SystemMatrix_free(Atemp);
 */  
308    
309                 if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);                             if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);            
310              }              }
311    */
312    
313                            
314              /*              /*
# Line 385  void Paso_Preconditioner_AMG_solve(Paso_ Line 385  void Paso_Preconditioner_AMG_solve(Paso_
385       time0=Esys_timer();       time0=Esys_timer();
386       Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);       Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
387       time0=Esys_timer()-time0;       time0=Esys_timer()-time0;
388       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);
389       /* end of presmoothing */       /* end of presmoothing */
390            
391       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 489  void Paso_Preconditioner_AMG_setStrongCo Line 489  void Paso_Preconditioner_AMG_setStrongCo
489           {           {
490          const double threshold = theta*max_offdiagonal;          const double threshold = theta*max_offdiagonal;
491              threshold_p[2*i+1]=threshold;              threshold_p[2*i+1]=threshold;
492          if (tau*main_row < sum_row) { /* no diagonal domainance */          if (tau*main_row < sum_row) { /* no diagonal dominance */
493                 threshold_p[2*i]=1;                 threshold_p[2*i]=1;
494             #pragma ivdep             #pragma ivdep
495             for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {             for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
# Line 638  void Paso_Preconditioner_AMG_setStrongCo Line 638  void Paso_Preconditioner_AMG_setStrongCo
638          rtmp_offset=-A->mainBlock->pattern->ptr[i];          rtmp_offset=-A->mainBlock->pattern->ptr[i];
639                    
640          threshold_p[2*i+1]=threshold;          threshold_p[2*i+1]=threshold;
641          if (tau*main_row < sum_row) { /* no diagonal domainance */          if (tau*main_row < sum_row) { /* no diagonal dominance */
642             threshold_p[2*i]=1;             threshold_p[2*i]=1;
643             rtmp_offset=-A->mainBlock->pattern->ptr[i];             rtmp_offset=-A->mainBlock->pattern->ptr[i];
644             #pragma ivdep             #pragma ivdep
# Line 711  void Paso_Preconditioner_AMG_setStrongCo Line 711  void Paso_Preconditioner_AMG_setStrongCo
711     }     }
712     TMPMEMFREE(threshold_p);     TMPMEMFREE(threshold_p);
713  }  }
714    
715  void Paso_Preconditioner_AMG_transposeStrongConnections(const dim_t n, const dim_t* degree_S, const index_t* offset_S, const index_t* S,  void Paso_Preconditioner_AMG_transposeStrongConnections(const dim_t n, const dim_t* degree_S, const index_t* offset_S, const index_t* S,
716                              const dim_t nT, dim_t* degree_ST, index_t* offset_ST,index_t* ST)                              const dim_t nT, dim_t* degree_ST, index_t* offset_ST,index_t* ST)
717  {  {
# Line 791  void Paso_Preconditioner_AMG_CIJPCoarsen Line 792  void Paso_Preconditioner_AMG_CIJPCoarsen
792            
793       }       }
794            
795        /* calculate the maximum value of naigbours following active strong connections:        /* calculate the maximum value of neighbours following active strong connections:
796          w2[i]=MAX(w[k]) with k in ST[i] or S[i] and (i,k) conenction is still active  */                w2[i]=MAX(w[k]) with k in ST[i] or S[i] and (i,k) connection is still active  */      
797        #pragma omp parallel for private(i, iptr)        #pragma omp parallel for private(i, iptr)
798        for (i=0; i<my_n; ++i) {        for (i=0; i<my_n; ++i) {
799       if (Status[i]>0) { /* status is still undefined */       if (Status[i]>0) { /* status is still undefined */
# Line 896  printf("found!\n"); Line 897  printf("found!\n");
897                    if (ST_flag[offset_ST[j]+kptr] >0) {                    if (ST_flag[offset_ST[j]+kptr] >0) {
898                   if (j< my_n ) {                   if (j< my_n ) {
899                      w[j]--;                      w[j]--;
900  printf("%d reduced by %d and %d \n",j, i,k);  printf("%d reduced by %d and %d\n",j, i,k);
901                   }                   }
902                   ST_flag[offset_ST[j]+kptr]=-1;                   ST_flag[offset_ST[j]+kptr]=-1;
903                    }                    }

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

  ViewVC Help
Powered by ViewVC 1.1.26