/[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 2309 by caltinay, Mon Mar 16 02:42:13 2009 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 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 30  Line 30 
30    
31  #include "PasoUtil.h"  #include "PasoUtil.h"
32  #include "Pattern_coupling.h"  #include "Pattern_coupling.h"
33    #include <limits.h>
34    
35    
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 93  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 104  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 coarsing 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 123  void Paso_Pattern_coup(Paso_SparseMatrix Line 120  void Paso_Pattern_coup(Paso_SparseMatrix
120               if (mis_marker[j]==IS_REMOVED)               if (mis_marker[j]==IS_REMOVED)
121                  sum+=A->val[iptr];                  sum+=A->val[iptr];
122             }             }
123             if(ABS(sum)<1.e-12)             if(ABS(sum)<1.e-25)
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 136  void Paso_Pattern_coup(Paso_SparseMatrix Line 133  void Paso_Pattern_coup(Paso_SparseMatrix
133       MEMFREE(diagptr);       MEMFREE(diagptr);
134  }  }
135    
136    /*
137     * Ruge-Stueben strength of connection mask.
138     *
139     */
140    void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
141    {
142      dim_t i,n,j;
143      index_t iptr;
144      double threshold,min_offdiagonal;
145      
146      Paso_Pattern *out=NULL;
147      
148      Paso_IndexList* index_list=NULL;
149    
150     index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
151       if (! Paso_checkPtr(index_list)) {
152            #pragma omp parallel for private(i) schedule(static)
153            for(i=0;i<A->pattern->numOutput;++i) {
154                 index_list[i].extension=NULL;
155                 index_list[i].n=0;
156            }
157        }
158      
159      
160      n=A->numRows;
161      if (A->pattern->type & PATTERN_FORMAT_SYM) {
162        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
163        return;
164      }
165        #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold,j) schedule(static)
166        for (i=0;i<n;++i) {
167          if(mis_marker[i]==IS_AVAILABLE) {
168            min_offdiagonal = DBL_MAX;
169            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170                if(A->pattern->index[iptr] != i){
171                    min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
172                }
173            }
174            
175            threshold = theta*min_offdiagonal;
176            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
177                j=A->pattern->index[iptr];
178                if(A->val[iptr]<=threshold) {
179                   if(j!=i) {
180                    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);
190        
191         /* clean up */
192       if (index_list!=NULL) {
193            #pragma omp parallel for private(i) schedule(static)
194            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
195         }
196        TMPMEMFREE(index_list);
197    
198        /*Paso_Pattern_mis(out,mis_marker);*/
199        Paso_Pattern_greedy(out,mis_marker);
200    }
201    
202    void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
203    {
204      dim_t i,j,n;
205      index_t iptr;
206      double diag,eps_Aii,val;
207      double* diags;
208    
209    
210      Paso_Pattern *out=NULL;
211      Paso_IndexList* index_list=NULL;
212    
213      n=A->numRows;  
214      diags=MEMALLOC(n,double);
215    
216      index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
217       if (! Paso_checkPtr(index_list)) {
218            #pragma omp parallel for private(i) schedule(static)
219            for(i=0;i<A->pattern->numOutput;++i) {
220                 index_list[i].extension=NULL;
221                 index_list[i].n=0;
222            }
223        }
224        
225      if (A->pattern->type & PATTERN_FORMAT_SYM) {
226        Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
227        return;
228      }
229    
230    
231        #pragma omp parallel for private(i,iptr) reduction(+:diag) schedule(static)
232          for (i=0;i<n;++i) {
233            diag = 0;
234            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
235                if(A->pattern->index[iptr] == i){
236                    diag+=A->val[iptr];
237                }
238            }
239            diags[i]=ABS(diag);
240          }
241    
242    
243        #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
244         for (i=0;i<n;++i) {
245           if (mis_marker[i]==IS_AVAILABLE) {
246            eps_Aii = theta*theta*diags[i];
247            val=0.;
248            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
249                j=A->pattern->index[iptr];
250                val=A->val[iptr];
251                if(j!= i) {
252                  if(val*val>=eps_Aii * diags[j]) {
253                   Paso_IndexList_insertIndex(&(index_list[i]),j);
254                  }
255                }
256            }
257           }
258         }
259    
260        out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
261        
262         /* clean up */
263        if (index_list!=NULL) {
264            #pragma omp parallel for private(i) schedule(static)
265            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
266         }
267    
268        TMPMEMFREE(index_list);
269        MEMFREE(diags);
270        
271        
272        /*Paso_Pattern_mis(out,mis_marker);*/
273        Paso_Pattern_greedy(out,mis_marker);
274    
275    }
276    
277    /* Greedy algorithm */
278    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
279    
280      dim_t i,j;
281      /*double sum;*/
282      index_t iptr;
283      bool_t passed=FALSE;
284      dim_t n=pattern->numOutput;
285    
286      if (pattern->type & PATTERN_FORMAT_SYM) {
287        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
288        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)
293       for (i=0;i<n;++i)
294            if(mis_marker[i]==IS_AVAILABLE)
295                        mis_marker[i]=IS_IN_SET;
296    
297    
298        for (i=0;i<n;++i) {
299          if (mis_marker[i]==IS_IN_SET) {
300            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
301                 j=pattern->index[iptr];
302                 mis_marker[j]=IS_REMOVED;
303            }
304          }
305        }
306        
307        
308        
309        for (i=0;i<n;i++) {
310            if (mis_marker[i]==IS_REMOVED) {
311               passed=TRUE;
312               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
313                  j=pattern->index[iptr];
314                    if (mis_marker[j]==IS_REMOVED) {
315                        passed=TRUE;
316                    }
317                    else {
318                        passed=FALSE;
319                        break;
320                    }
321                  }
322               }
323               if (passed) mis_marker[i]=IS_IN_SET;
324            }
325    
326         /* swap to TRUE/FALSE in mis_marker */
327         #pragma omp parallel for private(i) schedule(static)
328         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
329        
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
412  #undef IS_IN_SET  #undef IS_IN_SET

Legend:
Removed from v.2309  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26