/[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 2551 by gross, Thu Jul 23 09:19:15 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_coup(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    
231    
232      #pragma omp parallel for private(i,iptr,diag) schedule(static)      #pragma omp parallel for private(i,iptr) reduction(+:diag) schedule(static)
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    
276    }
277    
278    /* Greedy algorithm */
279    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
280    
281      dim_t i,j;
282      /*double sum;*/
283      index_t iptr;
284      bool_t passed=FALSE;
285      dim_t n=pattern->numOutput;
286    
287      if (pattern->type & PATTERN_FORMAT_SYM) {
288        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
289        return;
290      }
291      
292       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
293       #pragma omp parallel for private(i) schedule(static)
294       for (i=0;i<n;++i)
295            if(mis_marker[i]==IS_AVAILABLE)
296                        mis_marker[i]=IS_IN_SET;
297    
298    
299        for (i=0;i<n;++i) {
300          if (mis_marker[i]==IS_IN_SET) {
301            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
302                 j=pattern->index[iptr];
303                 mis_marker[j]=IS_REMOVED;
304            }
305          }
306        }
307        
308        
309        
310        for (i=0;i<n;i++) {
311            if (mis_marker[i]==IS_REMOVED) {
312               passed=TRUE;
313               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
314                  j=pattern->index[iptr];
315                    if (mis_marker[j]==IS_REMOVED) {
316                        passed=TRUE;
317                    }
318                    else {
319                        passed=FALSE;
320                        break;
321                    }
322                  }
323               }
324               if (passed) mis_marker[i]=IS_IN_SET;
325            }
326    
327         /* swap to TRUE/FALSE in mis_marker */
328         #pragma omp parallel for private(i) schedule(static)
329         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
330        
331  }  }
332    
333    
334    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
335    
336      dim_t i,j;
337      /*double sum;*/
338      index_t iptr;
339      index_t num_colors;
340      index_t* colorOf;
341      register index_t color;
342      bool_t passed=FALSE;
343      dim_t n=pattern->numOutput;
344    
345      
346      colorOf=MEMALLOC(n,index_t);
347    
348      if (pattern->type & PATTERN_FORMAT_SYM) {
349        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
350        return;
351      }
352      
353       Paso_Pattern_color(pattern,&num_colors,colorOf);
354      
355       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
356       #pragma omp parallel for private(i) schedule(static)
357       for (i=0;i<n;++i)
358            if(mis_marker[i]==IS_AVAILABLE)
359                        mis_marker[i]=IS_IN_SET;
360    
361       #pragma omp barrier
362       for (color=0;color<num_colors;++color) {
363        #pragma omp parallel for schedule(static) private(i,iptr,j)
364        for (i=0;i<n;++i) {
365         if (colorOf[i]==color) {  
366          if (mis_marker[i]==IS_IN_SET) {
367            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
368                 j=pattern->index[iptr];
369                 if (colorOf[j]<color)
370                  mis_marker[j]=IS_REMOVED;
371            }
372          }
373         }
374        }
375       }
376        
377        
378       #pragma omp barrier
379       for (color=0;color<num_colors;++color) {
380       #pragma omp parallel for schedule(static) private(i,iptr,j)
381        for (i=0;i<n;i++) {
382          if (colorOf[i]==color) {  
383            if (mis_marker[i]==IS_REMOVED) {
384               passed=TRUE;
385               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
386                  j=pattern->index[iptr];
387                   if (colorOf[j]<color && passed) {
388                    if (mis_marker[j]==IS_REMOVED) {
389                        passed=TRUE;
390                    }
391                    else {
392                        passed=FALSE;
393                        /*break;*/
394                    }
395                  }
396               }
397               }
398               if (passed) mis_marker[i]=IS_IN_SET;
399            }
400        }
401       }
402    
403         /* swap to TRUE/FALSE in mis_marker */
404         #pragma omp parallel for private(i) schedule(static)
405         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
406        
407        MEMFREE(colorOf);
408    }
409    
410    
411    

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

  ViewVC Help
Powered by ViewVC 1.1.26