/[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 3396 by gross, Fri Dec 3 08:24:16 2010 UTC revision 3402 by gross, Tue Dec 7 07:36:12 2010 UTC
# 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) {
55       if (in->AMG_C == NULL) {
56          return in->level;
57       } else {
58          return Paso_Preconditioner_LocalAMG_getMaxLevel(in->AMG_C);
59       }
60    }
61    double Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(const Paso_Preconditioner_LocalAMG * in) {
62       if (in->AMG_C == NULL) {
63          return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);
64       } else {
65          return Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(in->AMG_C);
66       }
67    }
68    dim_t Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(const Paso_Preconditioner_LocalAMG * in) {
69       if (in->AMG_C == NULL) {
70          return in->A_C->numRows;
71       } else {
72          return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);
73       }
74    }
75  /*****************************************************************  /*****************************************************************
76    
77     constructs AMG     constructs AMG
# Line 76  Paso_Preconditioner_LocalAMG* Paso_Preco Line 97  Paso_Preconditioner_LocalAMG* Paso_Preco
97                
98    */    */
99    if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {    if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {
100       if (verbose) printf("Paso_Solver: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",       if (verbose) printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
101      level,  options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size  );        level,  options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size  );  
102       return NULL;       return NULL;
103    }    }
# Line 109  Paso_Preconditioner_LocalAMG* Paso_Preco Line 130  Paso_Preconditioner_LocalAMG* Paso_Preco
130          */          */
131          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);
132          n_C=n-n_F;          n_C=n-n_F;
133          if (verbose) printf("Paso_Solver: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);          if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
134            
135          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 */
136             out = NULL;             out = NULL;
# Line 139  Paso_Preconditioner_LocalAMG* Paso_Preco Line 160  Paso_Preconditioner_LocalAMG* Paso_Preco
160    
161            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
162                
163            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 */            if (n_C != 0) {
164                   /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
165        
166              /* allocate helpers :*/              /* allocate helpers :*/
167              out->x_C=MEMALLOC(n_block*n_C,double);              out->x_C=MEMALLOC(n_block*n_C,double);
# Line 213  Paso_Preconditioner_LocalAMG* Paso_Preco Line 235  Paso_Preconditioner_LocalAMG* Paso_Preco
235                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);                      out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
236                      Paso_SparseMatrix_free(A_C);                      Paso_SparseMatrix_free(A_C);
237                      out->A_C->solver_package = PASO_MKL;                      out->A_C->solver_package = PASO_MKL;
238                      if (verbose) printf("Paso_Solver: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C);                      if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
239                    #else                    #else
240                      #ifdef UMFPACK                      #ifdef UMFPACK
241                         out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);                         out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
242                         Paso_SparseMatrix_free(A_C);                         Paso_SparseMatrix_free(A_C);
243                         out->A_C->solver_package = PASO_UMFPACK;                         out->A_C->solver_package = PASO_UMFPACK;
244                         if (verbose) printf("Paso_Solver: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C);                         if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
245                      #else                      #else
246                         out->A_C=A_C;                         out->A_C=A_C;
247                         out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);                         out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);
248                         out->A_C->solver_package = PASO_SMOOTHER;                         out->A_C->solver_package = PASO_SMOOTHER;
249                         if (verbose) printf("Paso_Solver: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C);                         if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
250                      #endif                      #endif
251                    #endif                    #endif
252                 } else {                 } else {
# Line 284  void Paso_Preconditioner_LocalAMG_solve( Line 306  void Paso_Preconditioner_LocalAMG_solve(
306            Paso_UMFPACK(amg->A_C, amg->x_C,amg->b_C, amg->refinements, SHOW_TIMING);            Paso_UMFPACK(amg->A_C, amg->x_C,amg->b_C, amg->refinements, SHOW_TIMING);
307            break;            break;
308             case (PASO_SMOOTHER):             case (PASO_SMOOTHER):
309                    Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);            Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);
310            break;            break;
311          }          }
312          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);
# Line 326  void Paso_Preconditioner_AMG_setStrongCo Line 348  void Paso_Preconditioner_AMG_setStrongCo
348    
349    
350        #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)
351        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {        
             
352       max_offdiagonal = 0.;       max_offdiagonal = 0.;
353       sum_row=0;       sum_row=0;
354       main_row=0;       main_row=0;
# Line 345  void Paso_Preconditioner_AMG_setStrongCo Line 366  void Paso_Preconditioner_AMG_setStrongCo
366       }       }
367       threshold = theta*max_offdiagonal;       threshold = theta*max_offdiagonal;
368       kdeg=0;       kdeg=0;
369        
370       if (tau*main_row < sum_row) { /* no diagonal domainance */       if (tau*main_row < sum_row) { /* no diagonal domainance */
371          #pragma ivdep          #pragma ivdep
372          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {

Legend:
Removed from v.3396  
changed lines
  Added in v.3402

  ViewVC Help
Powered by ViewVC 1.1.26