/[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 2122 by artak, Wed Dec 3 02:52:28 2008 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;
49    diagptr=MEMALLOC(n,index_t);    diagptr=MEMALLOC(n,index_t);
50    
51    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
52      Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
53      return;      return;
54    }    }
55        
# Line 68  void Paso_Pattern_coup(Paso_SparseMatrix Line 68  void Paso_Pattern_coup(Paso_SparseMatrix
68                                  sizeof(index_t),                                  sizeof(index_t),
69                                  Paso_comparIndex);                                  Paso_comparIndex);
70          if (where_p==NULL) {          if (where_p==NULL) {
71              Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");              Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72          } else {          } else {
73                  diagptr[i]+=(index_t)(where_p-index);                  diagptr[i]+=(index_t)(where_p-index);
74          }          }
# 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 103  void Paso_Pattern_coup(Paso_SparseMatrix Line 104  void Paso_Pattern_coup(Paso_SparseMatrix
104                      passed=FALSE;                      passed=FALSE;
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 137  void Paso_Pattern_coup(Paso_SparseMatrix Line 134  void Paso_Pattern_coup(Paso_SparseMatrix
134  }  }
135    
136  /*  /*
137     * Ruge-Stueben strength of connection mask.
138   *   *
139   * Return a strength of connection mask using the classical   */
140   * strength of connection measure by Ruge and Stuben.  void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
  *  
  * Specifically, an off-diagonal entry A[i.j] is a strong  
  * connection if:  
  *    
  *      -A[i,j] >= theta * max( -A[i,k] )   where k != i  
  *  
  * Otherwise, the connection is weak.  
  *    
  */    
 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
141  {  {
142    dim_t i;    dim_t i,n,j;
143    index_t iptr;    index_t iptr;
144    double threshold,max_offdiagonal;    double threshold,min_offdiagonal;
145    dim_t n=A->numRows;    
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) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
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) 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 = A->val[A->pattern->ptr[i]-index_offset];*/          min_offdiagonal = DBL_MAX;
         max_offdiagonal = -1.e+15;  
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                  max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
172              }              }
173          }          }
174            
175          threshold = theta*max_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              if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){              j=A->pattern->index[iptr];
178                      A->val[iptr]=0.;              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_coup(A,mis_marker,0.05);      /*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)
   
   
 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
203  {  {
204    dim_t i,j;    dim_t i,j,n;
205    index_t iptr;    index_t iptr;
206    double diag,eps_Aii,val;    double diag,eps_Aii,val;
207    dim_t n=A->numRows;    double* diags;
208    double* diags=MEMALLOC(n,double);  
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) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
226      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");
227      return;      return;
228    }    }
       
229    
230      #pragma omp parallel for private(i,iptr,diag) schedule(static)  
231        #pragma omp parallel for private(i,iptr) reduction(+:diag) schedule(static)
232        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
233          diag = 0;          diag = 0;
234          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) {
235              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] == i){
236                  diag+=A->val[iptr];                  diag+=A->val[iptr];
237              }              }
238          }          }
239          diags[i]=ABS(diag);          diags[i]=ABS(diag);
240        }        }
241        
242      
243      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
244       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
245         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
246          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
247          val=0.;          val=0.;
248          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) {
249              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
250              val=A->val[iptr];              val=A->val[iptr];
251              if(j!= i) {              if(j!= i) {
252                if(val*val>=eps_Aii * diags[j]) {                if(val*val>=eps_Aii * diags[j]) {
253                 A->val[iptr]=val;                 Paso_IndexList_insertIndex(&(index_list[i]),j);
               }  
               else {  
                A->val[iptr]=0.;  
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_coup(A,mis_marker,0.05);      /*Paso_Pattern_mis(out,mis_marker);*/
273       MEMFREE(diags);      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  #undef IS_AVAILABLE    dim_t i,j;
281  #undef IS_IN_SET    /*double sum;*/
282  #undef IS_REMOVED    index_t iptr;
283      bool_t passed=FALSE;
284      dim_t n=pattern->numOutput;
285    
286  void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {    if (pattern->type & PATTERN_FORMAT_SYM) {
287    index_t out=0, *mis_marker=NULL,i;      Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
288    dim_t n=A->numRows;      return;
289    mis_marker=TMPMEMALLOC(n,index_t);    }
290    if ( !Paso_checkPtr(mis_marker) ) {    
291      /* get coloring */     /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
292      #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
293      for (i = 0; i < n; ++i) {     for (i=0;i<n;++i)
294          colorOf[i]=-1;          if(mis_marker[i]==IS_AVAILABLE)
295          mis_marker[i]=-1;                      mis_marker[i]=IS_IN_SET;
     }  
296    
     while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {  
        Paso_Pattern_coup(A,mis_marker,0.05);  
297    
298         #pragma omp parallel for private(i) schedule(static)      for (i=0;i<n;++i) {
299         for (i = 0; i < n; ++i) {        if (mis_marker[i]==IS_IN_SET) {
300          if (mis_marker[i]) colorOf[i]=out;          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
301          mis_marker[i]=colorOf[i];               j=pattern->index[iptr];
302         }               mis_marker[j]=IS_REMOVED;
303         ++out;          }
304          }
305      }      }
306      TMPMEMFREE(mis_marker);      
307      *num_colors=out;      
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
412    #undef IS_IN_SET
413    #undef IS_REMOVED

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

  ViewVC Help
Powered by ViewVC 1.1.26