/[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 3258 by gross, Tue Sep 21 06:56:44 2010 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 143  Paso_Preconditioner_LocalAMG* Paso_Preco Line 143  Paso_Preconditioner_LocalAMG* Paso_Preco
143    /* identify independend set of rows/columns */    /* identify independend set of rows/columns */
144    mis_marker=TMPMEMALLOC(n,index_t);    mis_marker=TMPMEMALLOC(n,index_t);
145    counter=TMPMEMALLOC(n,index_t);    counter=TMPMEMALLOC(n,index_t);
146    if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) {    if ( !( Esys_checkPtr(mis_marker) || Esys_checkPtr(counter) || Esys_checkPtr(out)) ) {
147            
148            
149       out->post_sweeps=options->post_sweeps;       out->post_sweeps=options->post_sweeps;
# Line 210  Paso_Preconditioner_LocalAMG* Paso_Preco Line 210  Paso_Preconditioner_LocalAMG* Paso_Preco
210       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);
211    
212           /* Start Coarsening : */           /* Start Coarsening : */
213           time0=Paso_timer();           time0=Esys_timer();
214       Paso_Coarsening_Local(mis_marker, A_p, options->coarsening_threshold, options->coarsening_method);       Paso_Coarsening_Local(mis_marker, A_p, options->coarsening_threshold, options->coarsening_method);
215           time0=Paso_timer()-time0;           time0=Esys_timer()-time0;
216       if (timing) fprintf(stdout,"Level %d: timing: Coarsening: %e\n",level,time0);       if (timing) fprintf(stdout,"Level %d: timing: Coarsening: %e\n",level,time0);
217    
218          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
# Line 240  Paso_Preconditioner_LocalAMG* Paso_Preco Line 240  Paso_Preconditioner_LocalAMG* Paso_Preco
240              }              }
241          } else {          } else {
242            
243                if (Paso_noError()) {                if (Esys_noError()) {
244                                    
245                   /*#pragma omp parallel for private(i) schedule(static)                   /*#pragma omp parallel for private(i) schedule(static)
246                   for (i = 0; i < n; ++i) counter[i]=mis_marker[i];                   for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
# Line 249  Paso_Preconditioner_LocalAMG* Paso_Preco Line 249  Paso_Preconditioner_LocalAMG* Paso_Preco
249                                    
250                   out->mask_F=MEMALLOC(n,index_t);                   out->mask_F=MEMALLOC(n,index_t);
251                   out->rows_in_F=MEMALLOC(out->n_F,index_t);                   out->rows_in_F=MEMALLOC(out->n_F,index_t);
252                   if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->rows_in_F) ) ) {                   if (! (Esys_checkPtr(out->mask_F) || Esys_checkPtr(out->rows_in_F) ) ) {
253                      /* creates an index for F from mask */                      /* creates an index for F from mask */
254                      #pragma omp parallel for private(i) schedule(static)                      #pragma omp parallel for private(i) schedule(static)
255                      for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;                      for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
# Line 280  Paso_Preconditioner_LocalAMG* Paso_Preco Line 280  Paso_Preconditioner_LocalAMG* Paso_Preco
280                      level=1;                      level=1;
281                  }                  }
282                            
283                if ( Paso_noError()) {                if ( Esys_noError()) {
284                      out->n_C=n-out->n_F;                      out->n_C=n-out->n_F;
285                      out->rows_in_C=MEMALLOC(out->n_C,index_t);                      out->rows_in_C=MEMALLOC(out->n_C,index_t);
286                      out->mask_C=MEMALLOC(n,index_t);                      out->mask_C=MEMALLOC(n,index_t);
287                      if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {                      if (! (Esys_checkPtr(out->mask_C) || Esys_checkPtr(out->rows_in_C) ) ) {
288                           /* creates an index for C from mask */                           /* creates an index for C from mask */
289                           #pragma omp parallel for private(i) schedule(static)                           #pragma omp parallel for private(i) schedule(static)
290                           for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];                           for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
# Line 302  Paso_Preconditioner_LocalAMG* Paso_Preco Line 302  Paso_Preconditioner_LocalAMG* Paso_Preco
302                           }                           }
303                      }                      }
304                }                }
305                if ( Paso_noError()) {                if ( Esys_noError()) {
306                        /* get A_FF block: */                        /* get A_FF block: */
307                        /*                        /*
308                         out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F);                         out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F);
# Line 321  Paso_Preconditioner_LocalAMG* Paso_Preco Line 321  Paso_Preconditioner_LocalAMG* Paso_Preco
321                                    printf("##mis_marker[%d]=%d\n",i,mis_marker[i]);                                    printf("##mis_marker[%d]=%d\n",i,mis_marker[i]);
322                         }                         }
323                        */                                          */                  
324                        time0=Paso_timer();                        time0=Esys_timer();
325                        Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker);                        Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker);
326                        time0=Paso_timer()-time0;                        time0=Esys_timer()-time0;
327                        if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0);                        if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0);
328                    
329                        /*                        /*
# Line 332  Paso_Preconditioner_LocalAMG* Paso_Preco Line 332  Paso_Preconditioner_LocalAMG* Paso_Preco
332                        */                        */
333                                                
334                        /* get Prolongation and Restriction */                        /* get Prolongation and Restriction */
335                        time0=Paso_timer();                        time0=Esys_timer();
336                        out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker);                        out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker);
337                        time0=Paso_timer()-time0;                        time0=Esys_timer()-time0;
338                        if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0);                        if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0);
339                        /*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/                        /*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/
340                                                
# Line 343  Paso_Preconditioner_LocalAMG* Paso_Preco Line 343  Paso_Preconditioner_LocalAMG* Paso_Preco
343                        Paso_SparseMatrix_saveMM(out->P,filename);                        Paso_SparseMatrix_saveMM(out->P,filename);
344                        */                        */
345                                                
346                        time0=Paso_timer();                        time0=Esys_timer();
347                        out->R=Paso_SparseMatrix_getRestriction(out->P);                        out->R=Paso_SparseMatrix_getRestriction(out->P);
348                        time0=Paso_timer()-time0;                        time0=Esys_timer()-time0;
349                        if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0);                        if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0);
350                        /*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/                        /*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/
351                                                
# Line 355  Paso_Preconditioner_LocalAMG* Paso_Preco Line 355  Paso_Preconditioner_LocalAMG* Paso_Preco
355                        */                        */
356                                            
357                }                }
358                if ( Paso_noError()) {                if ( Esys_noError()) {
359                                            
360                      time0=Paso_timer();                      time0=Esys_timer();
361                                            
362                      Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);                      Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
363                                            
# Line 368  Paso_Preconditioner_LocalAMG* Paso_Preco Line 368  Paso_Preconditioner_LocalAMG* Paso_Preco
368                      Paso_SparseMatrix_free(Atemp);                      Paso_SparseMatrix_free(Atemp);
369                                            
370                      /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/                      /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/
371                      time0=Paso_timer()-time0;                      time0=Esys_timer()-time0;
372                      if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0);                      if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0);
373                                            
374                                                                                    
# Line 382  Paso_Preconditioner_LocalAMG* Paso_Preco Line 382  Paso_Preconditioner_LocalAMG* Paso_Preco
382                }                }
383    
384                /* allocate work arrays for AMG application */                /* allocate work arrays for AMG application */
385                if (Paso_noError()) {                if (Esys_noError()) {
386                           /*                           /*
387                            out->x_F=MEMALLOC(n_block*out->n_F,double);                            out->x_F=MEMALLOC(n_block*out->n_F,double);
388                           out->b_F=MEMALLOC(n_block*out->n_F,double);                           out->b_F=MEMALLOC(n_block*out->n_F,double);
# Line 390  Paso_Preconditioner_LocalAMG* Paso_Preco Line 390  Paso_Preconditioner_LocalAMG* Paso_Preco
390                           out->x_C=MEMALLOC(n_block*out->n_C,double);                           out->x_C=MEMALLOC(n_block*out->n_C,double);
391                           out->b_C=MEMALLOC(n_block*out->n_C,double);                           out->b_C=MEMALLOC(n_block*out->n_C,double);
392                
393                           /*if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {*/                           /*if (! (Esys_checkPtr(out->x_F) || Esys_checkPtr(out->b_F) || Esys_checkPtr(out->x_C) || Esys_checkPtr(out->b_C) ) ) {*/
394                           if ( ! ( Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {                           if ( ! ( Esys_checkPtr(out->x_C) || Esys_checkPtr(out->b_C) ) ) {
395                                                            
396                               /*                               /*
397                                #pragma omp parallel for private(i) schedule(static)                                #pragma omp parallel for private(i) schedule(static)
# Line 417  Paso_Preconditioner_LocalAMG* Paso_Preco Line 417  Paso_Preconditioner_LocalAMG* Paso_Preco
417    TMPMEMFREE(mis_marker);    TMPMEMFREE(mis_marker);
418    TMPMEMFREE(counter);    TMPMEMFREE(counter);
419    
420    if (Paso_noError()) {    if (Esys_noError()) {
421        if (verbose && level>0 && !out->coarsest_level) {        if (verbose && level>0 && !out->coarsest_level) {
422           printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);           printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);
423       }       }
# Line 470  void Paso_Preconditioner_LocalAMG_solve( Line 470  void Paso_Preconditioner_LocalAMG_solve(
470            
471       if (amg->coarsest_level) {       if (amg->coarsest_level) {
472                
473        time0=Paso_timer();        time0=Esys_timer();
474        /*If all unknown are eliminated then Jacobi is the best preconditioner*/        /*If all unknown are eliminated then Jacobi is the best preconditioner*/
475    
476        if (amg->n_F==0 || amg->n_F==amg->n) {        if (amg->n_F==0 || amg->n_F==amg->n) {
# Line 490  void Paso_Preconditioner_LocalAMG_solve( Line 490  void Paso_Preconditioner_LocalAMG_solve(
490         #endif         #endif
491         }         }
492                
493         time0=Paso_timer()-time0;         time0=Esys_timer()-time0;
494         if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0);         if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0);
495                
496       } else {       } else {
497          /* presmoothing */          /* presmoothing */
498           time0=Paso_timer();           time0=Esys_timer();
499       Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, pre_sweeps, FALSE);       Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, pre_sweeps, FALSE);
500           time0=Paso_timer()-time0;           time0=Esys_timer()-time0;
501           if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0);           if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0);
502            
503           /* end of presmoothing */           /* end of presmoothing */
504                    
505           time0=Paso_timer();           time0=Esys_timer();
506            #pragma omp parallel for private(i,j) schedule(static)            #pragma omp parallel for private(i,j) schedule(static)
507             for (i=0;i<amg->n;++i) {             for (i=0;i<amg->n;++i) {
508              for (j=0;j<amg->n_block;++j) {              for (j=0;j<amg->n_block;++j) {
# Line 516  void Paso_Preconditioner_LocalAMG_solve( Line 516  void Paso_Preconditioner_LocalAMG_solve(
516          /* b_c = R*r  */          /* b_c = R*r  */
517           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C);           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C);
518                    
519          time0=Paso_timer()-time0;          time0=Esys_timer()-time0;
520          if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0);          if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0);
521                    
522          /* x_C=AMG(b_C)     */          /* x_C=AMG(b_C)     */
523      Paso_Preconditioner_LocalAMG_solve(amg->AMG_of_Coarse,amg->x_C,amg->b_C);      Paso_Preconditioner_LocalAMG_solve(amg->AMG_of_Coarse,amg->x_C,amg->b_C);
524                    
525          time0=Paso_timer();          time0=Esys_timer();
526                    
527          /* x_0 = P*x_c */          /* x_0 = P*x_c */
528          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0);          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0);
# Line 538  void Paso_Preconditioner_LocalAMG_solve( Line 538  void Paso_Preconditioner_LocalAMG_solve(
538        /*postsmoothing*/        /*postsmoothing*/
539                
540        /*solve Ax=b with initial guess x */        /*solve Ax=b with initial guess x */
541        time0=Paso_timer();        time0=Esys_timer();
542        Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, post_sweeps, TRUE);        Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, post_sweeps, TRUE);
543        time0=Paso_timer()-time0;        time0=Esys_timer()-time0;
544        if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);        if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);
545    
546        /*end of postsmoothing*/        /*end of postsmoothing*/

Legend:
Removed from v.3258  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26