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

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

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

revision 2827 by artak, Thu Dec 10 02:09:43 2009 UTC revision 2828 by artak, Tue Dec 22 01:24:40 2009 UTC
# Line 51  void Paso_Solver_AMG_System_free(Paso_So Line 51  void Paso_Solver_AMG_System_free(Paso_So
51    
52  void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {  void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {
53       if (in!=NULL) {       if (in!=NULL) {
54          Paso_Solver_Jacobi_free(in->GS);          
55            if(in->Smoother->ID==PASO_JACOBI)
56                Paso_Solver_Jacobi_free(in->Smoother->Jacobi);
57            else if (in->Smoother->ID==PASO_GS)    
58                Paso_Solver_GS_free(in->Smoother->GS);
59            MEMFREE(in->Smoother);
60                
61          Paso_SparseMatrix_free(in->A_FC);          Paso_SparseMatrix_free(in->A_FC);
62          Paso_SparseMatrix_free(in->A_FF);          Paso_SparseMatrix_free(in->A_FF);
63          Paso_SparseMatrix_free(in->W_FC);          Paso_SparseMatrix_free(in->W_FC);
# Line 142  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 148  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
148       return NULL;       return NULL;
149    }    }
150    out=MEMALLOC(1,Paso_Solver_AMG);    out=MEMALLOC(1,Paso_Solver_AMG);
151      out->Smoother=MEMALLOC(1,Paso_Solver_Smoother);
152    /* identify independend set of rows/columns */    /* identify independend set of rows/columns */
153    mis_marker=TMPMEMALLOC(n,index_t);    mis_marker=TMPMEMALLOC(n,index_t);
154    counter=TMPMEMALLOC(n,index_t);    counter=TMPMEMALLOC(n,index_t);
# Line 161  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 168  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
168       out->b_F=NULL;       out->b_F=NULL;
169       out->x_C=NULL;       out->x_C=NULL;
170       out->b_C=NULL;       out->b_C=NULL;
      out->GS=NULL;  
171       out->A=Paso_SparseMatrix_getReference(A_p);       out->A=Paso_SparseMatrix_getReference(A_p);
172       out->solver=NULL;       out->solver=NULL;
173         out->Smoother->ID=options->smoother;
174         out->Smoother->Jacobi=NULL;
175         out->Smoother->GS=NULL;
176       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/
177       out->level=level;       out->level=level;
178       out->n=n;       out->n=n;
# Line 195  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 204  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
204           #else           #else
205              #ifdef UMFPACK              #ifdef UMFPACK
206              #else              #else
207                  out->GS=Paso_Solver_getJacobi(A_p);                if (options->smoother == PASO_JACOBI)
208                    out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
209                  else if (options->smoother == PASO_GS)
210                    out->Smoother->GS=Paso_Solver_getGS(A_p,verbose));
211              #endif              #endif
212           #endif           #endif
213                    
214       } else {       } else {
215           out->coarsest_level=FALSE;           out->coarsest_level=FALSE;
216           out->GS=Paso_Solver_getJacobi(A_p);          
217            if (options->smoother == PASO_JACOBI)
218                    out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
219            else if (options->smoother == PASO_GS)
220                    out->Smoother->GS=Paso_Solver_getGS(A_p,verbose);
221    
222           /* identify independend set of rows/columns */           /* identify independend set of rows/columns */
223           #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
# Line 224  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 240  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
240          }          }
241          else {          else {
242             /*Default coarseneing*/             /*Default coarseneing*/
             /*Paso_Pattern_RS_MI_Aggressive(A_p,mis_marker,options->coarsening_threshold);*/  
243              Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold);              Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold);
244              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/
245              /*Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_greedy_Agg(A_p,mis_marker,options->coarsening_threshold);*/
246              /*Paso_Pattern_greedy(A_p->pattern,mis_marker);*/              /*Paso_Pattern_greedy(A_p->pattern,mis_marker);*/
247              /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
248                            
# Line 286  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 301  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
301                                            
302                   }                   }
303                }                }
304                                /*
305                /* if(level==1) {                 if(level==1) {
306                     printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F);                     printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F);
307                     for (i = 0; i < n; ++i) {                     for (i = 0; i < n; ++i) {
308                      printf("##%d %d\n",i,mis_marker[i]);                      printf("##%d %d\n",i,mis_marker[i]);
# Line 359  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 374  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
374                        if (verbose) fprintf(stdout,"timing: getProlongation: %e\n",time0);                        if (verbose) fprintf(stdout,"timing: getProlongation: %e\n",time0);
375                                                
376                                                
377                                                /*
378                        /* sprintf(filename,"P_%d",level);                         sprintf(filename,"P_%d",level);
379                        Paso_SparseMatrix_saveMM(out->P,filename);                        Paso_SparseMatrix_saveMM(out->P,filename);
380                        */                        */
381                                                
# Line 492  void Paso_Solver_solveAMG(Paso_Solver_AM Line 507  void Paso_Solver_solveAMG(Paso_Solver_AM
507        /*Paso_Solver_solveJacobi(amg->GS,x,b);*/        /*Paso_Solver_solveJacobi(amg->GS,x,b);*/
508                
509        if (amg->n_F==0 || amg->n_F==amg->n) {        if (amg->n_F==0 || amg->n_F==amg->n) {
510          Paso_Solver_solveJacobi(amg->GS,x,b);          if(amg->Smoother->ID==PASO_JACOBI)
511                Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
512            else if (amg->Smoother->ID==PASO_GS)    
513                Paso_Solver_solveGS(amg->Smoother->GS,x,b);
514        }        }
515         else {         else {
516         #ifdef MKL         #ifdef MKL
# Line 503  void Paso_Solver_solveAMG(Paso_Solver_AM Line 521  void Paso_Solver_solveAMG(Paso_Solver_AM
521               Paso_UMFPACK1(&ptr,amg->A,x,b,verbose);               Paso_UMFPACK1(&ptr,amg->A,x,b,verbose);
522               amg->solver=(void*) ptr;               amg->solver=(void*) ptr;
523            #else                  #else      
524               Paso_Solver_solveJacobi(amg->GS,x,b);             if(amg->Smoother->ID==PASO_JACOBI)
525                 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
526              else if (amg->Smoother->ID==PASO_GS)    
527                Paso_Solver_solveGS(amg->Smoother->GS,x,b);
528           #endif           #endif
529         #endif         #endif
530         }         }
# Line 514  void Paso_Solver_solveAMG(Paso_Solver_AM Line 535  void Paso_Solver_solveAMG(Paso_Solver_AM
535       } else {       } else {
536          /* presmoothing */          /* presmoothing */
537           time0=Paso_timer();           time0=Paso_timer();
538           Paso_Solver_solveJacobi(amg->GS,x,b);           if(amg->Smoother->ID==PASO_JACOBI)
539                Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
540            else if (amg->Smoother->ID==PASO_GS)    
541                Paso_Solver_solveGS(amg->Smoother->GS,x,b);
542                    
543          /***********/          /***********/
544          if (pre_sweeps>1) {          if (pre_sweeps>1) {
# Line 529  void Paso_Solver_solveAMG(Paso_Solver_AM Line 553  void Paso_Solver_solveAMG(Paso_Solver_AM
553              /* Compute the residual r=r-Ax*/              /* Compute the residual r=r-Ax*/
554             Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);             Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
555             /* Go round again*/             /* Go round again*/
556             Paso_Solver_solveJacobi(amg->GS,x,r);            
557               if(amg->Smoother->ID==PASO_JACOBI)
558                Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,r);
559               else if (amg->Smoother->ID==PASO_GS)    
560                Paso_Solver_solveGS(amg->Smoother->GS,x,r);
561              
562             pre_sweeps-=1;             pre_sweeps-=1;
563          }          }
564           /***********/           /***********/
# Line 571  void Paso_Solver_solveAMG(Paso_Solver_AM Line 600  void Paso_Solver_solveAMG(Paso_Solver_AM
600                
601        /*r=b-Ax */        /*r=b-Ax */
602        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
603        Paso_Solver_solveJacobi(amg->GS,x0,r);        if(amg->Smoother->ID==PASO_JACOBI)
604                Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
605          else if (amg->Smoother->ID==PASO_GS)    
606                Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
607                
608                
609        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
# Line 585  void Paso_Solver_solveAMG(Paso_Solver_AM Line 617  void Paso_Solver_solveAMG(Paso_Solver_AM
617             for (i=0;i<amg->n;++i) r[i]=b[i];             for (i=0;i<amg->n;++i) r[i]=b[i];
618                        
619             Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);             Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
620             Paso_Solver_solveJacobi(amg->GS,x0,r);            
621               if(amg->Smoother->ID==PASO_JACOBI)
622                  Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
623               else if (amg->Smoother->ID==PASO_GS)    
624                  Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
625                  
626             #pragma omp parallel for private(i) schedule(static)             #pragma omp parallel for private(i) schedule(static)
627              for (i=0;i<amg->n;++i)  {              for (i=0;i<amg->n;++i)  {
628               x[i]+=x0[i];               x[i]+=x0[i];

Legend:
Removed from v.2827  
changed lines
  Added in v.2828

  ViewVC Help
Powered by ViewVC 1.1.26