/[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 2381 by artak, Tue Apr 14 03:46:59 2009 UTC revision 2652 by artak, Mon Sep 7 05:04:45 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 36  Line 36 
36  /***************************************************************/  /***************************************************************/
37    
38  #define IS_AVAILABLE -1  #define IS_AVAILABLE -1
39  #define IS_IN_SET -3  #define IS_IN_SET -3   /* Week connection */
40  #define IS_REMOVED -4  #define IS_REMOVED -4  /* strong */
41    
42  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43    
44    dim_t i,j;    dim_t i,j;
45    /*double threshold=0.05;*/    /*double sum;*/
   double sum;  
46    index_t iptr,*index,*where_p,*diagptr;    index_t iptr,*index,*where_p,*diagptr;
47    bool_t passed=FALSE;    bool_t passed=FALSE;
48    dim_t n=A->numRows;    dim_t n=A->numRows;
# Line 94  void Paso_Pattern_coup(Paso_SparseMatrix Line 93  void Paso_Pattern_coup(Paso_SparseMatrix
93        /*This loop cannot be parallelized, as order matters here.*/        /*This loop cannot be parallelized, as order matters here.*/
94      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
95          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_REMOVED) {
96               passed=TRUE;
97             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) {
98                j=A->pattern->index[iptr];                j=A->pattern->index[iptr];
99                if (mis_marker[j]==IS_IN_SET) {                if (mis_marker[j]==IS_IN_SET) {
# Line 105  void Paso_Pattern_coup(Paso_SparseMatrix Line 105  void Paso_Pattern_coup(Paso_SparseMatrix
105                      break;                      break;
106                  }                  }
107                }                }
108              }             }
109              if (passed) {             if (passed) mis_marker[i]=IS_IN_SET;
                mis_marker[i]=IS_IN_SET;  
                passed=FALSE;  
             }  
110          }          }
111      }      }
   
