/[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 3193 by gross, Tue Sep 21 06:56:44 2010 UTC revision 3351 by gross, Tue Nov 16 02:39:40 2010 UTC
# Line 18  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /* Author: artak@uq.edu.au                                */  /* Author: artak@uq.edu.au, l.gross@uq.edu.au                                */
22    
23  /**************************************************************/  /**************************************************************/
24    
25    #define SHOW_TIMING FALSE
26    
27  #include "Paso.h"  #include "Paso.h"
28  #include "Preconditioner.h"  #include "Preconditioner.h"
29  #include "Options.h"  #include "Options.h"
30  #include "PasoUtil.h"  #include "PasoUtil.h"
31  #include "UMFPACK.h"  #include "UMFPACK.h"
32  #include "MKL.h"  #include "MKL.h"
 #include "Coarsening.h"  
33    
34  /**************************************************************/  /**************************************************************/
35    
# Line 37  Line 38 
38  void Paso_Preconditioner_LocalAMG_free(Paso_Preconditioner_LocalAMG * in) {  void Paso_Preconditioner_LocalAMG_free(Paso_Preconditioner_LocalAMG * in) {
39       if (in!=NULL) {       if (in!=NULL) {
40      Paso_Preconditioner_LocalSmoother_free(in->Smoother);      Paso_Preconditioner_LocalSmoother_free(in->Smoother);
41        Paso_SparseMatrix_free(in->P);
42        Paso_SparseMatrix_free(in->R);
43        Paso_SparseMatrix_free(in->A_C);
44        Paso_Preconditioner_LocalAMG_free(in->AMG_C);
45        MEMFREE(in->r);
46        MEMFREE(in->x_C);
47        MEMFREE(in->b_C);
48    
49      /*=========================*/      
50          Paso_SparseMatrix_free(in->A_FC);      MEMFREE(in);
         Paso_SparseMatrix_free(in->A_FF);  
         Paso_SparseMatrix_free(in->W_FC);  
         Paso_SparseMatrix_free(in->A_CF);  
         Paso_SparseMatrix_free(in->P);  
         Paso_SparseMatrix_free(in->R);  
         Paso_SparseMatrix_free(in->A);  
         if(in->coarsest_level==TRUE) {  
         #ifdef MKL  
           Paso_MKL_free1(in->AOffset1);  
           Paso_SparseMatrix_free(in->AOffset1);  
           Paso_SparseMatrix_free(in->AUnrolled);  
         #else  
           #ifdef UMFPACK  
           Paso_UMFPACK1_free((Paso_UMFPACK_Handler*)(in->solver));  
           Paso_SparseMatrix_free(in->AUnrolled);  
           #endif  
         #endif  
         }  
         MEMFREE(in->rows_in_F);  
         MEMFREE(in->rows_in_C);  
         MEMFREE(in->mask_F);  
         MEMFREE(in->mask_C);  
         MEMFREE(in->x_F);  
         MEMFREE(in->b_F);  
         MEMFREE(in->x_C);  
         MEMFREE(in->b_C);  
         in->solver=NULL;  
         Paso_Preconditioner_LocalAMG_free(in->AMG_of_Coarse);  
         MEMFREE(in->b_C);  
         MEMFREE(in);  
51       }       }
52  }  }
53    
54  /**************************************************************/  /*****************************************************************
   
 /*   constructs the block-block factorization of  
   
         [ A_FF A_FC ]  
    A_p=  
         [ A_CF A_FF ]  
   
 to  
   
   [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FF ]    
   [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ]  
   
   
    where S=A_FF-ACF*invA_FF*A_FC within the shape of S  
   
    then AMG is applied to S again until S becomes empty  
55    
56  */     constructs AMG
57      
58    ******************************************************************/
59  Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {  Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
60    
61    Paso_Preconditioner_LocalAMG* out=NULL;    Paso_Preconditioner_LocalAMG* out=NULL;
   /*  
    Paso_Pattern* temp1=NULL;  
   Paso_Pattern* temp2=NULL;  
   */  
62    bool_t verbose=options->verbose;    bool_t verbose=options->verbose;
   bool_t timing=0;  
63        
64      Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;
65    const dim_t n=A_p->numRows;    const dim_t n=A_p->numRows;
66    const dim_t n_block=A_p->row_block_size;    const dim_t n_block=A_p->row_block_size;
67        index_t* split_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree=NULL;
68        dim_t n_F=0, n_C=0, i;
   index_t* mis_marker=NULL;    
   index_t* counter=NULL;  
   /*index_t iPtr,*index, *where_p;*/  
   dim_t i,j;  
   Paso_SparseMatrix * A_c=NULL;  
69    double time0=0;    double time0=0;
70    Paso_SparseMatrix * Atemp=NULL;    const double theta = options->coarsening_threshold;
71    double sparsity=0;    const double tau = options->diagonal_dominance_threshold;
     
   /*  
   double *temp,*temp_1;  
   double S;  
   index_t iptr;  
   */  
72        
   /*char filename[8];*/  
73        
74    /*      /*
75    sprintf(filename,"AMGLevel%d",level);        is the input matrix A suitable for coarsening
76            
   Paso_SparseMatrix_saveMM(A_p,filename);  
77    */    */
78        if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {
79    /*Make sure we have block sizes 1*/       if (verbose) printf("Paso_Solver: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
80    /*if (A_p->col_block_size>1) {      level,  options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size  );  
      Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1.");  
81       return NULL;       return NULL;
82    }    }
83    if (A_p->row_block_size>1) {       /* Start Coarsening : */
      Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires row block size 1.");  
      return NULL;  
   }*/  
   out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);  
   
     
   /* identify independend set of rows/columns */  
   mis_marker=TMPMEMALLOC(n,index_t);  
   counter=TMPMEMALLOC(n,index_t);  
   if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) {  
       
       
      out->post_sweeps=options->post_sweeps;  
      out->pre_sweeps=options->pre_sweeps;  
       
      out->AMG_of_Coarse=NULL;  
      out->A_FF=NULL;  
      out->A_FC=NULL;  
      out->A_CF=NULL;  
      out->W_FC=NULL;  
      out->P=NULL;  
      out->R=NULL;  
      out->rows_in_F=NULL;  
      out->rows_in_C=NULL;  
      out->mask_F=NULL;  
      out->mask_C=NULL;  
      out->x_F=NULL;  
      out->b_F=NULL;  
      out->x_C=NULL;  
      out->b_C=NULL;  
      out->A=Paso_SparseMatrix_getReference(A_p);  
      out->AUnrolled=NULL;  
      out->AOffset1=NULL;  
      out->solver=NULL;  
      out->Smoother=NULL;  
      out->level=level;  
      out->n=n;  
      out->n_F=n+1;  
      out->n_block=n_block;  
84    
85             split_marker=TMPMEMALLOC(n,index_t);
86       sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols);       counter=TMPMEMALLOC(n,index_t);
87             degree=TMPMEMALLOC(n, dim_t);
88       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);       S=TMPMEMALLOC(A_p->pattern->len, index_t);
89             if ( !( Esys_checkPtr(split_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree) || Esys_checkPtr(S) ) ) {
90             /*
91       if(sparsity>0.05) {            set splitting of unknows:
92        level=0;        
93       }           */
94             time0=Esys_timer();
95                 if (n_block>1) {
96       if (level==0 || n<=options->min_coarse_matrix_size) {             Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree, S, theta,tau);
97           out->coarsest_level=TRUE;       } else {
98           #ifdef MKL             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree, S, theta,tau);
99                    out->AUnrolled=Paso_SparseMatrix_unroll(A_p);       }
100                    out->AOffset1=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, out->AUnrolled->pattern,1,1, FALSE);       Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree, S, split_marker);
101                    #pragma omp parallel for private(i) schedule(static)       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
102                    for (i=0;i<out->A->len;++i) {      
103                         out->AOffset1->val[i]=out->AUnrolled->val[i];       if (Esys_noError() ) {
104                    }          #pragma omp parallel for private(i) schedule(static)
105           #else          for (i = 0; i < n; ++i) split_marker[i]= (split_marker[i] == PASO_AMG_IN_F);
106              #ifdef UMFPACK      
107                  out->AUnrolled=Paso_SparseMatrix_unroll(A_p);          /*
108                  /*Paso_SparseMatrix_saveMM(out->AUnrolled,"AUnroll.mat");             count number of unkowns to be eliminated:
109                  Paso_SparseMatrix_saveMM(A_p,"Aorg.mat");          */
110                  */          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);
111              #else          n_C=n-n_F;
112                out->Smoother=Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);          if (verbose) printf("Paso_Solver: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
113              #endif      
114           #endif          if ( n_F == 0 ) {  /*  is a nasty case. a direct solver should be used, return NULL */
115                       out = NULL;
116       } else {          } else {
117           out->coarsest_level=FALSE;             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);
118       out->Smoother=Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);             if (! Esys_checkPtr(out)) {
119              out->level = level;
120           /* Start Coarsening : */            out->n = n;
121           time0=Paso_timer();            out->n_F = n_F;
122       Paso_Coarsening_Local(mis_marker, A_p, options->coarsening_threshold, options->coarsening_method);            out->n_block = n_block;
123           time0=Paso_timer()-time0;            out->A_C = NULL;
124       if (timing) fprintf(stdout,"Level %d: timing: Coarsening: %e\n",level,time0);            out->P = NULL;  
125              out->R = NULL;          
126          #pragma omp parallel for private(i) schedule(static)            out->post_sweeps = options->post_sweeps;
127          for (i = 0; i < n; ++i) {            out->pre_sweeps  = options->pre_sweeps;
128         mis_marker[i]=(mis_marker[i]== PASO_COARSENING_IN_F);            out->r = NULL;
129         counter[i]=mis_marker[i];            out->x_C = NULL;
130          }            out->b_C = NULL;
131          out->n_F=Paso_Util_cumsum(n,counter);            out->AMG_C = NULL;
132                      out->Smoother=NULL;
133          if (out->n_F==0) {             }
134             out->coarsest_level=TRUE;             mask_C=TMPMEMALLOC(n,index_t);
135             level=0;             rows_in_F=TMPMEMALLOC(n_F,index_t);
136             if (verbose) {             Esys_checkPtr(mask_C);
137                 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/             Esys_checkPtr(rows_in_F);
138                 printf("AMG coarsening does not eliminate any of the unknowns, switching to Jacobi preconditioner.\n");             if ( Esys_noError() ) {
139             }  
140          }            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
141          else if (out->n_F==n) {        
142            out->coarsest_level=TRUE;            if ( n_F < n ) { /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
143             level=0;    
144             if (verbose) {              /* allocate helpers :*/
145                 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/              out->x_C=MEMALLOC(n_block*n_C,double);
146                 printf("AMG coarsening eliminates all of the unknowns, switching to Jacobi preconditioner.\n");              out->b_C=MEMALLOC(n_block*n_C,double);
147                          out->r=MEMALLOC(n_block*n,double);
148              }              
149          } else {              Esys_checkPtr(out->r);
150                    Esys_checkPtr(out->x_C);
151                if (Paso_noError()) {              Esys_checkPtr(out->b_C);
152                                
153                   /*#pragma omp parallel for private(i) schedule(static)              if ( Esys_noError() ) {
154                   for (i = 0; i < n; ++i) counter[i]=mis_marker[i];                 /* creates index for F:*/
155                   out->n_F=Paso_Util_cumsum(n,counter);                 #pragma omp parallel private(i)
156                   */                 {
157                                      #pragma omp for schedule(static)
158                   out->mask_F=MEMALLOC(n,index_t);                    for (i = 0; i < n; ++i) {
159                   out->rows_in_F=MEMALLOC(out->n_F,index_t);                   if  (split_marker[i]) rows_in_F[counter[i]]=i;
160                   if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->rows_in_F) ) ) {                    }
161                      /* creates an index for F from mask */                 }
162                      #pragma omp parallel for private(i) schedule(static)                 /*  create mask of C nodes with value >-1 gives new id */
163                      for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;                 i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);
164                      #pragma omp parallel for private(i) schedule(static)  
165                      for (i = 0; i < n; ++i) {                 #pragma omp parallel for private(i) schedule(static)
166                 if  (mis_marker[i]) {                 for (i = 0; i < n; ++i) {
167                                out->rows_in_F[counter[i]]=i;                    if  (split_marker[i]) {
168                                out->mask_F[i]=counter[i];                   mask_C[i]=-1;
169                         } else {                    } else {
170                                out->mask_F[i]=-1;                   mask_C[i]=counter[i];;
171                         }                    }
172                      }                 }
173                                       /*
174                   }                    get Restriction :  
175                }                 */                  
176                                 time0=Esys_timer();
177                /* if(level==1) {                 out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);
178                     printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F);                 if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);
179                     for (i = 0; i < n; ++i) {              }
180                      printf("##%d %d\n",i,!mis_marker[i]);              /*      
181                     }                 construct Prolongation operator as transposed of restriction operator:
182                  }              */
183                */              if ( Esys_noError()) {
184                                 time0=Esys_timer();
185                /*check whether coarsening process actually makes sense to continue.                 out->R=Paso_SparseMatrix_getTranspose(out->P);
186                if coarse matrix at least smaller by 30% then continue, otherwise we stop.*/                 if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
187                if ((out->n_F*100/n)<30) {              }      
188                      level=1;              /*
189                  }              construct coarse level matrix:
190                            */
191                if ( Paso_noError()) {              if ( Esys_noError()) {
192                      out->n_C=n-out->n_F;                 time0=Esys_timer();
193                      out->rows_in_C=MEMALLOC(out->n_C,index_t);                 Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
194                      out->mask_C=MEMALLOC(n,index_t);                 A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
195                      if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {                 Paso_SparseMatrix_free(Atemp);
196                           /* creates an index for C from mask */                 if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);            
197                           #pragma omp parallel for private(i) schedule(static)              }
198                           for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];  
199                           Paso_Util_cumsum(n,counter);              
200                           #pragma omp parallel for private(i) schedule(static)              /*
201                           for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;                 constructe courser level:
202                           #pragma omp parallel for private(i) schedule(static)                
203                           for (i = 0; i < n; ++i) {              */
204                  if  (! mis_marker[i]) {              if ( Esys_noError()) {
205                                        out->rows_in_C[counter[i]]=i;                 out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);
206                                        out->mask_C[i]=counter[i];              }
207                                     } else {              if ( Esys_noError()) {
208                                        out->mask_C[i]=-1;                 if ( out->AMG_C == NULL ) {
209                                     }                    out->reordering = options->reordering;
210                           }                    out->refinements = options->coarse_matrix_refinements;
211                      }                    /* no coarse level matrix has been constructed. use direct solver */
212                }                    #ifdef MKL
213                if ( Paso_noError()) {                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
214                        /* get A_FF block: */                      Paso_SparseMatrix_free(A_C);
215                        /*                      out->A_C->solver_package = PASO_MKL;
216                         out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F);                      if (verbose) printf("Paso_Solver: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C);
217                        out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);                    #else
218                        out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);                      #ifdef UMFPACK
219                        */                         out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
220                                                 Paso_SparseMatrix_free(A_C);
221                        /*Compute W_FC*/                         out->A_C->solver_package = PASO_UMFPACK;
222                        /*initialy W_FC=A_FC*/                         if (verbose) printf("Paso_Solver: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C);
223                        out->W_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);                      #else
224                                                 out->A_C=A_C;
225                        /*sprintf(filename,"W_FCbefore_%d",level);                         out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);
226                        Paso_SparseMatrix_saveMM(out->W_FC,filename);                         out->A_C->solver_package = PASO_SMOOTHER;
227                        */                         if (verbose) printf("Paso_Solver: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C);
228                        /* for (i = 0; i < n; ++i) {                      #endif
229                                    printf("##mis_marker[%d]=%d\n",i,mis_marker[i]);                    #endif
230                         }                 } else {
231                        */                                      /* finally we set some helpers for the solver step */
232                        time0=Paso_timer();                    out->A_C=A_C;
233                        Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker);                 }
234                        time0=Paso_timer()-time0;              }        
235                        if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0);            }
236                       }
237                        /*             TMPMEMFREE(mask_C);
238                        sprintf(filename,"W_FCafter_%d",level);             TMPMEMFREE(rows_in_F);
239                        Paso_SparseMatrix_saveMM(out->W_FC,filename);          }
240                        */       }
                         
                       /* get Prolongation and Restriction */  
                       time0=Paso_timer();  
                       out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker);  
                       time0=Paso_timer()-time0;  
                       if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0);  
                       /*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/  
                         
                       /*  
                       sprintf(filename,"P_%d",level);  
                       Paso_SparseMatrix_saveMM(out->P,filename);  
                       */  
                         
                       time0=Paso_timer();  
                       out->R=Paso_SparseMatrix_getRestriction(out->P);  
                       time0=Paso_timer()-time0;  
                       if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0);  
                       /*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/  
                         
                       /*  
                       sprintf(filename,"R_%d",level);  
                       Paso_SparseMatrix_saveMM(out->R,filename);  
                       */  
                       
               }  
               if ( Paso_noError()) {  
                       
                     time0=Paso_timer();  
                       
                     Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);  
                       
                     A_c=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);  
                       
                     /*A_c=Paso_SparseMatrix_loadMM_toCSR("A_C1.mtx");*/  
                       
                     Paso_SparseMatrix_free(Atemp);  
                       
                     /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/  
                     time0=Paso_timer()-time0;  
                     if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0);  
                       
                                           
                     /*Paso_Solver_getCoarseMatrix(A_c, A_p,out->R,out->P);*/  
                     /*                      
                     sprintf(filename,"A_C_%d",level);  
                     Paso_SparseMatrix_saveMM(A_c,filename);  
                     */  
                                           
             out->AMG_of_Coarse=Paso_Preconditioner_LocalAMG_alloc(A_c,level-1,options);  
               }  
   
               /* allocate work arrays for AMG application */  
               if (Paso_noError()) {  
                          /*  
                           out->x_F=MEMALLOC(n_block*out->n_F,double);  
                          out->b_F=MEMALLOC(n_block*out->n_F,double);  
                          */  
                          out->x_C=MEMALLOC(n_block*out->n_C,double);  
                          out->b_C=MEMALLOC(n_block*out->n_C,double);  
         
                          /*if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {*/  
                          if ( ! ( Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {  
                               
                              /*  
                               #pragma omp parallel for private(i) schedule(static)  
                              for (i = 0; i < out->n_F; ++i) {  
                                          out->x_F[i]=0.;  
                                          out->b_F[i]=0.;  
                               }  
                              */  
                               
                               #pragma omp parallel for private(i,j) schedule(static)  
                               for (i = 0; i < out->n_C; ++i) {  
                                      for(j=0;j<n_block;++j) {  
                                         out->x_C[i*n_block+j]=0.;  
                                         out->b_C[i*n_block+j]=0.;  
                                      }  
                               }  
                          }  
               }  
             Paso_SparseMatrix_free(A_c);  
          }  
      }  
