/[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 1914 by phornby, Thu Oct 23 08:31:22 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 28  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
 /*  
 #include "mpi_C.h"  
 #include "Paso.h"  
 #include "Pattern.h"  
 */  
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_MIS_NOW -2  #define IS_IN_SET -3   /* Week connection */
40  #define IS_IN_MIS -3  #define IS_REMOVED -4  /* strong */
 #define IS_CONNECTED_TO_MIS -4  
   
41    
42    void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43    
 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {  
   
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
44    dim_t i,j;    dim_t i,j;
45    double threshold=0.05;    /*double sum;*/
46    index_t iptr,*index,*where_p,diagptr;    index_t iptr,*index,*where_p,*diagptr;
47    bool_t flag;    bool_t passed=FALSE;
48    dim_t n=A->pattern->numOutput;    dim_t n=A->numRows;
49      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        
56       /* is there any vertex available ?*/     #pragma omp parallel for private(i) schedule(static)
57       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {     for (i=0;i<n;++i)
58            if(mis_marker[i]==IS_AVAILABLE)
59                        mis_marker[i]=IS_IN_SET;
60    
61        #pragma omp parallel for private(i,index,where_p) schedule(static)
62        for (i=0;i<n;++i) {
63             diagptr[i]=A->pattern->ptr[i];
64             index=&(A->pattern->index[diagptr[i]]);
65             where_p=(index_t*)bsearch(&i,
66                                    index,
67                                    A->pattern->ptr[i + 1]-A->pattern->ptr[i],
68                                    sizeof(index_t),
69                                    Paso_comparIndex);
70            if (where_p==NULL) {
71                Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72            } else {
73                    diagptr[i]+=(index_t)(where_p-index);
74            }
75        }
76        
77    
78    
79             #pragma omp parallel for private(i,iptr,flag) schedule(static)      /*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_AVAILABLE) {        if (mis_marker[i]==IS_IN_SET) {
82                   flag=IS_IN_MIS;          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
83                   diagptr=A->pattern->ptr[i];               j=A->pattern->index[iptr];
84                   index=&(A->pattern->index[diagptr]);               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85                   where_p=(index_t*)bsearch(&i,                  mis_marker[j]=IS_REMOVED;
86                                          index,               }
87                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],          }
88                                          sizeof(index_t),        }
89                                          Paso_comparIndex);      }
90                  if (where_p==NULL) {      
91                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");      
92                  } else {      
93                      diagptr+=(index_t)(where_p-index);        /*This loop cannot be parallelized, as order matters here.*/
94                  }      for (i=0;i<n;i++) {
95                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {          if (mis_marker[i]==IS_REMOVED) {
96                       j=A->pattern->index[iptr]-index_offset;             passed=TRUE;
97                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98                          flag=IS_AVAILABLE;                j=A->pattern->index[iptr];
99                          break;                if (mis_marker[j]==IS_IN_SET) {
100                       }                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101                   }                      passed=TRUE;
102                   mis_marker[i]=flag;                  }
103                    else {
104                        passed=FALSE;
105                        break;
106                  }                  }
             }  
             
               #pragma omp parallel for private(i,iptr) schedule(static)  
               for (i=0;i<n;i++) {  
                if (mis_marker[i]==IS_AVAILABLE) {  
                  diagptr=A->pattern->ptr[i];  
                  index=&(A->pattern->index[diagptr]);  
                  where_p=(index_t*)bsearch(&i,  
                                         index,  
                                         A->pattern->ptr[i + 1]-A->pattern->ptr[i],  
                                         sizeof(index_t),  
                                         Paso_comparIndex);  
                 if (where_p==NULL) {  
                     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");  
                 } else {  
                     diagptr+=(index_t)(where_p-index);  
                 }  
                  for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {  
                      j=A->pattern->index[iptr]-index_offset;  
                      if (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){  
                          mis_marker[i]=IS_IN_MIS;  
                      }  
                      else {  
                          mis_marker[i]=IS_CONNECTED_TO_MIS;  
                      }  
                  }  
                }  
107                }                }
108               }
109               if (passed) mis_marker[i]=IS_IN_SET;
110            }
111        }
112        /* 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'. */
114        /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115        for (i=0;i<n;i++) {
116            if (mis_marker[i]==IS_REMOVED) {
117               sum=0;
118               for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
119                 j=A->pattern->index[iptr];
120                 if (mis_marker[j]==IS_REMOVED)
121                    sum+=A->val[iptr];
122               }
123               if(ABS(sum)<1.e-25)
124                 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)
131       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
132        
133         MEMFREE(diagptr);
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    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    dim_t i,n,j;
   dim_t i,j;  
143    index_t iptr;    index_t iptr;
144    double threshold,min_offdiagonal;    double threshold,min_offdiagonal;
145    bool_t flag;    
146    dim_t n=A->pattern->numOutput;    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)
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    
 /* is there any vertex available ?*/  
      if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
230    
231       #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)      #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          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];          diag = 0;
234          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++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                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);                  diag+=A->val[iptr];
237              }              }
238          }          }
239            diags[i]=ABS(diag);
240          }
241    
242          threshold = theta*min_offdiagonal;  
243          mis_marker[i]=IS_CONNECTED_TO_MIS;      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
244          #pragma omp parallel for private(iptr) schedule(static)       for (i=0;i<n;++i) {
245          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {         if (mis_marker[i]==IS_AVAILABLE) {
246              if(-1.*(A->val[iptr-index_offset]) <= threshold){          eps_Aii = theta*theta*diags[i];
247                      mis_marker[i]=IS_IN_MIS;          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 */       /* swap to TRUE/FALSE in mis_marker */
327       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
328       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);       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_MIS_NOW  #undef IS_IN_SET
413  #undef IS_IN_MIS  #undef IS_REMOVED
 #undef IS_CONNECTED_TO_MIS  

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

  ViewVC Help
Powered by ViewVC 1.1.26