/[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 2802 by artak, Thu Dec 3 01:51:55 2009 UTC revision 2803 by artak, Thu Dec 3 05:35:20 2009 UTC
# Line 170  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 170  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
170       out->n=n;       out->n=n;
171       out->n_F=n+1;       out->n_F=n+1;
172       out->n_block=n_block;       out->n_block=n_block;
173         out->post_sweeps=options->post_sweeps;
174         out->pre_sweeps=options->pre_sweeps;
175            
176       sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols);       sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols);
177            
178       if (verbose) fprintf(stdout,"Stats: Sparsity of the Coarse Matrix with %d non-zeros (%d,%d) in level %d is %.6f\n",A_p->len,A_p->numRows,A_p->numCols,level,sparsity);       if (verbose) fprintf(stdout,"Stats: Sparsity of the Coarse Matrix with %d non-zeros (%d,%d) in level %d is %.6f\n",A_p->len,A_p->numRows,A_p->numCols,level,sparsity);
179            
180            
181       if(sparsity>0.01) {       /*if(sparsity>0.01) {
182        level=0;        level=0;
183       }       }
184             */
185                    
186       if (level==0 || n<=options->min_coarse_matrix_size) {       if (level==0 || n<=options->min_coarse_matrix_size) {
187           out->coarsest_level=TRUE;           out->coarsest_level=TRUE;
# Line 220  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 222  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
222          }          }
223          else {          else {
224             /*Default coarseneing*/             /*Default coarseneing*/
225              /*Paso_Pattern_RS_MI_Aggressive(A_p,mis_marker,options->coarsening_threshold);*/              Paso_Pattern_RS_MI_Aggressive(A_p,mis_marker,options->coarsening_threshold);
226              Paso_Pattern_RS_MI(A_p,mis_marker,options->coarsening_threshold);              /*Paso_Pattern_RS_MI(A_p,mis_marker,options->coarsening_threshold);*/
227              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/
228              /*Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);*/
229              /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
# Line 468  void Paso_Solver_solveAMG(Paso_Solver_AM Line 470  void Paso_Solver_solveAMG(Paso_Solver_AM
470       double *r=NULL, *x0=NULL;       double *r=NULL, *x0=NULL;
471       bool_t verbose=0;       bool_t verbose=0;
472            
473         dim_t post_sweeps=amg->post_sweeps;
474         dim_t pre_sweeps=amg->pre_sweeps;
475        
476       #ifdef UMFPACK       #ifdef UMFPACK
477            Paso_UMFPACK_Handler * ptr=NULL;            Paso_UMFPACK_Handler * ptr=NULL;
478       #endif       #endif
# Line 506  void Paso_Solver_solveAMG(Paso_Solver_AM Line 511  void Paso_Solver_solveAMG(Paso_Solver_AM
511          /* presmoothing */          /* presmoothing */
512           time0=Paso_timer();           time0=Paso_timer();
513           Paso_Solver_solveJacobi(amg->GS,x,b);           Paso_Solver_solveJacobi(amg->GS,x,b);
514            
515            /***********/
516            #pragma omp parallel for private(i) schedule(static)
517            for (i=0;i<amg->n;++i) r[i]=b[i];
518      
519            while(pre_sweeps>1) {
520               #pragma omp parallel for private(i) schedule(static)
521               for (i=0;i<amg->n;++i) r[i]+=b[i];
522              
523                /* Compute the residual b=b-Ax*/
524               Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
525               /* Go round again*/
526               Paso_Solver_solveJacobi(amg->GS,x,r);
527               pre_sweeps-=1;
528            }
529             /***********/
530            
531           time0=Paso_timer()-time0;           time0=Paso_timer()-time0;
532           if (verbose) fprintf(stdout,"timing: Presmooting: %e\n",time0);           if (verbose) fprintf(stdout,"timing: Presmooting: %e\n",time0);
533          /* end of presmoothing */           /* end of presmoothing */
           
534                    
535           time0=Paso_timer();           time0=Paso_timer();
536           #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
# Line 518  void Paso_Solver_solveAMG(Paso_Solver_AM Line 539  void Paso_Solver_solveAMG(Paso_Solver_AM
539           /*r=b-Ax*/           /*r=b-Ax*/
540           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
541                    
           
542          /* b_c <- R*r  */          /* b_c <- R*r  */
543           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);
544                    
# Line 538  void Paso_Solver_solveAMG(Paso_Solver_AM Line 558  void Paso_Solver_solveAMG(Paso_Solver_AM
558          for (i=0;i<amg->n;++i) x[i]+=x0[i];          for (i=0;i<amg->n;++i) x[i]+=x0[i];
559                    
560        /*postsmoothing*/        /*postsmoothing*/
561    
562        time0=Paso_timer();        time0=Paso_timer();
563        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
564        for (i=0;i<amg->n;++i) r[i]=b[i];        for (i=0;i<amg->n;++i) r[i]=b[i];
# Line 551  void Paso_Solver_solveAMG(Paso_Solver_AM Line 572  void Paso_Solver_solveAMG(Paso_Solver_AM
572        for (i=0;i<amg->n;++i)  {        for (i=0;i<amg->n;++i)  {
573         x[i]+=x0[i];         x[i]+=x0[i];
574        }        }
575            /***************/
576            while(post_sweeps>1) {
577              
578               #pragma omp parallel for private(i) schedule(static)
579               for (i=0;i<amg->n;++i) r[i]=b[i];
580              
581               Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
582               Paso_Solver_solveJacobi(amg->GS,x0,r);
583               #pragma omp parallel for private(i) schedule(static)
584                for (i=0;i<amg->n;++i)  {
585                 x[i]+=x0[i];
586                }
587               post_sweeps-=1;
588            }
589            /**************/
590                
591        time0=Paso_timer()-time0;        time0=Paso_timer()-time0;
592        if (verbose) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);        if (verbose) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);

Legend:
Removed from v.2802  
changed lines
  Added in v.2803

  ViewVC Help
Powered by ViewVC 1.1.26