241    }    }
   TMPMEMFREE(mis_marker);  
242    TMPMEMFREE(counter);    TMPMEMFREE(counter);
243      TMPMEMFREE(split_marker);
244      TMPMEMFREE(degree);
245      TMPMEMFREE(S);
246    
247    if (Paso_noError()) {    if (Esys_noError()) {
       if (verbose && level>0 && !out->coarsest_level) {  
          printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);  
      }  
248       return out;       return out;
249    } else  {    } else  {
250       Paso_Preconditioner_LocalAMG_free(out);       Paso_Preconditioner_LocalAMG_free(out);
# Line 428  Paso_Preconditioner_LocalAMG* Paso_Preco Line 252  Paso_Preconditioner_LocalAMG* Paso_Preco
252    }    }
253  }  }
254    
 /**************************************************************/  
255    
256  /* apply AMG precondition b-> x                                void Paso_Preconditioner_LocalAMG_solve(Paso_SparseMatrix* A, Paso_Preconditioner_LocalAMG * amg, double * x, double * b) {
257         const dim_t n = amg->n * amg->n_block;
258         double time0=0;
259         const dim_t post_sweeps=amg->post_sweeps;
260         const dim_t pre_sweeps=amg->pre_sweeps;
261    
262       in fact it solves       /* presmoothing */
263         time0=Esys_timer();
264         Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
265         time0=Esys_timer()-time0;
266         if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);
267         /* end of presmoothing */
268        
269         if (amg->n_F < amg->n) { /* is there work on the coarse level? */
270             time0=Esys_timer();
271         Paso_Copy(n, amg->r, b);                            /*  r <- b */
272         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
273         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */
274             time0=Esys_timer()-time0;
275         /* coarse level solve */
276         if ( amg->AMG_C == NULL) {
277            time0=Esys_timer();
278            /*  A_C is the coarsest level */
279            switch (amg->A_C->solver_package) {
280               case (PASO_MKL):
281              Paso_MKL(amg->A_C, amg->x_C,amg->b_C, amg->reordering, amg->refinements, SHOW_TIMING);
282              break;
283               case (PASO_UMFPACK):
284              Paso_UMFPACK(amg->A_C, amg->x_C,amg->b_C, amg->refinements, SHOW_TIMING);
285              break;
286               case (PASO_SMOOTHER):
287              Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->Smoother,amg->x_C,amg->b_C,pre_sweeps, FALSE);
288              break;
289            }
290            if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
291         } else {
292            Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */
293         }  
294         time0=time0+Esys_timer();
295         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */    
296        
297             /*postsmoothing*/
298          
299            /*solve Ax=b with initial guess x */
300            time0=Esys_timer();
301            Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
302            time0=Esys_timer()-time0;
303            if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
304            /*end of postsmoothing*/
305        
306         }
307         return;
308    }
309    
310    [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FC ]  [ x_F ]  = [b_F]  /* theta = threshold for strong connections */
311    [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ]  [ x_C ]  = [b_C]  /* tau = threshold for diagonal dominance */
312    
313   in the form  /*S_i={j \in N_i; i strongly coupled to j}
314    
315     b->[b_F,b_C]  in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
316     x_F=invA_FF*b_F  */
    b_C=b_C-A_CF*x_F  
    x_C=AMG(b_C)  
    b_F=b_F-A_FC*x_C  
    x_F=invA_FF*b_F  
    x<-[x_F,x_C]  
