/[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 3287 by gross, Tue Oct 19 00:42:22 2010 UTC revision 3303 by gross, Mon Oct 25 04:33:31 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    
# Line 61  Paso_Preconditioner_LocalAMG* Paso_Preco Line 61  Paso_Preconditioner_LocalAMG* Paso_Preco
61    Paso_Preconditioner_LocalAMG* out=NULL;    Paso_Preconditioner_LocalAMG* out=NULL;
62    bool_t verbose=options->verbose;    bool_t verbose=options->verbose;
63        
64    Paso_SparseMatrix* W_FC=NULL, *Atemp=NULL, *A_C=NULL;    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;    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;    dim_t n_F=0, n_C=0, i;
69    double time0=0;    double time0=0;
70    const double theta = options->coarsening_threshold;    const double theta = options->coarsening_threshold;
# Line 84  Paso_Preconditioner_LocalAMG* Paso_Preco Line 84  Paso_Preconditioner_LocalAMG* Paso_Preco
84    
85       split_marker=TMPMEMALLOC(n,index_t);       split_marker=TMPMEMALLOC(n,index_t);
86       counter=TMPMEMALLOC(n,index_t);       counter=TMPMEMALLOC(n,index_t);
87       if ( !( Esys_checkPtr(split_marker) || Esys_checkPtr(counter) ) ) {       degree=TMPMEMALLOC(n, dim_t);
88         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            set splitting of unknows:            set splitting of unknows:
92                
93           */           */
94       time0=Esys_timer();       time0=Esys_timer();
95       if (n_block>1) {       if (n_block>1) {
96          Paso_Preconditioner_AMG_RSCoarsening_Block(A_p, split_marker, theta,tau);             Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree, S, theta,tau);
97       } else {       } else {
98          Paso_Preconditioner_AMG_RSCoarsening(A_p, split_marker, theta,tau);             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree, S, theta,tau);
99       }       }
100         Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree, S, split_marker);
101       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
102            
103       if (Esys_noError() ) {       if (Esys_noError() ) {
# Line 149  Paso_Preconditioner_LocalAMG* Paso_Preco Line 152  Paso_Preconditioner_LocalAMG* Paso_Preco
152                 }                 }
153              }              }
154              /*              /*
155                    get Restriction : (can we do this in one go?)                      get Restriction :  
156              */              */                  
             time0=Esys_timer();  
             W_FC=Paso_SparseMatrix_getSubmatrix(A_p, n_F, n_C, rows_in_F, mask_C);  
             if (SHOW_TIMING) printf("timing: level %d: get Weights: %e\n",level, Esys_timer()-time0);  
               
             time0=Esys_timer();  
             Paso_SparseMatrix_updateWeights(A_p,W_FC,split_marker);  
             if (SHOW_TIMING) printf("timing: level %d: updateWeights: %e\n",level, Esys_timer()-time0);  
               
157              time0=Esys_timer();              time0=Esys_timer();
158              out->P=Paso_SparseMatrix_getProlongation(W_FC,split_marker);              out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);
159              if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);              if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);
160                    /*      
             Paso_SparseMatrix_free(W_FC);  
             /*        
161                 construct Prolongation operator as transposed of restriction operator:                 construct Prolongation operator as transposed of restriction operator:
162              */              */
163              time0=Esys_timer();              time0=Esys_timer();
164              out->R=Paso_SparseMatrix_getTranspose(out->P);              out->R=Paso_SparseMatrix_getTranspose(out->P);
165              if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);              if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
166                        
167              /*              /*
168              constrPaso AMG level 2: 4 unknowns are flagged for elimination.              construct coarse level matrix (can we do this in one call?)
             timing: Gauss-Seidel preparation: elemination : 8.096400e-05  
             AMG: level 2: 4 unknowns eliminated.  
             uct coarse level matrix (can we do this in one call?)  
169              */              */
170              time0=Esys_timer();              time0=Esys_timer();
171              Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);              Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
172              A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);              A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
173              Paso_SparseMatrix_free(Atemp);              Paso_SparseMatrix_free(Atemp);
             /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/  
