/[escript]/trunk/paso/src/Solver_AMG.c
ViewVC logotype

Diff of /trunk/paso/src/Solver_AMG.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2113 by artak, Mon Dec 1 06:27:14 2008 UTC revision 2360 by artak, Thu Apr 2 04:12:12 2009 UTC
# Line 71  to Line 71  to
71     then AMG is applied to S again until S becomes empty     then AMG is applied to S again until S becomes empty
72    
73  */  */
74  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level) {  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level, double couplingParam) {
75    Paso_Solver_AMG* out=NULL;    Paso_Solver_AMG* out=NULL;
76    dim_t n=A_p->numRows;    dim_t n=A_p->numRows;
77    dim_t n_block=A_p->row_block_size;    dim_t n_block=A_p->row_block_size;
# Line 112  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 112  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
112       /* identify independend set of rows/columns */       /* identify independend set of rows/columns */
113       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
114       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
115       /*Paso_Pattern_RS(A_p,mis_marker,0.25);*/       Paso_Pattern_coup(A_p,mis_marker,couplingParam);
      /*Paso_Pattern_Aggregiation(A_p,mis_marker,0.5);*/  
      Paso_Pattern_coup(A_p,mis_marker,0.05);  
116       if (Paso_noError()) {       if (Paso_noError()) {
117          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
118          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
# Line 140  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 138  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
138                }                }
139                                
140                /* Compute row-sum for getting rs(A_FF)^-1*/                /* Compute row-sum for getting rs(A_FF)^-1*/
141                #pragma omp parallel for private(i,iPtr,j) schedule(static)                #pragma omp parallel for private(i,iPtr,j,S) schedule(static)
142                for (i = 0; i < out->n_F; ++i) {                for (i = 0; i < out->n_F; ++i) {
143                  S=0;                  S=0;
144                  for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {                  for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
# Line 148  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 146  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
146                   if (mis_marker[j])                   if (mis_marker[j])
147                       S+=A_p->val[iPtr];                       S+=A_p->val[iPtr];
148                  }                  }
                 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsing process.*/  
                 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */  
                 if (ABS(S)>1.e-10)  
149                     out->inv_A_FF[i]=1./S;                     out->inv_A_FF[i]=1./S;
                 else  
                    out->inv_A_FF[i]=1.;  
150                }                }
151    
152             if( Paso_noError()) {             if( Paso_noError()) {
153                /* if there are no nodes in the coarse level there is no more work to do */                /* if there are no nodes in the coarse level there is no more work to do */
154                out->n_C=n-out->n_F;                out->n_C=n-out->n_F;
155                if (level<3) {                if (level>0 && out->n_F>500) {
156                 /*if (out->n_F>500) {*/                 /*if (out->n_F>500) {*/
157                     out->rows_in_C=MEMALLOC(out->n_C,index_t);                     out->rows_in_C=MEMALLOC(out->n_C,index_t);
158                     out->mask_C=MEMALLOC(n,index_t);                     out->mask_C=MEMALLOC(n,index_t);
# Line 191  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 184  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
184    
185                              schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);                              schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);
186                                                            
187                              fprintf(stderr,"Sparsity of Schure: %dx%d LEN %d Percentage %f\n",schur_withFillIn->pattern->numOutput,schur_withFillIn->pattern->numInput,schur_withFillIn->len,schur_withFillIn->len/(1.*schur_withFillIn->pattern->numOutput*schur_withFillIn->pattern->numInput));                              /*fprintf(stderr,"Sparsity of Schure: %dx%d LEN %d Percentage %f\n",schur_withFillIn->pattern->numOutput,schur_withFillIn->pattern->numInput,schur_withFillIn->len,schur_withFillIn->len/(1.*schur_withFillIn->pattern->numOutput*schur_withFillIn->pattern->numInput));*/
188                                                            
189                              /* copy values over*/                              /* copy values over*/
190                              #pragma omp parallel for private(i,iPtr,j,iPtr_s,index,where_p) schedule(static)                              #pragma omp parallel for private(i,iPtr,j,iPtr_s,index,where_p) schedule(static)
# Line 214  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 207  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
207                                                                    
208                              if (Paso_noError()) {                              if (Paso_noError()) {
209                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
210                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level+1);                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1,couplingParam);
211                                  Paso_SparseMatrix_free(schur);                                  Paso_SparseMatrix_free(schur);
212                              }                              }
213                              /* allocate work arrays for AMG application */                              /* allocate work arrays for AMG application */
# Line 254  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 247  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
247    if (Paso_noError()) {    if (Paso_noError()) {
248        if (verbose) {        if (verbose) {
249           printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);           printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);
250           if (level<3) {           if (level>0 && out->n_F>500) {
251           /*if (out->n_F<500) {*/          /* if (out->n_F<500) {*/
252              printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1);              printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1);
253           } else {           } else {
254              printf("timing: AMG: MIS: %e\n",time2);              printf("timing: AMG: MIS: %e\n",time2);
# Line 299  void Paso_Solver_solveAMG(Paso_Solver_AM Line 292  void Paso_Solver_solveAMG(Paso_Solver_AM
292       double *x0=MEMALLOC(amg->n,double);       double *x0=MEMALLOC(amg->n,double);
293       double time0=0;       double time0=0;
294            
295       if (amg->level==3) {       if (amg->level==0  || amg->n_F<=500) {
296       /*if (amg->n_F<=500) {*/       /*if (amg->n_F<=500) {*/
297        time0=Paso_timer();        time0=Paso_timer();
298                    

Legend:
Removed from v.2113  
changed lines
  Added in v.2360

  ViewVC Help
Powered by ViewVC 1.1.26