317    
318   should be called within a parallel region                                                void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,
319   barrier synconization should be performed to make sure that the input vector available                        dim_t *degree, index_t *S,
320                          const double theta, const double tau)
321    {
322       const dim_t n=A->numRows;
323       index_t iptr, i,j;
324       dim_t kdeg;
325       double max_offdiagonal, threshold, sum_row, main_row, fnorm;
326    
 */  
327    
328  void Paso_Preconditioner_LocalAMG_solve(Paso_Preconditioner_LocalAMG * amg, double * x, double * b) {        #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)
329       dim_t i,j;        for (i=0;i<n;++i) {
      double time0=0;  
      double *r=NULL, *x0=NULL;  
      bool_t timing=0;  
       
      dim_t post_sweeps=amg->post_sweeps;  
      dim_t pre_sweeps=amg->pre_sweeps;  
       
      #ifdef UMFPACK  
           Paso_UMFPACK_Handler * ptr=NULL;  
      #endif  
330                        
331       r=MEMALLOC(amg->n*amg->n_block,double);       max_offdiagonal = 0.;
332       x0=MEMALLOC(amg->n*amg->n_block,double);       sum_row=0;
333             main_row=0;
334       if (amg->coarsest_level) {       #pragma ivdep
335               for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
336        time0=Paso_timer();          j=A->pattern->index[iptr];
337        /*If all unknown are eliminated then Jacobi is the best preconditioner*/          fnorm=ABS(A->val[iptr]);
338            
339        if (amg->n_F==0 || amg->n_F==amg->n) {          if( j != i) {
340         Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother,x,b,1, FALSE);             max_offdiagonal = MAX(max_offdiagonal,fnorm);
341               sum_row+=fnorm;
342            } else {
343               main_row=fnorm;
344            }
345         }
346         threshold = theta*max_offdiagonal;
347         kdeg=0;
348         if (tau*main_row < sum_row) { /* no diagonal domainance */
349            #pragma ivdep
350            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
351               j=A->pattern->index[iptr];
352               if(ABS(A->val[iptr])>threshold && i!=j) {
353              S[A->pattern->ptr[i]+kdeg] = j;
354              kdeg++;
355               }
356            }
357         }
358         degree[i]=kdeg;
359        }        }
360         else {  
361         #ifdef MKL  }
362            Paso_MKL1(amg->AOffset1,x,b,timing);  
363         #else  /* theta = threshold for strong connections */
364            #ifdef UMFPACK  /* tau = threshold for diagonal dominance */
365               ptr=(Paso_UMFPACK_Handler *)(amg->solver);  /*S_i={j \in N_i; i strongly coupled to j}
366               Paso_UMFPACK1(&ptr,amg->AUnrolled,x,b,timing);  
367               amg->solver=(void*) ptr;  in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
368            #else        */
369            Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother,x,b,1, FALSE);  void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,
370           #endif                              dim_t *degree, index_t *S,
371         #endif                              const double theta, const double tau)
372         }  
373    {
374       const dim_t n_block=A->row_block_size;
375       const dim_t n=A->numRows;
376       index_t iptr, i,j, bi;
377       dim_t kdeg, max_deg;
378       register double max_offdiagonal, threshold, fnorm, sum_row, main_row;
379       double *rtmp;
380      
381      
382          #pragma omp parallel private(i,iptr,max_offdiagonal, kdeg, threshold,j, max_deg, fnorm, sum_row, main_row, rtmp)
383          {
384         max_deg=0;
385         #pragma omp for schedule(static)
386         for (i=0;i<n;++i) max_deg=MAX(max_deg, A->pattern->ptr[i+1]-A->pattern->ptr[i]);
387          
388         rtmp=TMPMEMALLOC(max_deg, double);
389                
390         time0=Paso_timer()-time0;       #pragma omp for schedule(static)
391         if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0);       for (i=0;i<n;++i) {
         
      } else {  
         /* presmoothing */  
          time0=Paso_timer();  
      Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, pre_sweeps, FALSE);  
          time0=Paso_timer()-time0;  
          if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0);  
