/[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 2474 by gross, Fri May 29 04:24:01 2009 UTC revision 2475 by artak, Wed Jun 17 01:48:46 2009 UTC
# Line 14  Line 14 
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: AMG preconditioner with reordering                 */  /* Paso: AMG preconditioner                                  */
18    
19  /**************************************************************/  /**************************************************************/
20    
# Line 24  Line 24 
24    
25  #include "Paso.h"  #include "Paso.h"
26  #include "Solver.h"  #include "Solver.h"
27    #include "Options.h"
28  #include "PasoUtil.h"  #include "PasoUtil.h"
29  #include "UMFPACK.h"  #include "UMFPACK.h"
30  #include "Pattern_coupling.h"  #include "Pattern_coupling.h"
# Line 71  to Line 72  to
72     then AMG is applied to S again until S becomes empty     then AMG is applied to S again until S becomes empty
73    
74  */  */
75  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level, double coarsening_threshold) {  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level, double coarsening_threshold, dim_t coarsening_method) {
76    Paso_Solver_AMG* out=NULL;    Paso_Solver_AMG* out=NULL;
77    dim_t n=A_p->numRows;    dim_t n=A_p->numRows;
78    dim_t n_block=A_p->row_block_size;    dim_t n_block=A_p->row_block_size;
# Line 82  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 83  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
83    Paso_SparseMatrix * schur=NULL;    Paso_SparseMatrix * schur=NULL;
84    Paso_SparseMatrix * schur_withFillIn=NULL;    Paso_SparseMatrix * schur_withFillIn=NULL;
85    double S=0;    double S=0;
86        
87    /*Make sure we have block sizes 1*/    /*Make sure we have block sizes 1*/
88    A_p->pattern->input_block_size=A_p->col_block_size;    A_p->pattern->input_block_size=A_p->col_block_size;
89    A_p->pattern->output_block_size=A_p->row_block_size;    A_p->pattern->output_block_size=A_p->row_block_size;
# Line 116  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 117  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
117       /* identify independend set of rows/columns */       /* identify independend set of rows/columns */
118       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
119       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
120       Paso_Pattern_coup(A_p,mis_marker,coarsening_threshold);      
121      /*Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);*/       /*Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);*/
122       /*Paso_Pattern_Aggregiation(A_p,mis_marker,coarsening_threshold);*/       /*fprintf(stderr,"\n Coarsening method == %d, %d \n",coarsening_method,PASO_RUGE_STUEBEN_COARSENING);*/
123       if (Paso_noError()) {       if (coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {
124          Paso_Pattern_coup(A_p,mis_marker,coarsening_threshold);
125         }
126         else if (coarsening_method == PASO_RUGE_STUEBEN_COARSENING) {
127          Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);
128         }
129         else if (coarsening_method == PASO_AGGREGATION_COARSENING) {
130          Paso_Pattern_Aggregiation(A_p,mis_marker,coarsening_threshold);
131         }
132         else {
133          /*Default coarseneing*/
134          Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);
135         }
136        
137        if (Paso_noError()) {
138          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
139          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
140          out->n=n;          out->n=n;
# Line 187  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 202  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
202                           if (Paso_noError()) {                           if (Paso_noError()) {
203                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
204                              /*find the pattern of the schur complement with fill in*/                              /*find the pattern of the schur complement with fill in*/
   
205                              schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);                              schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);
206                                                          
                             /*fprintf(stderr,"Sparsity of Schure: %dx%d LEN %d Percentage %f\n",schur_withFillIn->pattern->numOutput,schur_withFillIn->pattern->numInput,schur_withFillIn->len,schur_withFillIn->len/(1.*schur_withFillIn->pattern->numOutput*schur_withFillIn->pattern->numInput));*/  
                               
207                              /* copy values over*/                              /* copy values over*/
208                              #pragma omp parallel for private(i,iPtr,j,iPtr_s,index,where_p) schedule(static)                              #pragma omp parallel for private(i,iPtr,j,iPtr_s,index,where_p) schedule(static)
209                              for (i = 0; i < schur_withFillIn->numRows; ++i) {                              for (i = 0; i < schur_withFillIn->numRows; ++i) {
# Line 213  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 225  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
225                                                                    
226                              if (Paso_noError()) {                              if (Paso_noError()) {
227                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
228                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1,coarsening_threshold);                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1,coarsening_threshold,coarsening_method);
229                                  Paso_SparseMatrix_free(schur);                                  Paso_SparseMatrix_free(schur);
230                              }                              }
231                              /* allocate work arrays for AMG application */                              /* allocate work arrays for AMG application */
# Line 301  void Paso_Solver_solveAMG(Paso_Solver_AM Line 313  void Paso_Solver_solveAMG(Paso_Solver_AM
313       /*if (amg->n_F<=500) {*/       /*if (amg->n_F<=500) {*/
314        /*time0=Paso_timer();*/        /*time0=Paso_timer();*/
315                    
316          Paso_Solver_solveJacobi(amg->GS,x,b);          /*Paso_Solver_solveJacobi(amg->GS,x,b);*/
317                    
318          /* #pragma omp parallel for private(i) schedule(static)          /* #pragma omp parallel for private(i) schedule(static)
319          for (i=0;i<amg->n;++i) r[i]=b[i];          for (i=0;i<amg->n;++i) r[i]=b[i];
# Line 312  void Paso_Solver_solveAMG(Paso_Solver_AM Line 324  void Paso_Solver_solveAMG(Paso_Solver_AM
324           x[i]+=x0[i];           x[i]+=x0[i];
325          }          }
326          */          */
327          /*Paso_UMFPACK1(amg->A,x,b,0);*/          Paso_UMFPACK1(amg->A,x,b,0);
328                    
329         /*time0=Paso_timer()-time0;         /*
330         fprintf(stderr,"timing: DIRECT SOLVER: %e/\n",time0);*/         time0=Paso_timer()-time0;
331           fprintf(stderr,"timing: DIRECT SOLVER: %e/\n",time0);
332           */
333       } else {       } else {
334          /* presmoothing */          /* presmoothing */
335           Paso_Solver_solveJacobi(amg->GS,x,b);           Paso_Solver_solveJacobi(amg->GS,x,b);
336           /*          
337           #pragma omp parallel for private(i) schedule(static)           /* #pragma omp parallel for private(i) schedule(static)
338           for (i=0;i<amg->n;++i) r[i]=b[i];           for (i=0;i<amg->n;++i) r[i]=b[i];
339    
340            Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);            Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
341            Paso_Solver_solveGS(amg->GS,x0,r);            Paso_Solver_solveJacobi(amg->GS,x0,r);
342                        
343            #pragma omp parallel for private(i) schedule(static)            #pragma omp parallel for private(i) schedule(static)
344            for (i=0;i<amg->n;++i) {            for (i=0;i<amg->n;++i) {
345             x[i]+=x0[i];             x[i]+=x0[i];
346            }            }
347            */           */
348          /* end of presmoothing */          /* end of presmoothing */
349                    
350           #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
# Line 382  void Paso_Solver_solveAMG(Paso_Solver_AM Line 396  void Paso_Solver_solveAMG(Paso_Solver_AM
396       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
397       for (i=0;i<amg->n;++i) x[i]+=x0[i];       for (i=0;i<amg->n;++i) x[i]+=x0[i];
398            
399             /*
400       /*#pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
401       for (i=0;i<amg->n;++i) r[i]=b[i];       for (i=0;i<amg->n;++i) r[i]=b[i];
402        
403       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
404       Paso_Solver_solveJacobi(amg->GS,x0,r);       Paso_Solver_solveJacobi(amg->GS,x0,r);
405            
406       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
407       for (i=0;i<amg->n;++i) {       for (i=0;i<amg->n;++i) {
408        x[i]+=x0[i];        x[i]+=x0[i];
409       }*/       }
410             */
411       /*end of postsmoothing*/       /*end of postsmoothing*/
412            
413       }       }

Legend:
Removed from v.2474  
changed lines
  Added in v.2475

  ViewVC Help
Powered by ViewVC 1.1.26