174              if (SHOW_TIMING) printf("timing: level %d : getCoarseMatrix: %e\n",level,Esys_timer()-time0);              if (SHOW_TIMING) printf("timing: level %d : getCoarseMatrix: %e\n",level,Esys_timer()-time0);
175                            
176              /* allocate helpers :*/              /* allocate helpers :*/
# Line 239  Paso_Preconditioner_LocalAMG* Paso_Preco Line 229  Paso_Preconditioner_LocalAMG* Paso_Preco
229    }    }
230    TMPMEMFREE(counter);    TMPMEMFREE(counter);
231    TMPMEMFREE(split_marker);    TMPMEMFREE(split_marker);
232      TMPMEMFREE(degree);
233      TMPMEMFREE(S);
234    
235    if (Esys_noError()) {    if (Esys_noError()) {
236        if (verbose) printf("AMG: level %d: %d unknowns eliminated.\n",level, n_F);        if (verbose) printf("AMG: level %d: %d unknowns eliminated.\n",level, n_F);
# Line 266  void Paso_Preconditioner_LocalAMG_solve( Line 258  void Paso_Preconditioner_LocalAMG_solve(
258       if (amg->n_F < amg->n) { /* is there work on the coarse level? */       if (amg->n_F < amg->n) { /* is there work on the coarse level? */
259           time0=Esys_timer();           time0=Esys_timer();
260       Paso_Copy(n, amg->r, b);                            /*  r <- b */       Paso_Copy(n, amg->r, b);                            /*  r <- b */
261       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C); /*r=r-Ax*/
262           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */
263           time0=Esys_timer()-time0;           time0=Esys_timer()-time0;
264       /* coarse level solve */       /* coarse level solve */
# Line 289  void Paso_Preconditioner_LocalAMG_solve( Line 281  void Paso_Preconditioner_LocalAMG_solve(
281          Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */          Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */
282       }         }  
283       time0=time0+Esys_timer();       time0=time0+Esys_timer();
284       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,1.,x);       /* x = x + P*x_c */       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */    
285            
286           /*postsmoothing*/           /*postsmoothing*/
287                
# Line 307  void Paso_Preconditioner_LocalAMG_solve( Line 299  void Paso_Preconditioner_LocalAMG_solve(
299  /* theta = threshold for strong connections */  /* theta = threshold for strong connections */
300  /* tau = threshold for diagonal dominance */  /* tau = threshold for diagonal dominance */
301    
302  void Paso_Preconditioner_AMG_RSCoarsening(Paso_SparseMatrix* A, index_t* split_marker, const double theta, const double tau)  /*S_i={j \in N_i; i strongly coupled to j}
303    
304    in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
305    */
306    
307    void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,
308                          dim_t *degree, index_t *S,
309                          const double theta, const double tau)
310  {  {
311     const dim_t n=A->numRows;     const dim_t n=A->numRows;
312     dim_t *degree=NULL;   /* number of naighbours in strong connection graph */     index_t iptr, i,j;
    index_t *S=NULL, iptr, i,j;  
313     dim_t kdeg;     dim_t kdeg;
314     double max_offdiagonal, threshold, sum_row, main_row, fnorm;     double max_offdiagonal, threshold, sum_row, main_row, fnorm;
     
    degree=TMPMEMALLOC(n, dim_t);  
    S=TMPMEMALLOC(A->pattern->len, index_t);  
315    
316     if ( !( Esys_checkPtr(degree) || Esys_checkPtr(S) ) )  {  
       /*S_i={j \in N_i; i strongly coupled to j}  
     
      in the sense that |A_{ij}| >= theta * max_k |A_{ik}|  
       */  
317        #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)        #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)
318        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
319                        
# Line 355  void Paso_Preconditioner_AMG_RSCoarsenin Line 346  void Paso_Preconditioner_AMG_RSCoarsenin
346       }       }
347       degree[i]=kdeg;       degree[i]=kdeg;
348        }        }
349          
       Paso_Preconditioner_AMG_RSCoarsening_search(n, A->pattern->ptr, degree, S, split_marker);  
         
    }  
    TMPMEMFREE(degree);  
    TMPMEMFREE(S);  
350  }  }
351    
352  /* theta = threshold for strong connections */  /* theta = threshold for strong connections */
353  /* tau = threshold for diagonal dominance */  /* tau = threshold for diagonal dominance */
354  void Paso_Preconditioner_AMG_RSCoarsening_Block(Paso_SparseMatrix* A, index_t* split_marker, const double theta, const double tau)  /*S_i={j \in N_i; i strongly coupled to j}
355    
356    in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
357    */
358    void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,
359                                dim_t *degree, index_t *S,
360                                const double theta, const double tau)
361    
362  {  {
363     const dim_t n_block=A->row_block_size;     const dim_t n_block=A->row_block_size;
364     const dim_t n=A->numRows;     const dim_t n=A->numRows;
365     dim_t *degree=NULL;   /* number of naighbours in strong connection graph */     index_t iptr, i,j, bi;
    index_t *S=NULL, iptr, i,j, bi;  
366     dim_t kdeg, max_deg;     dim_t kdeg, max_deg;
367     register double max_offdiagonal, threshold, fnorm, sum_row, main_row;     register double max_offdiagonal, threshold, fnorm, sum_row, main_row;
368     double *rtmp;     double *rtmp;
369        
370        
    degree=TMPMEMALLOC(n, dim_t);  
    S=TMPMEMALLOC(A->pattern->len, index_t);  
     
    if ( !( Esys_checkPtr(degree) || Esys_checkPtr(S)  ) ) {  
       /*S_i={j \in N_i; i strongly coupled to j}  
     
       in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F  
       */  
371        #pragma omp parallel private(i,iptr,max_offdiagonal, kdeg, threshold,j, max_deg, fnorm, sum_row, main_row, rtmp)        #pragma omp parallel private(i,iptr,max_offdiagonal, kdeg, threshold,j, max_deg, fnorm, sum_row, main_row, rtmp)
372        {        {
373       max_deg=0;       max_deg=0;
# Line 430  void Paso_Preconditioner_AMG_RSCoarsenin Line 413  void Paso_Preconditioner_AMG_RSCoarsenin
413       }             }      
414       TMPMEMFREE(rtmp);       TMPMEMFREE(rtmp);
415        } /* end of parallel region */        } /* end of parallel region */
416        Paso_Preconditioner_AMG_RSCoarsening_search(n, A->pattern->ptr, degree, S, split_marker);  
    }  
    TMPMEMFREE(degree);  
    TMPMEMFREE(S);    
417  }    }  
418    
419  /* the runge stueben coarsening algorithm: */  /* the runge stueben coarsening algorithm: */
420  void Paso_Preconditioner_AMG_RSCoarsening_search(const dim_t n, const index_t* offset, const dim_t* degree, const index_t* S,  void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset,
421                             const dim_t* degree, const index_t* S,
422                           index_t*split_marker)                           index_t*split_marker)
423  {  {
424     index_t *lambda=NULL, j, *ST=NULL;     index_t *lambda=NULL, j, *ST=NULL;

Legend:
Removed from v.3287  
changed lines
  Added in v.3303

  ViewVC Help
Powered by ViewVC 1.1.26