/[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 3441 by gross, Fri Jan 14 01:09:09 2011 UTC revision 3442 by gross, Tue Jan 18 01:47:36 2011 UTC
# Line 46  void Paso_Preconditioner_AMG_free(Paso_P Line 46  void Paso_Preconditioner_AMG_free(Paso_P
46      MEMFREE(in->x_C);      MEMFREE(in->x_C);
47      MEMFREE(in->b_C);      MEMFREE(in->b_C);
48    
       
49      MEMFREE(in);      MEMFREE(in);
50       }       }
51  }  }
# Line 63  double Paso_Preconditioner_AMG_getCoarse Line 62  double Paso_Preconditioner_AMG_getCoarse
62       if (in->A_C == NULL) {       if (in->A_C == NULL) {
63          return 1.;          return 1.;
64       } else {       } else {
65          return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);          return Paso_SystemMatrix_getSparsity(in->A_C);
66       }       }
67        } else {        } else {
68          return Paso_Preconditioner_AMG_getCoarseLevelSparsity(in->AMG_C);          return Paso_Preconditioner_AMG_getCoarseLevelSparsity(in->AMG_C);
# Line 74  dim_t Paso_Preconditioner_AMG_getNumCoar Line 73  dim_t Paso_Preconditioner_AMG_getNumCoar
73        if (in->A_C == NULL) {        if (in->A_C == NULL) {
74       return 0;       return 0;
75        } else {        } else {
76       return in->A_C->numRows;       Paso_SystemMatrix_getTotalNumRows(in->A_C);
77        }        }
78     } else {     } else {
79       return Paso_Preconditioner_AMG_getNumCoarseUnknwons(in->AMG_C);       return Paso_Preconditioner_AMG_getNumCoarseUnknwons(in->AMG_C);
# Line 98  Paso_Preconditioner_AMG* Paso_Preconditi Line 97  Paso_Preconditioner_AMG* Paso_Preconditi
97    double time0=0;    double time0=0;
98    const double theta = options->coarsening_threshold;    const double theta = options->coarsening_threshold;
99    const double tau = options->diagonal_dominance_threshold;    const double tau = options->diagonal_dominance_threshold;
100      const double sparsity=Paso_SystemMatrix_getSparsity(A_p);
101      const dim_t total_n=Paso_SystemMatrix_getGlobalTotalNumRows(A_p);
102    
103        
104    /*    /*
105        is the input matrix A suitable for coarsening        is the input matrix A suitable for coarsening
106                
107    */    */
108    if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {    if ( (sparsity >= options->min_coarse_sparsity) ||
109           (total_n <= options->min_coarse_matrix_size) ||
110           (level > options->level_max) ) {
111    
112      if (verbose) {      if (verbose) {
113            /*            /*
# Line 113  Paso_Preconditioner_AMG* Paso_Preconditi Line 116  Paso_Preconditioner_AMG* Paso_Preconditi
116                        - 'SIZE' = min_coarse_matrix_size exceeded                        - 'SIZE' = min_coarse_matrix_size exceeded
117                        - 'LEVEL' = level_max exceeded                        - 'LEVEL' = level_max exceeded
118            */            */
119            printf("Paso_Preconditioner AMG: termination of coarsening by ");            printf("Paso_Preconditioner: AMG: termination of coarsening by ");
120    
121            if (A_p->pattern->len >= options->min_coarse_sparsity * n * n)            if (sparsity >= options->min_coarse_sparsity)
122                printf("SPAR ");                printf("SPAR ");
123    
124            if (n <= options->min_coarse_matrix_size)            if (total_n <= options->min_coarse_matrix_size)
125                printf("SIZE ");                printf("SIZE ");
126    
127            if (level > options->level_max)            if (level > options->level_max)
# Line 127  Paso_Preconditioner_AMG* Paso_Preconditi Line 130  Paso_Preconditioner_AMG* Paso_Preconditi
130            printf("\n");            printf("\n");
131    
132          printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",          printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
133                      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, sparsity, options->min_coarse_sparsity, total_n, options->min_coarse_matrix_size);  
134    
135      }      }
136    
# Line 150  Paso_Preconditioner_AMG* Paso_Preconditi Line 153  Paso_Preconditioner_AMG* Paso_Preconditi
153       } else {       } else {
154             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree_S, S, theta,tau);             Paso_Preconditioner_AMG_setStrongConnections(A_p, degree_S, S, theta,tau);
155       }       }
      Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);  
