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

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

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

revision 2450 by artak, Tue Jun 2 01:22:11 2009 UTC revision 2475 by artak, Wed Jun 17 01:48:46 2009 UTC
# Line 162  void Paso_Pattern_RS(Paso_SparseMatrix* Line 162  void Paso_Pattern_RS(Paso_SparseMatrix*
162      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
163      return;      return;
164    }    }
165        #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold,j) schedule(static)
     #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)  
166      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
167        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
168          min_offdiagonal = DBL_MAX;          min_offdiagonal = DBL_MAX;
# Line 172  void Paso_Pattern_RS(Paso_SparseMatrix* Line 171  void Paso_Pattern_RS(Paso_SparseMatrix*
171                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
172              }              }
173          }          }
174            
175          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
176          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
177              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
# Line 186  void Paso_Pattern_RS(Paso_SparseMatrix* Line 185  void Paso_Pattern_RS(Paso_SparseMatrix*
185         }         }
186        }        }
187            
188      
189      out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);      out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
190            
191       /* clean up */       /* clean up */
# Line 268  void Paso_Pattern_Aggregiation(Paso_Spar Line 268  void Paso_Pattern_Aggregiation(Paso_Spar
268      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
269      MEMFREE(diags);      MEMFREE(diags);
270            
271        
272      /*Paso_Pattern_mis(out,mis_marker);*/      /*Paso_Pattern_mis(out,mis_marker);*/
273      Paso_Pattern_greedy(out,mis_marker);      Paso_Pattern_greedy(out,mis_marker);
274    
# Line 287  void Paso_Pattern_greedy(Paso_Pattern* p Line 288  void Paso_Pattern_greedy(Paso_Pattern* p
288      return;      return;
289    }    }
290        
291       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
292     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
293     for (i=0;i<n;++i)     for (i=0;i<n;++i)
294          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
# Line 328  void Paso_Pattern_greedy(Paso_Pattern* p Line 330  void Paso_Pattern_greedy(Paso_Pattern* p
330  }  }
331    
332    
333    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
334    
335      dim_t i,j;
336      /*double sum;*/
337      index_t iptr;
338      index_t num_colors;
339      index_t* colorOf;
340      register index_t color;
341      bool_t passed=FALSE;
342      dim_t n=pattern->numOutput;
343    
344      
345      colorOf=MEMALLOC(n,index_t);
346    
347      if (pattern->type & PATTERN_FORMAT_SYM) {
348        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
349        return;
350      }
351      
352       Paso_Pattern_color(pattern,&num_colors,colorOf);
353      
354       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
355       #pragma omp parallel for private(i) schedule(static)
356       for (i=0;i<n;++i)
357            if(mis_marker[i]==IS_AVAILABLE)
358                        mis_marker[i]=IS_IN_SET;
359    
360       #pragma omp barrier
361       for (color=0;color<num_colors;++color) {
362        #pragma omp parallel for schedule(static) private(i,iptr,j)
363        for (i=0;i<n;++i) {
364         if (colorOf[i]==color) {  
365          if (mis_marker[i]==IS_IN_SET) {
366            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
367                 j=pattern->index[iptr];
368                 if (colorOf[j]<color)
369                  mis_marker[j]=IS_REMOVED;
370            }
371          }
372         }
373        }
374       }
375        
376        
377       #pragma omp barrier
378       for (color=0;color<num_colors;++color) {
379       #pragma omp parallel for schedule(static) private(i,iptr,j)
380        for (i=0;i<n;i++) {
381          if (colorOf[i]==color) {  
382            if (mis_marker[i]==IS_REMOVED) {
383               passed=TRUE;
384               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
385                  j=pattern->index[iptr];
386                   if (colorOf[j]<color && passed) {
387                    if (mis_marker[j]==IS_REMOVED) {
388                        passed=TRUE;
389                    }
390                    else {
391                        passed=FALSE;
392                        /*break;*/
393                    }
394                  }
395               }
396               }
397               if (passed) mis_marker[i]=IS_IN_SET;
398            }
399        }
400       }
401    
402         /* swap to TRUE/FALSE in mis_marker */
403         #pragma omp parallel for private(i) schedule(static)
404         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
405        
406        MEMFREE(colorOf);
407    }
408    
409    
410    
411  #undef IS_AVAILABLE  #undef IS_AVAILABLE

Legend:
Removed from v.2450  
changed lines
  Added in v.2475

  ViewVC Help
Powered by ViewVC 1.1.26