/[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 3314 by gross, Tue Oct 26 22:22:20 2010 UTC revision 3351 by gross, Tue Nov 16 02:39:40 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        
              /* creates index for F:*/  
             #pragma omp parallel for private(i) schedule(static)  
             for (i = 0; i < n; ++i) {  
                if  (split_marker[i]) rows_in_F[counter[i]]=i;  
             }            
             /*  create mask of C nodes with value >-1 gives new id */  
             i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);  
   
             #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];;  
                }  
             }  
             /*  
                   get Restriction :    
             */                    
             time0=Esys_timer();  
             out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);  
             if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);  
     /*        
                construct Prolongation operator as transposed of restriction operator:  
             */  
             time0=Esys_timer();  
             out->R=Paso_SparseMatrix_getTranspose(out->P);  
             if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);  
                       
             /*  
             construct coarse level matrix:  
             */  
             time0=Esys_timer();  
             Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);  
             A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);  
             Paso_SparseMatrix_free(Atemp);  
             if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);  
               
144              /* allocate helpers :*/              /* allocate helpers :*/
145              out->x_C=MEMALLOC(n_block*n_C,double);              out->x_C=MEMALLOC(n_block*n_C,double);
146              out->b_C=MEMALLOC(n_block*n_C,double);              out->b_C=MEMALLOC(n_block*n_C,double);
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                if ( Esys_noError() ) {
154                   /* creates index for F:*/
155                   #pragma omp parallel private(i)
156                   {
157                      #pragma omp for schedule(static)
158                      for (i = 0; i < n; ++i) {
159                     if  (split_marker[i]) rows_in_F[counter[i]]=i;
160                      }
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 :  
175                   */                  
176                   time0=Esys_timer();
177                   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);
179                }
180                /*      
181                   construct Prolongation operator as transposed of restriction operator:
182                */
183                if ( Esys_noError()) {
184                   time0=Esys_timer();
185                   out->R=Paso_SparseMatrix_getTranspose(out->P);
186                   if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
187                }      
188                /*
189                construct coarse level matrix:
190                */
191                if ( Esys_noError()) {
192                   time0=Esys_timer();
193                   Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
194                   A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
195                   Paso_SparseMatrix_free(Atemp);
196                   if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);            
197                }
198    
199                            
200              /*              /*
201                 constructe courser level:                 constructe courser level:
202                                
203              */              */
204              out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);              if ( Esys_noError()) {
205                               out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);
206                }
207              if ( Esys_noError()) {              if ( Esys_noError()) {
208                 if ( out->AMG_C == NULL ) {                 if ( out->AMG_C == NULL ) {
209                    out->reordering = options->reordering;                    out->reordering = options->reordering;
# Line 198  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 417  void Paso_Preconditioner_AMG_RungeStuebe Line 432  void Paso_Preconditioner_AMG_RungeStuebe
432                           const dim_t* degree, const index_t* S,                           const dim_t* degree, const index_t* S,
433                           index_t*split_marker)                           index_t*split_marker)
434  {  {
435     index_t *lambda=NULL, j, *ST=NULL;     const bool_t usePanel=FALSE;
436     dim_t i,k, p, q, *degreeT=NULL, itmp;    
437       index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;
438       dim_t i,k, p, q, *degreeT=NULL, len_panel, len_panel_new;
439       register index_t j, itmp;
440        
441     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 */
442        
443     lambda=TMPMEMALLOC(n, index_t);     lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);
444     degreeT=TMPMEMALLOC(n, dim_t);     degreeT=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degreeT);
445     ST=TMPMEMALLOC(offset[n], index_t);     ST=TMPMEMALLOC(offset[n], index_t);  Esys_checkPtr(ST);
446       if (usePanel) {
447          notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);
448          panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);
449       }
450      
451        
452     if (! ( Esys_checkPtr(lambda) || Esys_checkPtr(degreeT) || Esys_checkPtr(ST) ) ) {    
453       if (Esys_noError() ) {
454        /* initialize  split_marker and split_marker :*/        /* initialize  split_marker and split_marker :*/
455        /* 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 */
456        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
# Line 464  void Paso_Preconditioner_AMG_RungeStuebe Line 488  void Paso_Preconditioner_AMG_RungeStuebe
488          lambda[i]=itmp;          lambda[i]=itmp;
489       }       }
490        }        }
491          if (usePanel) {
492         #pragma omp parallel for private(i) schedule(static)
493         for (i=0;i<n;++i) notInPanel[i]=TRUE;
494          }
495        /* start search :*/        /* start search :*/
496        i=Paso_Util_arg_max(n,lambda);        i=Paso_Util_arg_max(n,lambda);
497        while (lambda[i]>-1) { /* is there any undecided unknowns? */        while (lambda[i]>-1) { /* is there any undecided unknown? */
498        
499       /* the unknown i is moved to C */       if (usePanel) {
500       split_marker[i]=PASO_AMG_IN_C;          len_panel=0;
501       lambda[i]=-1;  /* lambda fro unavailable unknowns is set to -1 */          do {
502                   /* the unknown i is moved to C */
503       /* all undecided unknown strongly coupled to i are moved to F */             split_marker[i]=PASO_AMG_IN_C;
504       for (p=0; p<degreeT[i]; ++p) {             lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
505          j=ST[offset[i]+p];            
506                   /* all undecided unknown strongly coupled to i are moved to F */
507          if (split_marker[j]==PASO_AMG_UNDECIDED) {             for (p=0; p<degreeT[i]; ++p) {
508                    j=ST[offset[i]+p];
509             split_marker[j]=PASO_AMG_IN_F;            
510             lambda[j]=-1;            if (split_marker[j]==PASO_AMG_UNDECIDED) {
511                      
512             for (q=0; q<degreeT[j]; ++q) {               split_marker[j]=PASO_AMG_IN_F;
513            k=ST[offset[j]+q];               lambda[j]=-1;
514            if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;              
515                 for (q=0; q<degreeT[j]; ++q) {
516                k=ST[offset[j]+q];
517                if (split_marker[k]==PASO_AMG_UNDECIDED) {
518                   lambda[k]++;
519                   if (notInPanel[k]) {
520                      notInPanel[k]=FALSE;
521                      panel[len_panel]=k;
522                      len_panel++;
523                   }
524    
525                }    /* the unknown i is moved to C */
526                split_marker[i]=PASO_AMG_IN_C;
527                lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
528                 }
529                
530              }
531               }
532               for (p=0; p<degree[i]; ++p) {
533              j=S[offset[i]+p];
534              if (split_marker[j]==PASO_AMG_UNDECIDED) {
535                 lambda[j]--;
536                 if (notInPanel[j]) {
537                notInPanel[j]=FALSE;
538                panel[len_panel]=j;
539                len_panel++;
540                 }
541              }
542             }             }
543    
544               /* consolidate panel */
545               /* remove lambda[q]=-1 */
546               lambda_max=-1;
547               i=-1;
548               len_panel_new=0;
549               for (q=0; q<len_panel; q++) {
550                 k=panel[q];
551                 lambda_k=lambda[k];
552                 if (split_marker[k]==PASO_AMG_UNDECIDED) {
553                panel[len_panel_new]=k;
554                len_panel_new++;
555    
556                if (lambda_max == lambda_k) {
557                   if (k<i) i=k;
558                } else if (lambda_max < lambda_k) {
559                   lambda_max =lambda_k;
560                   i=k;
561                }
562                 }
563               }
564               len_panel=len_panel_new;
565            } while (len_panel>0);    
566         } else {
567            /* the unknown i is moved to C */
568            split_marker[i]=PASO_AMG_IN_C;
569            lambda[i]=-1;  /* lambda from unavailable unknowns is set to -1 */
570            
571            /* all undecided unknown strongly coupled to i are moved to F */
572            for (p=0; p<degreeT[i]; ++p) {
573               j=ST[offset[i]+p];
574               if (split_marker[j]==PASO_AMG_UNDECIDED) {
575            
576              split_marker[j]=PASO_AMG_IN_F;
577              lambda[j]=-1;
578            
579              for (q=0; q<degreeT[j]; ++q) {
580                 k=ST[offset[j]+q];
581                 if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;
582              }
583    
584               }
585          }          }
586            for (p=0; p<degree[i]; ++p) {
587               j=S[offset[i]+p];
588               if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;
589            }
590            
591       }       }
      for (p=0; p<degree[i]; ++p) {  
         j=S[offset[i]+p];  
         if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;  
      }  
       
592       i=Paso_Util_arg_max(n,lambda);       i=Paso_Util_arg_max(n,lambda);
593        }        }
594     }     }
595     TMPMEMFREE(lambda);     TMPMEMFREE(lambda);
596     TMPMEMFREE(ST);     TMPMEMFREE(ST);
597     TMPMEMFREE(degreeT);     TMPMEMFREE(degreeT);
598       TMPMEMFREE(panel);
599       TMPMEMFREE(notInPanel);
600  }  }

Legend:
Removed from v.3314  
changed lines
  Added in v.3351

  ViewVC Help
Powered by ViewVC 1.1.26