/[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 3440 by gross, Fri Jan 14 00:04:53 2011 UTC revision 3441 by gross, Fri Jan 14 01:09:09 2011 UTC
# Line 14  Line 14 
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: AMG preconditioner                                  */  /* Paso: AMG preconditioner  (local version)                  */
18    
19  /**************************************************************/  /**************************************************************/
20    
# Line 35  Line 35 
35    
36  /* free all memory used by AMG                                */  /* free all memory used by AMG                                */
37    
38  void Paso_Preconditioner_LocalAMG_free(Paso_Preconditioner_LocalAMG * in) {  void Paso_Preconditioner_AMG_free(Paso_Preconditioner_AMG * in) {
39       if (in!=NULL) {       if (in!=NULL) {
40      Paso_Preconditioner_LocalSmoother_free(in->Smoother);      Paso_Preconditioner_Smoother_free(in->Smoother);
41      Paso_SparseMatrix_free(in->P);      Paso_SystemMatrix_free(in->P);
42      Paso_SparseMatrix_free(in->R);      Paso_SystemMatrix_free(in->R);
43      Paso_SparseMatrix_free(in->A_C);      Paso_SystemMatrix_free(in->A_C);
44      Paso_Preconditioner_LocalAMG_free(in->AMG_C);      Paso_Preconditioner_AMG_free(in->AMG_C);
45      MEMFREE(in->r);      MEMFREE(in->r);
46      MEMFREE(in->x_C);      MEMFREE(in->x_C);
47      MEMFREE(in->b_C);      MEMFREE(in->b_C);
# Line 51  void Paso_Preconditioner_LocalAMG_free(P Line 51  void Paso_Preconditioner_LocalAMG_free(P
51       }       }
52  }  }
53    
54  index_t Paso_Preconditioner_LocalAMG_getMaxLevel(const Paso_Preconditioner_LocalAMG * in) {  index_t Paso_Preconditioner_AMG_getMaxLevel(const Paso_Preconditioner_AMG * in) {
55     if (in->AMG_C == NULL) {     if (in->AMG_C == NULL) {
56        return in->level;        return in->level;
57     } else {     } else {
58        return Paso_Preconditioner_LocalAMG_getMaxLevel(in->AMG_C);        return Paso_Preconditioner_AMG_getMaxLevel(in->AMG_C);
59     }     }
60  }  }
61  double Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(const Paso_Preconditioner_LocalAMG * in) {  double Paso_Preconditioner_AMG_getCoarseLevelSparsity(const Paso_Preconditioner_AMG * in) {
62        if (in->AMG_C == NULL) {        if (in->AMG_C == NULL) {
63       if (in->A_C == NULL) {       if (in->A_C == NULL) {
64          return 1.;          return 1.;
# Line 66  double Paso_Preconditioner_LocalAMG_getC Line 66  double Paso_Preconditioner_LocalAMG_getC
66          return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);          return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);
67       }       }
68        } else {        } else {
69          return Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(in->AMG_C);          return Paso_Preconditioner_AMG_getCoarseLevelSparsity(in->AMG_C);
70        }        }
71  }  }
72  dim_t Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(const Paso_Preconditioner_LocalAMG * in) {  dim_t Paso_Preconditioner_AMG_getNumCoarseUnknwons(const Paso_Preconditioner_AMG * in) {
73     if (in->AMG_C == NULL) {     if (in->AMG_C == NULL) {
74        if (in->A_C == NULL) {        if (in->A_C == NULL) {
75       return 0;       return 0;
# Line 77  dim_t Paso_Preconditioner_LocalAMG_getNu Line 77  dim_t Paso_Preconditioner_LocalAMG_getNu
77       return in->A_C->numRows;       return in->A_C->numRows;
78        }        }
79     } else {     } else {
80       return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);       return Paso_Preconditioner_AMG_getNumCoarseUnknwons(in->AMG_C);
81     }     }
82  }  }
83  /*****************************************************************  /*****************************************************************
# Line 85  dim_t Paso_Preconditioner_LocalAMG_getNu Line 85  dim_t Paso_Preconditioner_LocalAMG_getNu
85     constructs AMG     constructs AMG
86        
87  ******************************************************************/  ******************************************************************/
88  Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {  Paso_Preconditioner_AMG* Paso_Preconditioner_AMG_alloc(Paso_SystemMatrix *A_p,dim_t level,Paso_Options* options) {
89    
90    Paso_Preconditioner_LocalAMG* out=NULL;    Paso_Preconditioner_AMG* out=NULL;
91    bool_t verbose=options->verbose;    bool_t verbose=options->verbose;
92        
93    Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;    Paso_SystemMatrix *Atemp=NULL, *A_C=NULL;
94    const dim_t n=A_p->numRows;    const dim_t n=A_p->numRows;
95    const dim_t n_block=A_p->row_block_size;    const dim_t n_block=A_p->row_block_size;
96    index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree_S=NULL;    index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree_S=NULL;
# Line 171  Paso_Preconditioner_LocalAMG* Paso_Preco Line 171  Paso_Preconditioner_LocalAMG* Paso_Preco
171          if ( n_F == 0 ) {  /*  is a nasty case. a direct solver should be used, return NULL */          if ( n_F == 0 ) {  /*  is a nasty case. a direct solver should be used, return NULL */
172             out = NULL;             out = NULL;
173          } else {          } else {
174             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);             out=MEMALLOC(1,Paso_Preconditioner_AMG);
175             if (! Esys_checkPtr(out)) {             if (! Esys_checkPtr(out)) {
176            out->level = level;            out->level = level;
177            out->n = n;            out->n = n;
# Line 240  Paso_Preconditioner_LocalAMG* Paso_Preco Line 240  Paso_Preconditioner_LocalAMG* Paso_Preco
240              */              */
241              if ( Esys_noError()) {              if ( Esys_noError()) {
242                 time0=Esys_timer();                 time0=Esys_timer();
243                 out->R=Paso_SparseMatrix_getTranspose(out->P);                 out->R=Paso_SystemMatrix_getTranspose(out->P);
244                 if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);                 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
245              }                    }      
246              /*              /*
247              construct coarse level matrix:              construct coarse level matrix:
248              */              */
249              if ( Esys_noError()) {              if ( Esys_noError()) {
250                 time0=Esys_timer();                 time0=Esys_timer();
251                 Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);                 Atemp=Paso_SystemMatrix_MatrixMatrix(A_p,out->P);
252                 A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);                 A_C=Paso_SystemMatrix_MatrixMatrix(out->R,Atemp);
253                 Paso_SparseMatrix_free(Atemp);                 Paso_SystemMatrix_free(Atemp);
254                 if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);                             if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);            
255              }              }
256    
# Line 260  Paso_Preconditioner_LocalAMG* Paso_Preco Line 260  Paso_Preconditioner_LocalAMG* Paso_Preco
260                                
261              */              */
262              if ( Esys_noError()) {              if ( Esys_noError()) {
263                 out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);                 out->AMG_C=Paso_Preconditioner_AMG_alloc(A_C,level+1,options);
264              }              }
265              if ( Esys_noError()) {              if ( Esys_noError()) {
266                 if ( out->AMG_C == NULL ) {                 if ( out->AMG_C == NULL ) {
# Line 268  Paso_Preconditioner_LocalAMG* Paso_Preco Line 268  Paso_Preconditioner_LocalAMG* Paso_Preco
268                    out->refinements = options->coarse_matrix_refinements;                    out->refinements = options->coarse_matrix_refinements;
269                    /* no coarse level matrix has been constructed. use direct solver */                    /* no coarse level matrix has been constructed. use direct solver */
270                    #ifdef MKL                    #ifdef MKL
271                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);                      out->A_C=Paso_SystemMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
272                      Paso_SparseMatrix_free(A_C);                      Paso_SystemMatrix_free(A_C);
273                      out->A_C->solver_package = PASO_MKL;                      out->A_C->solver_package = PASO_MKL;
274                      if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);                      if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
275                    #else                    #else
276                      #ifdef UMFPACK                      #ifdef UMFPACK
277                         out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);                         out->A_C=Paso_SystemMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
278                         Paso_SparseMatrix_free(A_C);                         Paso_SystemMatrix_free(A_C);
279                         out->A_C->solver_package = PASO_UMFPACK;                         out->A_C->solver_package = PASO_UMFPACK;
280                         if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);                         if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
281                      #else                      #else
# Line 305  Paso_Preconditioner_LocalAMG* Paso_Preco Line 305  Paso_Preconditioner_LocalAMG* Paso_Preco
305    if (Esys_noError()) {    if (Esys_noError()) {
306       return out;       return out;
307    } else  {    } else  {
308       Paso_Preconditioner_LocalAMG_free(out);       Paso_Preconditioner_AMG_free(out);
309       return NULL;       return NULL;
310    }    }
311  }  }
312    
313    
314  void Paso_Preconditioner_LocalAMG_solve(Paso_SparseMatrix* A, Paso_Preconditioner_LocalAMG * amg, double * x, double * b) {  void Paso_Preconditioner_AMG_solve(Paso_SystemMatrix* A, Paso_Preconditioner_AMG * amg, double * x, double * b) {
315       const dim_t n = amg->n * amg->n_block;       const dim_t n = amg->n * amg->n_block;
316       double time0=0;       double time0=0;
317       const dim_t post_sweeps=amg->post_sweeps;       const dim_t post_sweeps=amg->post_sweeps;
# Line 327  void Paso_Preconditioner_LocalAMG_solve( Line 327  void Paso_Preconditioner_LocalAMG_solve(
327       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? */
328           time0=Esys_timer();           time0=Esys_timer();
329       Paso_Copy(n, amg->r, b);                            /*  r <- b */       Paso_Copy(n, amg->r, b);                            /*  r <- b */
330       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
331       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C);  /* b_c = R*r  */
332           time0=Esys_timer()-time0;           time0=Esys_timer()-time0;
333       /* coarse level solve */       /* coarse level solve */
334       if ( amg->AMG_C == NULL) {       if ( amg->AMG_C == NULL) {
# Line 347  void Paso_Preconditioner_LocalAMG_solve( Line 347  void Paso_Preconditioner_LocalAMG_solve(
347          }          }
348          if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);          if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
349       } else {       } else {
350          Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */          Paso_Preconditioner_AMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C)     */
351       }         }  
352       time0=time0+Esys_timer();       time0=time0+Esys_timer();
353       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */    
354            
355           /*postsmoothing*/           /*postsmoothing*/
356                
# Line 373  void Paso_Preconditioner_LocalAMG_solve( Line 373  void Paso_Preconditioner_LocalAMG_solve(
373  in the sense that |A_{ij}| >= theta * max_k |A_{ik}|  in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
374  */  */
375    
376  void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,  void Paso_Preconditioner_AMG_setStrongConnections(Paso_SystemMatrix* A,
377                        dim_t *degree_S, index_t *S,                        dim_t *degree_S, index_t *S,
378                        const double theta, const double tau)                        const double theta, const double tau)
379  {  {
# Line 424  void Paso_Preconditioner_AMG_setStrongCo Line 424  void Paso_Preconditioner_AMG_setStrongCo
424    
425  in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F  in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
426  */  */
427  void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,  void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SystemMatrix* A,
428                              dim_t *degree_S, index_t *S,                              dim_t *degree_S, index_t *S,
429                              const double theta, const double tau)                              const double theta, const double tau)
430    

Legend:
Removed from v.3440  
changed lines
  Added in v.3441

  ViewVC Help
Powered by ViewVC 1.1.26