/[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 2685 by artak, Tue Sep 15 03:05:23 2009 UTC revision 2686 by artak, Tue Sep 29 03:39:36 2009 UTC
# Line 56  void Paso_Pattern_YS(Paso_SparseMatrix* Line 56  void Paso_Pattern_YS(Paso_SparseMatrix*
56     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
57     for (i=0;i<n;++i)     for (i=0;i<n;++i)
58          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
59                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_REMOVED;
60    
61      #pragma omp parallel for private(i,index,where_p) schedule(static)      #pragma omp parallel for private(i,index,where_p) schedule(static)
62      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
63           diagptr[i]=A->pattern->ptr[i];           diagptr[i]=A->pattern->ptr[i];
64           index=&(A->pattern->index[diagptr[i]]);           index=&(A->pattern->index[A->pattern->ptr[i]]);
65           where_p=(index_t*)bsearch(&i,           where_p=(index_t*)bsearch(&i,
66                                  index,                                  index,
67                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],
# Line 78  void Paso_Pattern_YS(Paso_SparseMatrix* Line 78  void Paso_Pattern_YS(Paso_SparseMatrix*
78    
79      /*This loop cannot be parallelized, as order matters here.*/      /*This loop cannot be parallelized, as order matters here.*/
80      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
81        if (mis_marker[i]==IS_IN_SET) {        if (mis_marker[i]==IS_REMOVED) {
82          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) {
83               j=A->pattern->index[iptr];               j=A->pattern->index[iptr];
84               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85                  mis_marker[j]=IS_REMOVED;                  mis_marker[j]=IS_IN_SET;
86               }               }
87          }          }
88        }        }
# Line 92  void Paso_Pattern_YS(Paso_SparseMatrix* Line 92  void Paso_Pattern_YS(Paso_SparseMatrix*
92            
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_IN_SET) {
96             passed=TRUE;             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_REMOVED) {
100                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101                      passed=TRUE;                      passed=TRUE;
102                  }                  }
# Line 106  void Paso_Pattern_YS(Paso_SparseMatrix* Line 106  void Paso_Pattern_YS(Paso_SparseMatrix*
106                  }                  }
107                }                }
108             }             }
109             if (passed) mis_marker[i]=IS_IN_SET;             if (passed) mis_marker[i]=IS_REMOVED;
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.*/
# Line 141  void Paso_Pattern_RS(Paso_SparseMatrix* Line 141  void Paso_Pattern_RS(Paso_SparseMatrix*
141  {  {
142    dim_t i,n,j;    dim_t i,n,j;
143    index_t iptr;    index_t iptr;
144    double threshold,min_offdiagonal;    double threshold,max_offdiagonal;
145        
146    Paso_Pattern *out=NULL;    Paso_Pattern *out=NULL;
147        
# 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,max_offdiagonal,threshold,j) 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;          max_offdiagonal = DBL_MIN;
169          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) {
170              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
171                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);                  max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
172              }              }
173          }          }
174                    
175          threshold = theta*min_offdiagonal;          threshold = theta*max_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) {
                if(j!=i) {  
179                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
180                  /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/                  Paso_IndexList_insertIndex(&(index_list[j]),i);
                 }  
181              }              }
182          }          }
183         }         }
# Line 293  void Paso_Pattern_greedy(Paso_Pattern* p Line 291  void Paso_Pattern_greedy(Paso_Pattern* p
291     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
292     for (i=0;i<n;++i)     for (i=0;i<n;++i)
293          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
294                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_REMOVED;
295    
296    
297      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
298        if (mis_marker[i]==IS_IN_SET) {        if (mis_marker[i]==IS_REMOVED) {
299          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
300               j=pattern->index[iptr];               j=pattern->index[iptr];
301               mis_marker[j]=IS_REMOVED;               mis_marker[j]=IS_IN_SET;
302          }          }
303        }        }
304      }      }
# Line 308  void Paso_Pattern_greedy(Paso_Pattern* p Line 306  void Paso_Pattern_greedy(Paso_Pattern* p
306            
307            
308      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
309          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_IN_SET) {
310             passed=TRUE;             passed=TRUE;
311             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
312                j=pattern->index[iptr];                j=pattern->index[iptr];
313                  if (mis_marker[j]==IS_REMOVED) {                  if (mis_marker[j]==IS_IN_SET) {
314                      passed=TRUE;                      passed=TRUE;
315                  }                  }
316                  else {                  else {
# Line 321  void Paso_Pattern_greedy(Paso_Pattern* p Line 319  void Paso_Pattern_greedy(Paso_Pattern* p
319                  }                  }
320                }                }
321             }             }
322             if (passed) mis_marker[i]=IS_IN_SET;             if (passed) mis_marker[i]=IS_REMOVED;
323          }          }
324    
325       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
# Line 410  void Paso_Pattern_greedy_color(Paso_Patt Line 408  void Paso_Pattern_greedy_color(Paso_Patt
408  #undef IS_AVAILABLE  #undef IS_AVAILABLE
409  #undef IS_IN_SET  #undef IS_IN_SET
410  #undef IS_REMOVED  #undef IS_REMOVED
411    
412    
413    
414    
415    
416    
417    
418    
419    
420    

Legend:
Removed from v.2685  
changed lines
  Added in v.2686

  ViewVC Help
Powered by ViewVC 1.1.26