/[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 3403 by gross, Tue Dec 7 08:13:51 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         if (in->A_C == NULL) {
64            return 1.;
65         } else {
66            return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);
67         }
68          } else {
69            return Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(in->AMG_C);
70          }
71    }
72    dim_t Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(const Paso_Preconditioner_LocalAMG * in) {
73       if (in->AMG_C == NULL) {
74          if (in->A_C == NULL) {
75         return 0;
76          } else {
77         return in->A_C->numRows;
78          }
79       } else {
80         return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);
81       }
82    }
83  /*****************************************************************  /*****************************************************************
84    
85     constructs AMG     constructs AMG
# Line 76  Paso_Preconditioner_LocalAMG* Paso_Preco Line 105  Paso_Preconditioner_LocalAMG* Paso_Preco
105                
106    */    */
107    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) ) {
108       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",
109      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  );  
110       return NULL;       return NULL;
111    }    }
# Line 109  Paso_Preconditioner_LocalAMG* Paso_Preco Line 138  Paso_Preconditioner_LocalAMG* Paso_Preco
138          */          */
139          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);
140          n_C=n-n_F;          n_C=n-n_F;
141          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);
142            
143          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 */
144             out = NULL;             out = NULL;
# Line 139  Paso_Preconditioner_LocalAMG* Paso_Preco Line 168  Paso_Preconditioner_LocalAMG* Paso_Preco
168    
169            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);
170                
171            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) {
172                   /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
173        
174              /* allocate helpers :*/              /* allocate helpers :*/
175              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 243  Paso_Preconditioner_LocalAMG* Paso_Preco
243                      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);
244                      Paso_SparseMatrix_free(A_C);                      Paso_SparseMatrix_free(A_C);
245                      out->A_C->solver_package = PASO_MKL;                      out->A_C->solver_package = PASO_MKL;
246                      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);
247                    #else                    #else
248                      #ifdef UMFPACK                      #ifdef UMFPACK
249                         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);
250                         Paso_SparseMatrix_free(A_C);                         Paso_SparseMatrix_free(A_C);
251                         out->A_C->solver_package = PASO_UMFPACK;                         out->A_C->solver_package = PASO_UMFPACK;
252                         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);
253                      #else                      #else
254                         out->A_C=A_C;                         out->A_C=A_C;
255                         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);
256                         out->A_C->solver_package = PASO_SMOOTHER;                         out->A_C->solver_package = PASO_SMOOTHER;
257                         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);
258                      #endif                      #endif
259                    #endif                    #endif
260                 } else {                 } else {
# Line 284  void Paso_Preconditioner_LocalAMG_solve( Line 314  void Paso_Preconditioner_LocalAMG_solve(
314            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);
315            break;            break;
316             case (PASO_SMOOTHER):             case (PASO_SMOOTHER):
317                    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);
318            break;            break;
319          }          }
320          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 356  void Paso_Preconditioner_AMG_setStrongCo
356    
357    
358        #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)
359        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {        
             
360       max_offdiagonal = 0.;       max_offdiagonal = 0.;
361       sum_row=0;       sum_row=0;
362       main_row=0;       main_row=0;
# Line 345  void Paso_Preconditioner_AMG_setStrongCo Line 374  void Paso_Preconditioner_AMG_setStrongCo
374       }       }
375       threshold = theta*max_offdiagonal;       threshold = theta*max_offdiagonal;
376       kdeg=0;       kdeg=0;
377        
378       if (tau*main_row < sum_row) { /* no diagonal domainance */       if (tau*main_row < sum_row) { /* no diagonal domainance */
379          #pragma ivdep          #pragma ivdep
380          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.3403

  ViewVC Help
Powered by ViewVC 1.1.26