392            
393           /* end of presmoothing */          max_offdiagonal = 0.;
394                    sum_row=0;
395           time0=Paso_timer();          main_row=0;
396            #pragma omp parallel for private(i,j) schedule(static)          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
397             for (i=0;i<amg->n;++i) {             j=A->pattern->index[iptr];
398              for (j=0;j<amg->n_block;++j) {             fnorm=0;
399                r[i*amg->n_block+j]=b[i*amg->n_block+j];               #pragma ivdep
400              }             for(bi=0;bi<n_block*n_block;++bi) fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
401             }             fnorm=sqrt(fnorm);
402                       rtmp[iptr-A->pattern->ptr[i]]=fnorm;
403           /*r=b-Ax*/             if( j != i) {
404           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);            max_offdiagonal = MAX(max_offdiagonal,fnorm);
405                      sum_row+=fnorm;
406          /* b_c = R*r  */             } else {
407           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C);            main_row=fnorm;
408                       }
409          time0=Paso_timer()-time0;          }
410          if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0);          threshold = theta*max_offdiagonal;
           
         /* x_C=AMG(b_C)     */  
     Paso_Preconditioner_LocalAMG_solve(amg->AMG_of_Coarse,amg->x_C,amg->b_C);  
           
         time0=Paso_timer();  
           
         /* x_0 = P*x_c */  
         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0);  
           
         /* x=x+x0 */  
          #pragma omp parallel for private(i,j) schedule(static)  
            for (i=0;i<amg->n;++i) {  
             for (j=0;j<amg->n_block;++j) {  
               x[i*amg->n_block+j]+=x0[i*amg->n_block+j];    
             }  
            }  
           
       /*postsmoothing*/  
