/[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 2704 by artak, Thu Oct 1 05:38:33 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 57  void Paso_Pattern_coup(Paso_SparseMatrix Line 56  void Paso_Pattern_coup(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 79  void Paso_Pattern_coup(Paso_SparseMatrix Line 78  void Paso_Pattern_coup(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 93  void Paso_Pattern_coup(Paso_SparseMatrix Line 92  void Paso_Pattern_coup(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;
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 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_REMOVED;
                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 145  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 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,max_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;          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);
181              }              }
182          }          }
183         }         }
184        }        }
185            
186      
187      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);
188            
189       /* clean up */       /* clean up */
# Line 198  void Paso_Pattern_RS(Paso_SparseMatrix* Line 193  void Paso_Pattern_RS(Paso_SparseMatrix*
193       }       }
194      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
195    
196      Paso_Pattern_mis(out,mis_marker);      /*Paso_Pattern_mis(out,mis_marker);*/
197        Paso_Pattern_greedy(out,mis_marker);
198        Paso_Pattern_free(out);
199  }  }
200    
201  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 222  void Paso_Pattern_Aggregiation(Paso_Spar
222      }      }
223            
224    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
225      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");
226      return;      return;
227    }    }
228    
# Line 234  void Paso_Pattern_Aggregiation(Paso_Spar Line 231  void Paso_Pattern_Aggregiation(Paso_Spar
231        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
232          diag = 0;          diag = 0;
233          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) {
234              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] == i){
235                  diag+=A->val[iptr];                  diag+=A->val[iptr];
236              }              }
237          }          }
# Line 246  void Paso_Pattern_Aggregiation(Paso_Spar Line 243  void Paso_Pattern_Aggregiation(Paso_Spar
243       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
244         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
245          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
         val=0.;  
246          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) {
247              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
248              val=A->val[iptr];              val=A->val[iptr];
249              if(j!= i) {                if((val*val)>=(eps_Aii*diags[j])) {
               if(val*val>=eps_Aii * diags[j]) {  
250                 Paso_IndexList_insertIndex(&(index_list[i]),j);                 Paso_IndexList_insertIndex(&(index_list[i]),j);
251                }                }
             }  
252          }          }
253         }         }
254       }       }
# Line 270  void Paso_Pattern_Aggregiation(Paso_Spar Line 264  void Paso_Pattern_Aggregiation(Paso_Spar
264      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
265      MEMFREE(diags);      MEMFREE(diags);
266            
267      Paso_Pattern_mis(out,mis_marker);      
268        /*Paso_Pattern_mis(out,mis_marker);*/
269        Paso_Pattern_greedy(out,mis_marker);
270        Paso_Pattern_free(out);
271    
272    }
273    
274    /* Greedy algorithm */
275    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
276    
277      dim_t i,j;
278      /*double sum;*/
279      index_t iptr;
280      bool_t passed=FALSE;
281      dim_t n=pattern->numOutput;
282    
283      if (pattern->type & PATTERN_FORMAT_SYM) {
284        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
285        return;
286      }
287      
288       /* We do not need this loop if we set IS_REMOVED=IS_AVAILABLE. */
289       #pragma omp parallel for private(i) schedule(static)
290       for (i=0;i<n;++i)
291            if(mis_marker[i]==IS_AVAILABLE)
292                        mis_marker[i]=IS_IN_SET;
293    
294    
295        for (i=0;i<n;++i) {
296          if (mis_marker[i]==IS_IN_SET) {
297            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
298                 j=pattern->index[iptr];
299                 mis_marker[j]=IS_REMOVED;
300            }
301          }
302        }
303        
304        
305        
306        for (i=0;i<n;i++) {
307            if (mis_marker[i]==IS_REMOVED) {
308               passed=TRUE;
309               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
310                  j=pattern->index[iptr];
311                    if (mis_marker[j]==IS_REMOVED) {
312                        passed=TRUE;
313                    }
314                    else {
315                        passed=FALSE;
316                        break;
317                    }
318                  }
319               if (passed) mis_marker[i]=IS_IN_SET;
320               }
321            }
322    
323         /* swap to TRUE/FALSE in mis_marker */
324         #pragma omp parallel for private(i) schedule(static)
325         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_REMOVED);
326        
327  }  }
328    
329    
330    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
331    
332      dim_t i,j;
333      /*double sum;*/
334      index_t iptr;
335      index_t num_colors;
336      index_t* colorOf;
337      register index_t color;
338      bool_t passed=FALSE;
339      dim_t n=pattern->numOutput;
340    
341      
342      colorOf=MEMALLOC(n,index_t);
343    
344      if (pattern->type & PATTERN_FORMAT_SYM) {
345        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
346        return;
347      }
348      
349       Paso_Pattern_color(pattern,&num_colors,colorOf);
350      
351       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
352       #pragma omp parallel for private(i) schedule(static)
353       for (i=0;i<n;++i)
354            if(mis_marker[i]==IS_AVAILABLE)
355                        mis_marker[i]=IS_IN_SET;
356    
357       #pragma omp barrier
358       for (color=0;color<num_colors;++color) {
359        #pragma omp parallel for schedule(static) private(i,iptr,j)
360        for (i=0;i<n;++i) {
361         if (colorOf[i]==color) {  
362          if (mis_marker[i]==IS_IN_SET) {
363            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
364                 j=pattern->index[iptr];
365                 if (colorOf[j]<color)
366                  mis_marker[j]=IS_REMOVED;
367            }
368          }
369         }
370        }
371       }
372        
373        
374       #pragma omp barrier
375       for (color=0;color<num_colors;++color) {
376       #pragma omp parallel for schedule(static) private(i,iptr,j)
377        for (i=0;i<n;i++) {
378          if (colorOf[i]==color) {  
379            if (mis_marker[i]==IS_REMOVED) {
380               passed=TRUE;
381               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
382                  j=pattern->index[iptr];
383                   if (colorOf[j]<color && passed) {
384                    if (mis_marker[j]==IS_REMOVED) {
385                        passed=TRUE;
386                    }
387                    else {
388                        passed=FALSE;
389                        /*break;*/
390                    }
391                  }
392               }
393               if (passed) mis_marker[i]=IS_IN_SET;
394               }
395              
396            }
397        }
398       }
399    
400         /* swap to TRUE/FALSE in mis_marker */
401         #pragma omp parallel for private(i) schedule(static)
402         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
403        
404        MEMFREE(colorOf);
405    }
406    
407  #undef IS_AVAILABLE  #undef IS_AVAILABLE
408  #undef IS_IN_SET  #undef IS_IN_SET
409  #undef IS_REMOVED  #undef IS_REMOVED
410    
411    
412    
413    
414    
415    
416    
417    
418    
419    

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

  ViewVC Help
Powered by ViewVC 1.1.26