112      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113      /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */      /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114      #pragma omp parallel for private(i,iptr,j,sum) schedule(static)      /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
116          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_REMOVED) {
117             sum=0;             sum=0;
# Line 128  void Paso_Pattern_coup(Paso_SparseMatrix Line 124  void Paso_Pattern_coup(Paso_SparseMatrix
124               mis_marker[i]=IS_IN_SET;               mis_marker[i]=IS_IN_SET;
125          }          }
126      }      }
127        */
128    
129       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
130       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
# Line 166  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 176  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];
178              if(A->val[iptr]<=threshold) {              if(A->val[iptr]<=threshold) {
179                 if(j!=i) {                 if(j!=i) {
180                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
181                    /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
182                  }                  }
183              }              }
184          }          }
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 198  void Paso_Pattern_RS(Paso_SparseMatrix* Line 195  void Paso_Pattern_RS(Paso_SparseMatrix*
195       }       }
196      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
197    
198      Paso_Pattern_mis(out,mis_marker);      /*Paso_Pattern_mis(out,mis_marker);*/
199        Paso_Pattern_greedy(out,mis_marker);
200        Paso_Pattern_free(out);
201  }  }
202    
203  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
# Line 225  void Paso_Pattern_Aggregiation(Paso_Spar Line 224  void Paso_Pattern_Aggregiation(Paso_Spar
224      }      }
225            
226    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
227      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
228      return;      return;
229    }    }
230    
# Line 234  void Paso_Pattern_Aggregiation(Paso_Spar Line 233  void Paso_Pattern_Aggregiation(Paso_Spar
233        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
234          diag = 0;          diag = 0;
235          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) {
236              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] == i){
237                  diag+=A->val[iptr];                  diag+=A->val[iptr];
238              }              }
239          }          }
# Line 270  void Paso_Pattern_Aggregiation(Paso_Spar Line 269  void Paso_Pattern_Aggregiation(Paso_Spar
269      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
270      MEMFREE(diags);      MEMFREE(diags);
271            
272      Paso_Pattern_mis(out,mis_marker);      
273        /*Paso_Pattern_mis(out,mis_marker);*/
274        Paso_Pattern_greedy(out,mis_marker);
275        Paso_Pattern_free(out);
276    
277    }
278    
279    /* Greedy algorithm */
280    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
281    
282      dim_t i,j;
283      /*double sum;*/
284      index_t iptr;
285      bool_t passed=FALSE;
286      dim_t n=pattern->numOutput;
287    
288      if (pattern->type & PATTERN_FORMAT_SYM) {
289        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
290        return;
291      }
292      
293       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
294       #pragma omp parallel for private(i) schedule(static)
295       for (i=0;i<n;++i)
296            if(mis_marker[i]==IS_AVAILABLE)
297                        mis_marker[i]=IS_IN_SET;
298    
299    
300        for (i=0;i<n;++i) {
301          if (mis_marker[i]==IS_IN_SET) {
302            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
303                 j=pattern->index[iptr];
304                 mis_marker[j]=IS_REMOVED;
305            }
306          }
307        }
308        
309        
310        
311        for (i=0;i<n;i++) {
312            if (mis_marker[i]==IS_REMOVED) {
313               passed=TRUE;
314               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
315                  j=pattern->index[iptr];
316                    if (mis_marker[j]==IS_REMOVED) {
317                        passed=TRUE;
318                    }
319                    else {
320                        passed=FALSE;
321                        break;
322                    }
323                  }
324               }
325               if (passed) mis_marker[i]=IS_IN_SET;
326            }
327    
328         /* swap to TRUE/FALSE in mis_marker */
329         #pragma omp parallel for private(i) schedule(static)
330         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
331        
332  }  }
333    
334    
335    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
336    
337      dim_t i,j;
338      /*double sum;*/
339      index_t iptr;
340      index_t num_colors;
341      index_t* colorOf;
342      register index_t color;
343      bool_t passed=FALSE;
344      dim_t n=pattern->numOutput;
345    
346      
347      colorOf=MEMALLOC(n,index_t);
348    
349      if (pattern->type & PATTERN_FORMAT_SYM) {
350        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
351        return;
352      }
353      
354       Paso_Pattern_color(pattern,&num_colors,colorOf);
355      
356       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
357       #pragma omp parallel for private(i) schedule(static)
358       for (i=0;i<n;++i)
359            if(mis_marker[i]==IS_AVAILABLE)
360                        mis_marker[i]=IS_IN_SET;
361    
362       #pragma omp barrier
363       for (color=0;color<num_colors;++color) {
364        #pragma omp parallel for schedule(static) private(i,iptr,j)
365        for (i=0;i<n;++i) {
366         if (colorOf[i]==color) {  
367          if (mis_marker[i]==IS_IN_SET) {
368            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
369                 j=pattern->index[iptr];
370                 if (colorOf[j]<color)
371                  mis_marker[j]=IS_REMOVED;
372            }
373          }
374         }
375        }
376       }
377        
378        
379       #pragma omp barrier
380       for (color=0;color<num_colors;++color) {
381       #pragma omp parallel for schedule(static) private(i,iptr,j)
382        for (i=0;i<n;i++) {
383          if (colorOf[i]==color) {  
384            if (mis_marker[i]==IS_REMOVED) {
385               passed=TRUE;
386               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
387                  j=pattern->index[iptr];
388                   if (colorOf[j]<color && passed) {
389                    if (mis_marker[j]==IS_REMOVED) {
390                        passed=TRUE;
391                    }
392                    else {
393                        passed=FALSE;
394                        /*break;*/
395                    }
396                  }
397               }
398               }
399               if (passed) mis_marker[i]=IS_IN_SET;
400            }
401        }
402       }
403    
404         /* swap to TRUE/FALSE in mis_marker */
405         #pragma omp parallel for private(i) schedule(static)
406         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
407        
408        MEMFREE(colorOf);
409    }
410    
411  #undef IS_AVAILABLE  #undef IS_AVAILABLE
412  #undef IS_IN_SET  #undef IS_IN_SET

Legend:
Removed from v.2381  
changed lines
  Added in v.2652

  ViewVC Help
Powered by ViewVC 1.1.26