/[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 2712 by artak, Wed Oct 7 00:22:43 2009 UTC revision 2726 by artak, Wed Oct 21 23:50:05 2009 UTC
# Line 116  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 116  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
116    Paso_SparseMatrix * schur_withFillIn=NULL;    Paso_SparseMatrix * schur_withFillIn=NULL;
117    double S=0;    double S=0;
118        
119      
120     /* char filename[8];
121      sprintf(filename,"AMGLevel%d",level);
122      
123     Paso_SparseMatrix_saveMM(A_p,filename);
124      */
125      
126    /*Make sure we have block sizes 1*/    /*Make sure we have block sizes 1*/
127    if (A_p->col_block_size>1) {    if (A_p->col_block_size>1) {
128       Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1.");       Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1.");
# Line 150  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 157  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
157       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/
158       out->level=level;       out->level=level;
159       out->n=n;       out->n=n;
160       out->n_F=0;       out->n_F=n+1;
161       out->n_block=n_block;       out->n_block=n_block;
162            
163       if (level==0 || n<=options->min_coarse_matrix_size) {       if (level==0 || n<=options->min_coarse_matrix_size) {
# Line 188  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 195  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
195             /*Default coarseneing*/             /*Default coarseneing*/
196              Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);              Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);
197              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/
198               /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/              /*Paso_Pattern_greedy_diag(A_p,mis_marker,options->coarsening_threshold);*/
199                /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
200          }          }
201                    
202          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
203          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
204            
205          out->n_F=Paso_Util_cumsum(n,counter);          out->n_F=Paso_Util_cumsum(n,counter);
206                    
207          if (out->n_F==n) {          if (out->n_F==0) {
208             out->coarsest_level=TRUE;             out->coarsest_level=TRUE;
209               level=0;
210             if (verbose) {             if (verbose) {
211                 printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");                 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
212                   printf("AMG coarsening does not eliminate any of the unknowns, switching to Jacobi preconditioner.\n");
213             }             }
214            }
215            else if (out->n_F==n) {
216              out->coarsest_level=TRUE;
217               level=0;
218               if (verbose) {
219                   /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
220                   printf("AMG coarsening eliminates all of the unknowns, switching to Jacobi preconditioner.\n");
221              
222                }
223          } else {          } else {
224            
225                if (Paso_noError()) {                if (Paso_noError()) {
# Line 390  void Paso_Solver_solveAMG(Paso_Solver_AM Line 410  void Paso_Solver_solveAMG(Paso_Solver_AM
410                
411        time0=Paso_timer();        time0=Paso_timer();
412        /*If all unknown are eliminated then Jacobi is the best preconditioner*/        /*If all unknown are eliminated then Jacobi is the best preconditioner*/
413        if (amg->n_F==amg->n) {        if (amg->n_F==0 || amg->n_F==amg->n) {
414          Paso_Solver_solveJacobi(amg->GS,x,b);          Paso_Solver_solveJacobi(amg->GS,x,b);
415        }        }
416         else {         else {
# Line 425  void Paso_Solver_solveAMG(Paso_Solver_AM Line 445  void Paso_Solver_solveAMG(Paso_Solver_AM
445           /*r=b-Ax*/           /*r=b-Ax*/
446           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
447                
448          /* b->[b_F,b_C]     */          /* r->[b_F,b_C]     */
449          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
450          for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]];          for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]];
451                    
# Line 449  void Paso_Solver_solveAMG(Paso_Solver_AM Line 469  void Paso_Solver_solveAMG(Paso_Solver_AM
469          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F);          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F);
470          /* x_F=invA_FF*b_F  */          /* x_F=invA_FF*b_F  */
471          Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);          Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);
472            
473          /* x<-[x_F,x_C]     */          /* x<-[x_F,x_C]     */
   
474          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
475          for (i=0;i<amg->n;++i) {          for (i=0;i<amg->n;++i) {
476              if (amg->mask_C[i]>-1) {              if (amg->mask_C[i]>-1) {

Legend:
Removed from v.2712  
changed lines
  Added in v.2726

  ViewVC Help
Powered by ViewVC 1.1.26