/[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 1975 by artak, Thu Nov 6 03:07:02 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    
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
44    dim_t i,j;    dim_t i,j;
   /*double threshold=0.05;*/  
   index_t iptr,*index,*where_p,diagptr;  
   bool_t fail;  
   dim_t n=A->numRows;  
45    /*double sum;*/    /*double sum;*/
46      index_t iptr,*index,*where_p,*diagptr;
47      bool_t passed=FALSE;
48      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       if (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,index,where_p,diagptr,j) 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                   diagptr=A->pattern->ptr[i];          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
83                   index=&(A->pattern->index[diagptr]);               j=A->pattern->index[iptr];
84                   where_p=(index_t*)bsearch(&i,               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85                                          index,                  mis_marker[j]=IS_REMOVED;
86                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],               }
87                                          sizeof(index_t),          }
88                                          Paso_comparIndex);        }
89                  if (where_p==NULL) {      }
90                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");      
91                  } else {      
92                      diagptr+=(index_t)(where_p-index);      
93                  }        /*This loop cannot be parallelized, as order matters here.*/
94                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {      for (i=0;i<n;i++) {
95                       j=A->pattern->index[iptr]-index_offset;          if (mis_marker[i]==IS_REMOVED) {
96                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {             passed=TRUE;
97                          mis_marker[j]=IS_REMOVED;             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98                         }                j=A->pattern->index[iptr];
99                   }                if (mis_marker[j]==IS_IN_SET) {
100                    if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101                        passed=TRUE;
102                  }                  }
103              }                  else {
104                                    passed=FALSE;
105              #pragma omp parallel for private(i,sum,iptr) schedule(static)                      break;
106              for (i=0;i<n;++i)                  }
107                  if(mis_marker[i]==IS_AVAILABLE)                }
108                             mis_marker[i]=IS_IN_SET;             }
109                           if (passed) mis_marker[i]=IS_IN_SET;
110                #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)          }
111                for (i=0;i<n;i++) {      }
112                 if (mis_marker[i]==IS_REMOVED) {      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113                   fail=FALSE;      /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114                   diagptr=A->pattern->ptr[i];      /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115                   index=&(A->pattern->index[diagptr]);      for (i=0;i<n;i++) {
116                   where_p=(index_t*)bsearch(&i,          if (mis_marker[i]==IS_REMOVED) {
117                                          index,             sum=0;
118                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
119                                          sizeof(index_t),               j=A->pattern->index[iptr];
120                                          Paso_comparIndex);               if (mis_marker[j]==IS_REMOVED)
121                  if (where_p==NULL) {                  sum+=A->val[iptr];
122                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");             }
123                  } else {             if(ABS(sum)<1.e-25)
124                      diagptr+=(index_t)(where_p-index);               mis_marker[i]=IS_IN_SET;
                 }  
                  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 (mis_marker[j]==IS_IN_SET && (A->val[iptr]/A->val[diagptr])<-threshold){  
                          fail=TRUE;  
                          #ifndef _OPENMP    
                          break;  
                          #endif  
                      }  
                  }  
                  if(!fail) {  
                     mis_marker[i]=IS_IN_SET;  
                  }  
                }  
               }  
                 
          /*   #pragma omp parallel for private(i,sum,iptr) schedule(static)  
             for (i=0;i<n;++i)  
                 if(mis_marker[i]==IS_IN_SET)  
                     {  
                         sum=0.;  
                         for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {  
                             sum+=A->val[iptr];  
                         }  
                         if(ABS(sum)<ZERO)  
                            mis_marker[i]=IS_REMOVED;  
                     }  
          */    
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_SET);       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;  
143    index_t iptr;    index_t iptr;
144    double threshold,min_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)
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    
 /* is there any vertex available ?*/  
      if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
209    
210       #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)    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) {        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-index_offset] != 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_IN_SET;      #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(A->val[iptr-index_offset] < threshold){          eps_Aii = theta*theta*diags[i];
247                      mis_marker[i]=IS_REMOVED;          val=0.;
248                      #ifndef _OPENMP              for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
249                       break;              j=A->pattern->index[iptr];
250                      #endif              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_SET);       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    
   
 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
 {  
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
335    dim_t i,j;    dim_t i,j;
336      /*double sum;*/
337    index_t iptr;    index_t iptr;
338    double diag,eps_Aii,val;    index_t num_colors;
339    dim_t n=A->numRows;    index_t* colorOf;
340    double* diags=MEMALLOC(n,double);    register index_t color;
341      bool_t passed=FALSE;
342      dim_t n=pattern->numOutput;
343    
344        
345    if (A->pattern->type & PATTERN_FORMAT_SYM) {    colorOf=MEMALLOC(n,index_t);
346      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");  
347      if (pattern->type & PATTERN_FORMAT_SYM) {
348        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
349      return;      return;
350    }    }
351          
352     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {     Paso_Pattern_color(pattern,&num_colors,colorOf);
353      #pragma omp parallel for private(i,iptr,diag) schedule(static)    
354        for (i=0;i<n;++i) {     /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
355          diag = 0;     #pragma omp parallel for private(i) schedule(static)
356          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {     for (i=0;i<n;++i)
357              if(A->pattern->index[iptr-index_offset] != i){          if(mis_marker[i]==IS_AVAILABLE)
358                  diag+=A->val[iptr-index_offset];                      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          }          }
         diags[i]=ABS(diag);  
371        }        }
372         }
373        }
374       }
375            
376      #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)      
377        for (i=0;i<n;++i) {     #pragma omp barrier
378          eps_Aii = theta*theta*diags[i];     for (color=0;color<num_colors;++color) {
379          mis_marker[i]=IS_REMOVED;     #pragma omp parallel for schedule(static) private(i,iptr,j)
380        for (i=0;i<n;i++) {
381          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {        if (colorOf[i]==color) {  
382              j=A->pattern->index[iptr-index_offset];          if (mis_marker[i]==IS_REMOVED) {
383              val=A->val[iptr-index_offset];             passed=TRUE;
384              if(j!= i && val*val>=eps_Aii * diags[j]){             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
385                  mis_marker[i]=IS_IN_SET;                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      /* swap to TRUE/FALSE in mis_marker */  
402         /* swap to TRUE/FALSE in mis_marker */
403       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
404       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
405            
406       MEMFREE(diags);      MEMFREE(colorOf);
407  }  }
408    
409    
410    
411  #undef IS_AVAILABLE  #undef IS_AVAILABLE
412  #undef IS_IN_SET  #undef IS_IN_SET
413  #undef IS_REMOVED  #undef IS_REMOVED

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

  ViewVC Help
Powered by ViewVC 1.1.26