411                
412        /*solve Ax=b with initial guess x */          kdeg=0;
413        time0=Paso_timer();          if (tau*main_row < sum_row) { /* no diagonal domainance */
414        Paso_Preconditioner_LocalSmoother_solve(amg->A, amg->Smoother, x, b, post_sweeps, TRUE);             #pragma ivdep
415        time0=Paso_timer()-time0;             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
416        if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);            j=A->pattern->index[iptr];
417              if(rtmp[iptr-A->pattern->ptr[i]] > threshold && i!=j) {
418                 S[A->pattern->ptr[i]+kdeg] = j;
419                 kdeg++;
420              }
421               }
422            }
423            degree[i]=kdeg;
424         }      
425         TMPMEMFREE(rtmp);
426          } /* end of parallel region */
427    
428        /*end of postsmoothing*/  }  
       
      }  
      MEMFREE(r);  
      MEMFREE(x0);  
429    
430       return;  /* the runge stueben coarsening algorithm: */
431    void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset,
432                             const dim_t* degree, const index_t* S,
433                             index_t*split_marker)
434    {
435       const bool_t usePanel=FALSE;
436      
437       index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;
438       dim_t i,k, p, q, *degreeT=NULL, len_panel, len_panel_new;
439       register index_t j, itmp;
440      
441       if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
442      
443       lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);
444       degreeT=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degreeT);
445       ST=TMPMEMALLOC(offset[n], index_t);  Esys_checkPtr(ST);
446       if (usePanel) {
447          notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);
448          panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);
449       }
450      
451      
452      
453       if (Esys_noError() ) {
454          /* initialize  split_marker and split_marker :*/
455          /* those unknows which are not influenced go into F, the rest is available for F or C */
456          #pragma omp parallel for private(i) schedule(static)
457          for (i=0;i<n;++i) {
458         degreeT[i]=0;
459         if (degree[i]>0) {
460            lambda[i]=0;
461            split_marker[i]=PASO_AMG_UNDECIDED;
462         } else {
463            split_marker[i]=PASO_AMG_IN_F;
464            lambda[i]=-1;
465         }
466          }
467          /* create transpose :*/
468          for (i=0;i<n;++i) {
469            for (p=0; p<degree[i]; ++p) {
470               j=S[offset[i]+p];
471               ST[offset[j]+degreeT[j]]=i;
472               degreeT[j]++;
473            }
474          }
475          /* lambda[i] = |undecided k in ST[i]| + 2 * |F-unknown in ST[i]| */
476          #pragma omp parallel for private(i, j, itmp) schedule(static)
477          for (i=0;i<n;++i) {
478         if (split_marker[i]==PASO_AMG_UNDECIDED) {
479            itmp=lambda[i];
480            for (p=0; p<degreeT[i]; ++p) {
481               j=ST[offset[i]+p];
482               if (split_marker[j]==PASO_AMG_UNDECIDED) {
483              itmp++;
484               } else {  /* at this point there are no C points */
485              itmp+=2;
486               }
487            }
488            lambda[i]=itmp;
489         }
490          }
491          if (usePanel) {
492         #pragma omp parallel for private(i) schedule(static)
493         for (i=0;i<n;++i) notInPanel[i]=TRUE;
494          }
495          /* start search :*/
496          i=Paso_Util_arg_max(n,lambda);
497          while (lambda[i]>-1) { /* is there any undecided unknown? */
498    
499         if (usePanel) {
500            len_panel=0;
501            do {
502               /* the unknown i is moved to C */
503               split_marker[i]=PASO_AMG_IN_C;
504               lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
505              
506               /* all undecided unknown strongly coupled to i are moved to F */
507               for (p=0; p<degreeT[i]; ++p) {
508              j=ST[offset[i]+p];
509              
510              if (split_marker[j]==PASO_AMG_UNDECIDED) {
511                
512                 split_marker[j]=PASO_AMG_IN_F;
513                 lambda[j]=-1;
514                
515                 for (q=0; q<degreeT[j]; ++q) {
516                k=ST[offset[j]+q];
517                if (split_marker[k]==PASO_AMG_UNDECIDED) {
518                   lambda[k]++;
519                   if (notInPanel[k]) {
520                      notInPanel[k]=FALSE;
521                      panel[len_panel]=k;
522                      len_panel++;
523                   }
524    
525                }    /* the unknown i is moved to C */
526                split_marker[i]=PASO_AMG_IN_C;
527                lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
528                 }
529                
530              }
531               }
532               for (p=0; p<degree[i]; ++p) {
533              j=S[offset[i]+p];
534              if (split_marker[j]==PASO_AMG_UNDECIDED) {
535                 lambda[j]--;
536                 if (notInPanel[j]) {
537                notInPanel[j]=FALSE;
538                panel[len_panel]=j;
539                len_panel++;
540                 }
541              }
542               }
543    
544               /* consolidate panel */
545               /* remove lambda[q]=-1 */
546               lambda_max=-1;
547               i=-1;
548               len_panel_new=0;
549               for (q=0; q<len_panel; q++) {
550                 k=panel[q];
551                 lambda_k=lambda[k];
552                 if (split_marker[k]==PASO_AMG_UNDECIDED) {
553                panel[len_panel_new]=k;
554                len_panel_new++;
555    
556                if (lambda_max == lambda_k) {
557                   if (k<i) i=k;
558                } else if (lambda_max < lambda_k) {
559                   lambda_max =lambda_k;
560                   i=k;
561                }
562                 }
563               }
564               len_panel=len_panel_new;
565            } while (len_panel>0);    
566         } else {
567            /* the unknown i is moved to C */
568            split_marker[i]=PASO_AMG_IN_C;
569            lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
570            
571            /* all undecided unknown strongly coupled to i are moved to F */
572            for (p=0; p<degreeT[i]; ++p) {
573               j=ST[offset[i]+p];
574               if (split_marker[j]==PASO_AMG_UNDECIDED) {
575            
576              split_marker[j]=PASO_AMG_IN_F;
577              lambda[j]=-1;
578            
579              for (q=0; q<degreeT[j]; ++q) {
580                 k=ST[offset[j]+q];
581                 if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;
582              }
583    
584               }
585            }
586            for (p=0; p<degree[i]; ++p) {
587               j=S[offset[i]+p];
588               if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;
589            }
590            
591         }
592         i=Paso_Util_arg_max(n,lambda);
593          }
594       }
595       TMPMEMFREE(lambda);
596       TMPMEMFREE(ST);
597       TMPMEMFREE(degreeT);
598       TMPMEMFREE(panel);
599       TMPMEMFREE(notInPanel);
600  }  }

Legend:
Removed from v.3193  
changed lines
  Added in v.3351

  ViewVC Help
Powered by ViewVC 1.1.26