/[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 2507 by artak, Thu Jul 2 03:14:35 2009 UTC revision 2509 by artak, Thu Jul 2 04:51:51 2009 UTC
# Line 85  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 85  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
85    Paso_SparseMatrix * schur=NULL;    Paso_SparseMatrix * schur=NULL;
86    Paso_SparseMatrix * schur_withFillIn=NULL;    Paso_SparseMatrix * schur_withFillIn=NULL;
87    double S=0;    double S=0;
   double time0;  
88        
89    /*Make sure we have block sizes 1*/    /*Make sure we have block sizes 1*/
90    A_p->pattern->input_block_size=A_p->col_block_size;    A_p->pattern->input_block_size=A_p->col_block_size;
# Line 121  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 120  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
120       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
121       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
122            
      time0=Paso_timer();  
   
123       if (coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {       if (coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {
124        Paso_Pattern_coup(A_p,mis_marker,coarsening_threshold);        Paso_Pattern_coup(A_p,mis_marker,coarsening_threshold);
125       }       }
# Line 136  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 133  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
133        /*Default coarseneing*/        /*Default coarseneing*/
134        Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);        Paso_Pattern_RS(A_p,mis_marker,coarsening_threshold);
135       }       }
136        
      time0=Paso_timer()-time0;  
      if (verbose) fprintf(stderr,"timing: Coarseneing: %e\n",time0);  
137      if (Paso_noError()) {      if (Paso_noError()) {
138          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
139          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
# Line 208  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 203  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
203                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
204                              /*find the pattern of the schur complement with fill in*/                              /*find the pattern of the schur complement with fill in*/
205            
                             time0=Paso_timer();  
206                              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);
                             time0=Paso_timer()-time0;  
                             if (verbose) fprintf(stderr,"timing: Fill_in pattern computation: %e\n",time0);  
       
                             time0=Paso_timer();  
207                              /* copy values over*/                              /* copy values over*/
208                              #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)
209                              for (i = 0; i < schur_withFillIn->numRows; ++i) {                              for (i = 0; i < schur_withFillIn->numRows; ++i) {
# Line 231  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 221  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
221                                      schur_withFillIn->val[iPtr]=schur->val[iPtr_s+(index_t)(where_p-index)];                                      schur_withFillIn->val[iPtr]=schur->val[iPtr_s+(index_t)(where_p-index)];
222                                  }                                  }
223                                }                                }
224                              }                              }                                
                             time0=Paso_timer()-time0;  
                             if (verbose) fprintf(stderr,"timing: Copying values for Fill in: %e\n",time0);  
                                   
225                              if (Paso_noError()) {                              if (Paso_noError()) {
226                                  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);
227                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1,coarsening_threshold,coarsening_method);                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1,coarsening_threshold,coarsening_method);
# Line 314  void Paso_Solver_solveAMG(Paso_Solver_AM Line 301  void Paso_Solver_solveAMG(Paso_Solver_AM
301       double *r=MEMALLOC(amg->n,double);       double *r=MEMALLOC(amg->n,double);
302       double *x0=MEMALLOC(amg->n,double);       double *x0=MEMALLOC(amg->n,double);
303       double time0=0;       double time0=0;
304       bool_t verbose=1;       bool_t verbose=0;
305            
306       #ifdef MKL       #ifdef MKL
307       Paso_SparseMatrix *temp=NULL;       Paso_SparseMatrix *temp=NULL;
# Line 323  void Paso_Solver_solveAMG(Paso_Solver_AM Line 310  void Paso_Solver_solveAMG(Paso_Solver_AM
310            
311       if (amg->level==0  || amg->n_C<=500) {       if (amg->level==0  || amg->n_C<=500) {
312                
      /*if (amg->n_F<=500) {*/  
313        time0=Paso_timer();        time0=Paso_timer();
         
           
        /* Paso_Solver_solveJacobi(amg->GS,x,b);*/  
           
         /* #pragma omp parallel for private(i) schedule(static)  
         for (i=0;i<amg->n;++i) r[i]=b[i];  
         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);  
         Paso_Solver_solveGS(amg->GS,x0,r);  
         #pragma omp parallel for private(i) schedule(static)  
         for (i=0;i<amg->n;++i) {  
          x[i]+=x0[i];  
         }  
         */  
         
314         #ifdef MKL         #ifdef MKL
         /*dim_t iptr;  
         temp=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, amg->A->pattern,1,1);  
         #pragma omp parallel for private(i,iptr) schedule(static)  
         for (i=0;i<amg->A->numRows;++i) {  
          for (iptr=amg->A->pattern->ptr[i];iptr<amg->A->pattern->ptr[i+1]; ++iptr)  
              temp->val[iptr]=amg->A->val[iptr];  
         }*/  
315          temp=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, amg->A->pattern,1,1);          temp=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, amg->A->pattern,1,1);
316          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
317          for (i=0;i<amg->A->len;++i) {          for (i=0;i<amg->A->len;++i) {
# Line 371  void Paso_Solver_solveAMG(Paso_Solver_AM Line 336  void Paso_Solver_solveAMG(Paso_Solver_AM
336           Paso_Solver_solveJacobi(amg->GS,x,b);           Paso_Solver_solveJacobi(amg->GS,x,b);
337           time0=Paso_timer()-time0;           time0=Paso_timer()-time0;
338           if (verbose) fprintf(stderr,"timing: Presmooting: %e\n",time0);           if (verbose) fprintf(stderr,"timing: Presmooting: %e\n",time0);
          /* #pragma omp parallel for private(i) schedule(static)  
          for (i=0;i<amg->n;++i) r[i]=b[i];  
   
           Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);  
           Paso_Solver_solveJacobi(amg->GS,x0,r);  
             
           #pragma omp parallel for private(i) schedule(static)  
           for (i=0;i<amg->n;++i) {  
            x[i]+=x0[i];  
           }  
          */  
339          /* end of presmoothing */          /* end of presmoothing */
340                    
341                    
# Line 444  void Paso_Solver_solveAMG(Paso_Solver_AM Line 398  void Paso_Solver_solveAMG(Paso_Solver_AM
398            
399       time0=Paso_timer()-time0;       time0=Paso_timer()-time0;
400       if (verbose) fprintf(stderr,"timing: Postsmoothing: %e\n",time0);       if (verbose) fprintf(stderr,"timing: Postsmoothing: %e\n",time0);
401       /*  
      #pragma omp parallel for private(i) schedule(static)  
      for (i=0;i<amg->n;++i) r[i]=b[i];  
       
      Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);  
      Paso_Solver_solveJacobi(amg->GS,x0,r);  
       
      #pragma omp parallel for private(i) schedule(static)  
      for (i=0;i<amg->n;++i) {  
       x[i]+=x0[i];  
      }  
      */  
402       /*end of postsmoothing*/       /*end of postsmoothing*/
403            
404       }       }

Legend:
Removed from v.2507  
changed lines
  Added in v.2509

  ViewVC Help
Powered by ViewVC 1.1.26