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

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

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

revision 2783 by artak, Mon Nov 23 06:00:37 2009 UTC revision 2784 by artak, Thu Nov 26 05:09:14 2009 UTC
# Line 72  Paso_SparseMatrix* Paso_SparseMatrix_get Line 72  Paso_SparseMatrix* Paso_SparseMatrix_get
72      out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE);      out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE);
73            
74      k=0;      k=0;
75    
76      for (i=0;i<out->numRows;++i) {      for (i=0;i<out->numRows;++i) {
77        if (mis_marker[i]) {        if (mis_marker[i]) {
78          wptr=W->pattern->ptr[k];          wptr=W->pattern->ptr[k];
# Line 123  Paso_SparseMatrix* Paso_SparseMatrix_get Line 124  Paso_SparseMatrix* Paso_SparseMatrix_get
124      }      }
125        
126    
     /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/  
127      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
128            for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {            for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {
129               j=P->pattern->index[iptr];               j=P->pattern->index[iptr];
# Line 225  void Paso_SparseMatrix_updateWeights(Pas Line 225  void Paso_SparseMatrix_updateWeights(Pas
225         /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/         /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/
226        }        }
227     }     }
228            #pragma omp parallel for private(i,iPtr,j) schedule(static)
229        for (i = 0; i < W_FC->numRows; ++i) {        for (i = 0; i < W_FC->numRows; ++i) {
230              for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {              for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
231                   j=W_FC->pattern->index[iPtr];                   j=W_FC->pattern->index[iPtr];
# Line 246  void Paso_SparseMatrix_updateWeights(Pas Line 246  void Paso_SparseMatrix_updateWeights(Pas
246        TMPMEMFREE(beta_den);        TMPMEMFREE(beta_den);
247  }  }
248    
 /*A_c=R*A*P taking into account coarseneing */  
 /*  
 void Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A_c, Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {  
   
   index_t iptrA_c,iptrR,iptrP,iptrA;  
   dim_t i,j,k,l,m;  
   double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;  
  */  
   /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/  
  /*  
   for(i = 0; i < A_c->numRows; i++) {  
      for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {  
       j = A_c->pattern->index[iptrA_c];  
       second_sum=0;  
       for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {  
         k = R->pattern->index[iptrR];  
         first_sum=0;  
         for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {  
          l= A->pattern->index[iptrA];  
  */        /*This loop has to be replaced by bsearch */  
   /*       for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {  
            m=P->pattern->index[iptrP];  
            if(m==j) {  
             p_lj=P->val[iptrP];  
             break;  
            }  
          }  
          a_kl=A->val[iptrA];  
          first_sum+=a_kl*p_lj;  
         }  
         r_ik=R->val[iptrR];  
         second_sum+=r_ik*first_sum;  
         if(i==0) {  
         printf("row %d, col %d, k=%d, firstsum = %e, r_ik=%e \n",i,j,k,first_sum,r_ik);  
         }  
       }  
       if(i==0) {  
         printf("row %d, col %d, secondsum = %e \n",i,j,second_sum);  
       }  
       A_c->val[iptrA_c]=second_sum;  
      }  
     }  
   
 }  
   */  
249        
250  /*A_c=R*A*P taking into account coarseneing */  /*A_c=R*A*P taking into account coarseneing */
251  /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {  /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {
# Line 314  Paso_SparseMatrix* Paso_SparseMatrix_Mat Line 269  Paso_SparseMatrix* Paso_SparseMatrix_Mat
269    Paso_SparseMatrix *out=NULL;    Paso_SparseMatrix *out=NULL;
270    double sum,b_lj=0;    double sum,b_lj=0;
271        
272      double time0=0;
273      bool_t verbose=0;
274      
275      time0=Paso_timer();
276      
277    outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);    outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
278    out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);    out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
279        
280      time0=Paso_timer()-time0;
281      if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0);
282        
283      time0=Paso_timer();
284      
285      #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static)
286    for(i = 0; i < out->numRows; i++) {    for(i = 0; i < out->numRows; i++) {
287      for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {      for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
288        j = out->pattern->index[iptrC];        j = out->pattern->index[iptrC];
289        sum=0;        sum=0;
       b_lj=0;  
290        for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {        for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
291              k=A->pattern->index[iptrA];              k=A->pattern->index[iptrA];
292              /*can be replaced by bsearch */              /*can be replaced by bsearch */
293                b_lj=0;
294              for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {              for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
295                 l=B->pattern->index[iptrB];                 l=B->pattern->index[iptrB];
296                 if(l==j) {                 if(l==j) {
# Line 338  Paso_SparseMatrix* Paso_SparseMatrix_Mat Line 303  Paso_SparseMatrix* Paso_SparseMatrix_Mat
303        out->val[iptrC]=sum;        out->val[iptrC]=sum;
304     }     }
305    }    }
306      
307      time0=Paso_timer()-time0;
308      if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0);
309                
310    Paso_Pattern_free(outpattern);    Paso_Pattern_free(outpattern);
311    
# Line 354  Paso_SparseMatrix* Paso_Solver_getCoarse Line 322  Paso_SparseMatrix* Paso_Solver_getCoarse
322    Paso_Pattern* outpattern=NULL;    Paso_Pattern* outpattern=NULL;
323    Paso_SparseMatrix *A_c=NULL;    Paso_SparseMatrix *A_c=NULL;
324        
325      double time0=0;
326      bool_t verbose=0;
327      
328      time0=Paso_timer();
329      
330    temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);    temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
331    outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);    outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
332    A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);    A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
333        
334      time0=Paso_timer()-time0;
335      if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0);
336        
337    /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/    /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
338    
339      time0=Paso_timer();
340      
341      #pragma omp parallel for private(i,iptrA_c,j,second_sum,iptrR,k,first_sum,p_lj,iptrP,m,a_kl,r_ik) schedule(static)
342    for(i = 0; i < A_c->numRows; i++) {    for(i = 0; i < A_c->numRows; i++) {
343       for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {       for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
344        j = A_c->pattern->index[iptrA_c];        j = A_c->pattern->index[iptrA_c];
# Line 387  Paso_SparseMatrix* Paso_Solver_getCoarse Line 365  Paso_SparseMatrix* Paso_Solver_getCoarse
365        }        }
366        A_c->val[iptrA_c]=second_sum;        A_c->val[iptrA_c]=second_sum;
367       }       }
368      }    }
369      time0=Paso_timer()-time0;
370      if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0);
371            
372      Paso_Pattern_free(outpattern);    Paso_Pattern_free(outpattern);
373      Paso_Pattern_free(temp);    Paso_Pattern_free(temp);
374            
375  return A_c;  return A_c;
376  }  }

Legend:
Removed from v.2783  
changed lines
  Added in v.2784

  ViewVC Help
Powered by ViewVC 1.1.26