/[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 3885 by gross, Wed Apr 4 22:12:26 2012 UTC revision 3886 by gross, Thu Apr 5 00:50:30 2012 UTC
# Line 101  Paso_Preconditioner_AMG* Paso_Preconditi Line 101  Paso_Preconditioner_AMG* Paso_Preconditi
101    
102    const dim_t n_block=A_p->row_block_size;    const dim_t n_block=A_p->row_block_size;
103    index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F;    index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F;
104    dim_t i, n_F, n_C, F_flag, *F_set=NULL, global_n_C=0, global_n_F=0;    dim_t i, my_n_F, my_n_C, n_C, F_flag, *F_set=NULL, global_n_C=0, global_n_F=0, n_F;
105    double time0=0;    double time0=0;
106    const double theta = options->coarsening_threshold;    const double theta = options->coarsening_threshold;
107    const double tau = options->diagonal_dominance_threshold;    const double tau = options->diagonal_dominance_threshold;
# Line 201  Paso_Preconditioner_AMG* Paso_Preconditi Line 201  Paso_Preconditioner_AMG* Paso_Preconditi
201          /*          /*
202             count number of unkowns to be eliminated:             count number of unkowns to be eliminated:
203          */          */
204          n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);          my_n_F=Paso_Util_cumsum_maskedTrue(my_n,counter, F_marker);
205              /* collect n_F values on all processes, a direct solver should          n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);
206                  be used if any n_F value is 0 */              /* collect my_n_F values on all processes, a direct solver should
207                    be used if any my_n_F value is 0 */
208              F_set = TMPMEMALLOC(A_p->mpi_info->size, dim_t);              F_set = TMPMEMALLOC(A_p->mpi_info->size, dim_t);
209          #ifdef ESYS_MPI          #ifdef ESYS_MPI
210              MPI_Allgather(&n_F, 1, MPI_INT, F_set, 1, MPI_INT, A_p->mpi_info->comm);              MPI_Allgather(&my_n_F, 1, MPI_INT, F_set, 1, MPI_INT, A_p->mpi_info->comm);
211          #endif          #endif
212              global_n_F=0;              global_n_F=0;
213              F_flag = 1;              F_flag = 1;
214              for (i=0; i<A_p->mpi_info->size; i++) {              for (i=0; i<A_p->mpi_info->size; i++) {
215                  global_n_F+=F_set[i];                  global_n_F+=F_set[i];
216                  if (F_set[i] == 0) F_flag = 0;                  if (F_set[i] == 0) F_flag = 0;
 printf(" p [%d]=%d\n",i, F_set[i]);  
217              }              }
 printf(" global n = %d\n",global_n);  
218              TMPMEMFREE(F_set);              TMPMEMFREE(F_set);
219    
220          n_C=n-n_F;          my_n_C=my_n-my_n_F;
221              global_n_C=global_n-global_n_F;              global_n_C=global_n-global_n_F;
222          if (verbose) printf("Paso_Preconditioner: AMG (non-local) level %d: %d unknowns are flagged for elimination. %d left.\n",level,global_n_F,global_n_C);          if (verbose) printf("Paso_Preconditioner: AMG (non-local) level %d: %d unknowns are flagged for elimination. %d left.\n",level,global_n_F,global_n_C);
223    
# Line 230  printf(" global n = %d\n",global_n); Line 229  printf(" global n = %d\n",global_n);
229             out=MEMALLOC(1,Paso_Preconditioner_AMG);             out=MEMALLOC(1,Paso_Preconditioner_AMG);
230             if (! Esys_checkPtr(out)) {             if (! Esys_checkPtr(out)) {
231            out->level = level;            out->level = level;
           out->n = n;  
           out->n_F = n_F;  
           out->n_block = n_block;  
232            out->A_C = NULL;            out->A_C = NULL;
233            out->P = NULL;              out->P = NULL;  
234            out->R = NULL;                      out->R = NULL;          
# Line 253  printf(" global n = %d\n",global_n); Line 249  printf(" global n = %d\n",global_n);
249            out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), 0, verbose);            out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), 0, verbose);
250                
251            if (global_n_C != 0) {            if (global_n_C != 0) {
252                 /* if nothing has been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */              /*  create mask of C nodes with value >-1, gives new id */
253                n_C=Paso_Util_cumsum_maskedFalse(n, mask_C, F_marker);
254                /* if nothing has been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
255        
256              /* allocate helpers :*/              /* allocate helpers :*/
257              out->x_C=MEMALLOC(n_block*n_C,double);              out->x_C=MEMALLOC(n_block*my_n_C,double);
258              out->b_C=MEMALLOC(n_block*n_C,double);              out->b_C=MEMALLOC(n_block*my_n_C,double);
259              out->r=MEMALLOC(n_block*n,double);              out->r=MEMALLOC(n_block*my_n,double);
260                            
261              Esys_checkPtr(out->r);              Esys_checkPtr(out->r);
262              Esys_checkPtr(out->x_C);              Esys_checkPtr(out->x_C);
# Line 273  printf(" global n = %d\n",global_n); Line 271  printf(" global n = %d\n",global_n);
271                   if  (F_marker[i]) rows_in_F[counter[i]]=i;                   if  (F_marker[i]) rows_in_F[counter[i]]=i;
272                    }                    }
273                 }                 }
                /*  create mask of C nodes with value >-1, gives new id */  
                i=Paso_Util_cumsum_maskedFalse(n, mask_C, F_marker);  
