/[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 3322 by gross, Wed Oct 27 01:20:27 2010 UTC revision 3323 by gross, Thu Oct 28 09:53:46 2010 UTC
# Line 76  Paso_Preconditioner_LocalAMG* Paso_Preco Line 76  Paso_Preconditioner_LocalAMG* Paso_Preco
76                
77    */    */
78    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) ) {
79       if (verbose) printf("Paso: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",       if (verbose) printf("Paso_Solver: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
80      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  );  
81       return NULL;       return NULL;
82    }    }
# Line 109  Paso_Preconditioner_LocalAMG* Paso_Preco Line 109  Paso_Preconditioner_LocalAMG* Paso_Preco
109          */          */
110          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);          n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);
111          n_C=n-n_F;          n_C=n-n_F;
112          if (verbose) printf("Paso: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);          if (verbose) printf("Paso_Solver: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
113            
114          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 */
115             out = NULL;             out = NULL;
116          } else {          } else {
117             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);             out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);
118             mask_C=TMPMEMALLOC(n,index_t);             if (! Esys_checkPtr(out)) {
            rows_in_F=TMPMEMALLOC(n_F,index_t);  
            if ( !( Esys_checkPtr(mask_C) || Esys_checkPtr(rows_in_F) || Esys_checkPtr(out)) ) {  
119            out->level = level;            out->level = level;
120            out->n = n;            out->n = n;
121            out->n_F = n_F;            out->n_F = n_F;
# Line 131  Paso_Preconditioner_LocalAMG* Paso_Preco Line 129  Paso_Preconditioner_LocalAMG* Paso_Preco
129            out->x_C = NULL;            out->x_C = NULL;
130            out->b_C = NULL;            out->b_C = NULL;
131            out->AMG_C = NULL;            out->AMG_C = NULL;
132              out->Smoother=NULL;
133               }
134               mask_C=TMPMEMALLOC(n,index_t);
135               rows_in_F=TMPMEMALLOC(n_F,index_t);
136               Esys_checkPtr(mask_C);
137               Esys_checkPtr(rows_in_F);
138               if ( Esys_noError() ) {
139    
140            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);
141                    
142            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_F < n ) { /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
143        
144              /* allocate helpers :*/              /* allocate helpers :*/
# Line 141  Paso_Preconditioner_LocalAMG* Paso_Preco Line 147  Paso_Preconditioner_LocalAMG* Paso_Preco
147              out->r=MEMALLOC(n_block*n,double);              out->r=MEMALLOC(n_block*n,double);
148                            
149              Esys_checkPtr(out->r);              Esys_checkPtr(out->r);
             Esys_checkPtr(out->Smoother);  
150              Esys_checkPtr(out->x_C);              Esys_checkPtr(out->x_C);
151              Esys_checkPtr(out->b_C);              Esys_checkPtr(out->b_C);
152                            
153              /* creates index for F:*/              if ( Esys_noError() ) {
154              #pragma omp parallel for private(i) schedule(static)                 /* creates index for F:*/
155              for (i = 0; i < n; ++i) {                 #pragma omp parallel private(i)
156                 if  (split_marker[i]) rows_in_F[counter[i]]=i;                 {
157              }                              #pragma omp for schedule(static)
158              /*  create mask of C nodes with value >-1 gives new id */                    for (i = 0; i < n; ++i) {
159              i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);                   if  (split_marker[i]) rows_in_F[counter[i]]=i;
160                      }
             #pragma omp parallel for private(i) schedule(static)  
             for (i = 0; i < n; ++i) {  
                if  (split_marker[i]) {  
                   mask_C[i]=-1;  
                } else {  
                   mask_C[i]=counter[i];;  
161                 }                 }
162              }                 /*  create mask of C nodes with value >-1 gives new id */
163              /*                 i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);
164    
165                   #pragma omp parallel for private(i) schedule(static)
166                   for (i = 0; i < n; ++i) {
167                      if  (split_marker[i]) {
168                     mask_C[i]=-1;
169                      } else {
170                     mask_C[i]=counter[i];;
171                      }
172                   }
173                   /*
174                    get Restriction :                      get Restriction :  
175              */                                   */                  
176              time0=Esys_timer();                 time0=Esys_timer();
177              out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);                 out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);
178              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);
179      /*                    }
180                /*      
181                 construct Prolongation operator as transposed of restriction operator:                 construct Prolongation operator as transposed of restriction operator:
182              */              */
183              if ( Esys_noError()) {              if ( Esys_noError()) {
# Line 203  Paso_Preconditioner_LocalAMG* Paso_Preco Line 213  Paso_Preconditioner_LocalAMG* Paso_Preco
213                      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);
214                      Paso_SparseMatrix_free(A_C);                      Paso_SparseMatrix_free(A_C);
215                      out->A_C->solver_package = PASO_MKL;                      out->A_C->solver_package = PASO_MKL;
216                      if (verbose) printf("Paso: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C);                      if (verbose) printf("Paso_Solver: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C);
217                    #else                    #else
218                      #ifdef UMFPACK                      #ifdef UMFPACK
219                         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);
220                         Paso_SparseMatrix_free(A_C);                         Paso_SparseMatrix_free(A_C);
221                         out->A_C->solver_package = PASO_UMFPACK;                         out->A_C->solver_package = PASO_UMFPACK;
222                         if (verbose) printf("Paso: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C);                         if (verbose) printf("Paso_Solver: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C);
223                      #else                      #else
224                         out->A_C=A_C;                         out->A_C=A_C;
225                         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);
226                         out->A_C->solver_package = PASO_SMOOTHER;                         out->A_C->solver_package = PASO_SMOOTHER;
227                         if (verbose) printf("Paso: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C);                         if (verbose) printf("Paso_Solver: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C);
228                      #endif                      #endif
229                    #endif                    #endif
230                 } else {                 } else {
# Line 424  void Paso_Preconditioner_AMG_RungeStuebe Line 434  void Paso_Preconditioner_AMG_RungeStuebe
434  {  {
435     index_t *lambda=NULL, j, *ST=NULL;     index_t *lambda=NULL, j, *ST=NULL;
436     dim_t i,k, p, q, *degreeT=NULL, itmp;     dim_t i,k, p, q, *degreeT=NULL, itmp;
437       double time0=0;
438        
439     if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */     if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
440        
# Line 432  void Paso_Preconditioner_AMG_RungeStuebe Line 443  void Paso_Preconditioner_AMG_RungeStuebe
443     ST=TMPMEMALLOC(offset[n], index_t);     ST=TMPMEMALLOC(offset[n], index_t);
444        
445     if (! ( Esys_checkPtr(lambda) || Esys_checkPtr(degreeT) || Esys_checkPtr(ST) ) ) {     if (! ( Esys_checkPtr(lambda) || Esys_checkPtr(degreeT) || Esys_checkPtr(ST) ) ) {
446    time0=Esys_timer();
447    
448        /* initialize  split_marker and split_marker :*/        /* initialize  split_marker and split_marker :*/
449        /* those unknows which are not influenced go into F, the rest is available for F or C */        /* those unknows which are not influenced go into F, the rest is available for F or C */
450        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
# Line 469  void Paso_Preconditioner_AMG_RungeStuebe Line 482  void Paso_Preconditioner_AMG_RungeStuebe
482          lambda[i]=itmp;          lambda[i]=itmp;
483       }       }
484        }        }
485          /* printf("timing: %e\n",Esys_timer()-time0); */
486          time0=Esys_timer();
487        /* start search :*/        /* start search :*/
488        i=Paso_Util_arg_max(n,lambda);        i=Paso_Util_arg_max(n,lambda);
489        while (lambda[i]>-1) { /* is there any undecided unknowns? */        while (lambda[i]>-1) { /* is there any undecided unknowns? */
# Line 502  void Paso_Preconditioner_AMG_RungeStuebe Line 516  void Paso_Preconditioner_AMG_RungeStuebe
516       i=Paso_Util_arg_max(n,lambda);       i=Paso_Util_arg_max(n,lambda);
517        }        }
518     }     }
519       /* printf("timing: %e\n",Esys_timer()-time0); */
520     TMPMEMFREE(lambda);     TMPMEMFREE(lambda);
521     TMPMEMFREE(ST);     TMPMEMFREE(ST);
522     TMPMEMFREE(degreeT);     TMPMEMFREE(degreeT);

Legend:
Removed from v.3322  
changed lines
  Added in v.3323

  ViewVC Help
Powered by ViewVC 1.1.26