156    
157    /*MPI:
158         Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);
159    */
160        
161           /* in BoomerAMG interpolation is used FF connectiovity is required :*/           /* in BoomerAMG interpolation is used FF connectiovity is required :*/
162    /*MPI:
163           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)           if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
164                               Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);                                 Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);  
165    */
166    
167       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);       options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
168            
169       if (Esys_noError() ) {       if (Esys_noError() ) {
# Line 194  Paso_Preconditioner_AMG* Paso_Preconditi Line 203  Paso_Preconditioner_AMG* Paso_Preconditi
203             Esys_checkPtr(rows_in_F);             Esys_checkPtr(rows_in_F);
204             if ( Esys_noError() ) {             if ( Esys_noError() ) {
205    
206            out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);            out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
207                
208            if (n_C != 0) {            if (n_C != 0) {
209                 /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */                 /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
# Line 232  Paso_Preconditioner_AMG* Paso_Preconditi Line 241  Paso_Preconditioner_AMG* Paso_Preconditi
241                    get Prolongation :                        get Prolongation :    
242                 */                                   */                  
243                 time0=Esys_timer();                 time0=Esys_timer();
244    /*MPI:
245                 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);                 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);
246    */
247                 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);
248              }              }
249              /*                    /*      
# Line 240  Paso_Preconditioner_AMG* Paso_Preconditi Line 251  Paso_Preconditioner_AMG* Paso_Preconditi
251              */              */
252              if ( Esys_noError()) {              if ( Esys_noError()) {
253                 time0=Esys_timer();                 time0=Esys_timer();
254    /*MPI:
255                 out->R=Paso_SystemMatrix_getTranspose(out->P);                 out->R=Paso_SystemMatrix_getTranspose(out->P);
256    */
257                 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);                 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
258              }                    }      
259              /*              /*
# Line 248  Paso_Preconditioner_AMG* Paso_Preconditi Line 261  Paso_Preconditioner_AMG* Paso_Preconditi
261              */              */
262              if ( Esys_noError()) {              if ( Esys_noError()) {
263                 time0=Esys_timer();                 time0=Esys_timer();
264    /*MPI:
265                 Atemp=Paso_SystemMatrix_MatrixMatrix(A_p,out->P);                 Atemp=Paso_SystemMatrix_MatrixMatrix(A_p,out->P);
266                 A_C=Paso_SystemMatrix_MatrixMatrix(out->R,Atemp);                 A_C=Paso_SystemMatrix_MatrixMatrix(out->R,Atemp);
267                  
268                 Paso_SystemMatrix_free(Atemp);                 Paso_SystemMatrix_free(Atemp);
269    */
270    
271                 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);            
272              }              }
273    
# Line 280  Paso_Preconditioner_AMG* Paso_Preconditi Line 297  Paso_Preconditioner_AMG* Paso_Preconditi
297                         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);
298                      #else                      #else
299                         out->A_C=A_C;                         out->A_C=A_C;
300                         out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);                         out->A_C->solver_p=Paso_Preconditioner_Smoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);
301                         out->A_C->solver_package = PASO_SMOOTHER;                         out->A_C->solver_package = PASO_SMOOTHER;
302                         if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);                         if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
303                      #endif                      #endif
# Line 319  void Paso_Preconditioner_AMG_solve(Paso_ Line 336  void Paso_Preconditioner_AMG_solve(Paso_
336    
337       /* presmoothing */       /* presmoothing */
338       time0=Esys_timer();       time0=Esys_timer();
339       Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);       Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
340       time0=Esys_timer()-time0;       time0=Esys_timer()-time0;
341       if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);       if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);
342       /* end of presmoothing */       /* end of presmoothing */
# Line 342  void Paso_Preconditioner_AMG_solve(Paso_ Line 359  void Paso_Preconditioner_AMG_solve(Paso_
359            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);
360            break;            break;
361             case (PASO_SMOOTHER):             case (PASO_SMOOTHER):
362            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_Smoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);
363            break;            break;
364          }          }
365          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 356  void Paso_Preconditioner_AMG_solve(Paso_ Line 373  void Paso_Preconditioner_AMG_solve(Paso_
373                
374          /*solve Ax=b with initial guess x */          /*solve Ax=b with initial guess x */
375          time0=Esys_timer();          time0=Esys_timer();
376          Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);          Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
377          time0=Esys_timer()-time0;          time0=Esys_timer()-time0;
378          if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);          if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
379          /*end of postsmoothing*/          /*end of postsmoothing*/
# Line 696  void Paso_Preconditioner_AMG_enforceFFCo Line 713  void Paso_Preconditioner_AMG_enforceFFCo
713              }              }
714         }         }
715  }  }
716    
717    void Paso_Preconditioner_AMG_CIJPCoarsening( )
718    {
719      const dim_t my_n;
720      const dim_t overlap_n;
721      const dim_t n= my_n + overlap_n;
722       /* set local lambda + overlap */
723       #pragma omp parallel for private(i)
724       for (i=0; i<n ++i) {
725           w[i]=degree_ST[i];
726       }
727       for (i=0; i<my_n; i++) {
728          w2[i]=random;
729       }
730    
731    
732    Paso_Coupler_add(w, 1., w2,
733       /* adjust w accross processors */
734      
735       /*  */
736       global_n_C=0;
737       global_n_F=..;
738      
739       while (global_n_C + global_n_F < global_n) {
740    
741          
742          is_in_D[i]=FALSE;
743          /*  test  local connectivit*/
744          for i ..
745             for k in S_i:
746                 w2[i]=MAX(w2[i],w[k]);
747                 w2[k]=MAX(w2[k],w[i]);
748          /* adjust overlaps by MAX */
749          ...
750          /* points with w[i]>w2[i] become C nodes */
751          
752              
753          
754          
755          
756    
757    
758          global_n_C=?
759          global_n_F=?
760    
761       }
762      
763    }

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

  ViewVC Help
Powered by ViewVC 1.1.26