274                 /*                 /*
275                    get Prolongation :                        get Prolongation :    
276                 */                                   */                  
# Line 372  void Paso_Preconditioner_AMG_solve(Paso_ Line 368  void Paso_Preconditioner_AMG_solve(Paso_
368       if (SHOW_TIMING) printf("timing: level %d: Presmoothing: %e\n",amg->level, time0);       if (SHOW_TIMING) printf("timing: level %d: Presmoothing: %e\n",amg->level, time0);
369       /* end of presmoothing */       /* end of presmoothing */
370            
371       if (amg->n_F < amg->n) { /* is there work on the coarse level? */       time0=Esys_timer();
372           time0=Esys_timer();  
373         Paso_Copy(n, amg->r, b);                            /*  r <- b */
374         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
375         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */
376    
377       Paso_Copy(n, amg->r, b);                            /*  r <- b */       time0=Esys_timer()-time0;
378       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/       /* coarse level solve */
379       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */       if ( amg->AMG_C == NULL) {
   
          time0=Esys_timer()-time0;  
      /* coarse level solve */  
      if ( amg->AMG_C == NULL) {  
380          time0=Esys_timer();          time0=Esys_timer();
381          /*  A_C is the coarsest level */          /*  A_C is the coarsest level */
382          Paso_Preconditioner_AMG_mergeSolve(amg);          Paso_Preconditioner_AMG_mergeSolve(amg);
383    
384          if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);          if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
385       } else {       } else {
386          Paso_Preconditioner_AMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */          Paso_Preconditioner_AMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */
387       }       }
388    
389       time0=time0+Esys_timer();      time0=time0+Esys_timer();
390       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */          Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */    
391    
392           /*postsmoothing*/      /*postsmoothing*/
393                
394          /*solve Ax=b with initial guess x */      /*solve Ax=b with initial guess x */
395          time0=Esys_timer();      time0=Esys_timer();
396          Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);      Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
397          time0=Esys_timer()-time0;      time0=Esys_timer()-time0;
398          if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);      if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
399          /*end of postsmoothing*/      return;
      }  
   
      return;  
400  }  }
401    
402  /* theta = threshold for strong connections */  /* theta = threshold for strong connections */
# Line 1117  void Paso_Preconditioner_AMG_mergeSolve( Line 1109  void Paso_Preconditioner_AMG_mergeSolve(
1109    Paso_SparseMatrix *A_D, *A_temp;    Paso_SparseMatrix *A_D, *A_temp;
1110    double* x=NULL;    double* x=NULL;
1111    double* b=NULL;    double* b=NULL;
1112    index_t rank = A->mpi_info->rank;    const index_t rank = A->mpi_info->rank;
1113    index_t size = A->mpi_info->size;    const index_t size = A->mpi_info->size;
1114    index_t i, n, p, n_block;    const n_block = A->mainBlock->row_block_size;
1115    
1116      index_t i, n, p;
1117    index_t *counts, *offset, *dist;    index_t *counts, *offset, *dist;
1118    #ifdef ESYS_MPI    #ifdef ESYS_MPI
1119    index_t count;    index_t count;
1120    #endif    #endif
   n_block = amg->n_block;  
1121    A_D = Paso_Preconditioner_AMG_mergeSystemMatrix(A);    A_D = Paso_Preconditioner_AMG_mergeSystemMatrix(A);
1122    
1123    /* First, gather x and b into rank 0 */    /* First, gather x and b into rank 0 */

Legend:
Removed from v.3885  
changed lines
  Added in v.3886

  ViewVC Help
Powered by ViewVC 1.1.26