/[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 2726 by artak, Wed Oct 21 23:50:05 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_F -3   /* in F (strong) */
40  #define IS_REMOVED -4  #define IS_IN_C -4  /* in C (weak) */
41    
42  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  void Paso_Pattern_YS(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 57  void Paso_Pattern_coup(Paso_SparseMatrix Line 56  void Paso_Pattern_coup(Paso_SparseMatrix
56     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
57     for (i=0;i<n;++i)     for (i=0;i<n;++i)
58          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
59                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_IN_C;
60    
61      #pragma omp parallel for private(i,index,where_p) schedule(static)      /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
62      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
63           diagptr[i]=A->pattern->ptr[i];           diagptr[i]=A->pattern->ptr[i];
64           index=&(A->pattern->index[diagptr[i]]);           index=&(A->pattern->index[A->pattern->ptr[i]]);
65           where_p=(index_t*)bsearch(&i,           where_p=(index_t*)bsearch(&i,
66                                  index,                                  index,
67                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],
# Line 79  void Paso_Pattern_coup(Paso_SparseMatrix Line 78  void Paso_Pattern_coup(Paso_SparseMatrix
78    
79      /*This loop cannot be parallelized, as order matters here.*/      /*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_IN_SET) {        if (mis_marker[i]==IS_IN_C) {
82          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) {
83               j=A->pattern->index[iptr];               j=A->pattern->index[iptr];
84               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85                  mis_marker[j]=IS_REMOVED;                  mis_marker[j]=IS_IN_F;
86               }               }
87          }          }
88        }        }
# Line 93  void Paso_Pattern_coup(Paso_SparseMatrix Line 92  void Paso_Pattern_coup(Paso_SparseMatrix
92            
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_IN_F) {
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_C) {
100                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101                      passed=TRUE;                      passed=TRUE;
102                  }                  }
# Line 105  void Paso_Pattern_coup(Paso_SparseMatrix Line 105  void Paso_Pattern_coup(Paso_SparseMatrix
105                      break;                      break;
106                  }                  }
107                }                }
             }  
             if (passed) {  
                mis_marker[i]=IS_IN_SET;  
                passed=FALSE;  
             }  
         }  
     }  
   
     /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/  
     /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */  
     #pragma omp parallel for private(i,iptr,j,sum) schedule(static)  
     for (i=0;i<n;i++) {  
         if (mis_marker[i]==IS_REMOVED) {  
            sum=0;  
            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
              j=A->pattern->index[iptr];  
              if (mis_marker[j]==IS_REMOVED)  
                 sum+=A->val[iptr];  
108             }             }
109             if(ABS(sum)<1.e-25)             if (passed) mis_marker[i]=IS_IN_C;
              mis_marker[i]=IS_IN_SET;  
110          }          }
111      }      }
112    
   
113       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
114       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
115       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_F);
116            
117       MEMFREE(diagptr);       MEMFREE(diagptr);
118  }  }
119    
120    
121  /*  /*
122   * Ruge-Stueben strength of connection mask.   * Ruge-Stueben strength of connection mask.
123   *   *
# Line 145  void Paso_Pattern_RS(Paso_SparseMatrix* Line 126  void Paso_Pattern_RS(Paso_SparseMatrix*
126  {  {
127    dim_t i,n,j;    dim_t i,n,j;
128    index_t iptr;    index_t iptr;
129    double threshold,min_offdiagonal;    double threshold,max_offdiagonal;
130        
131    Paso_Pattern *out=NULL;    Paso_Pattern *out=NULL;
132        
# Line 166  void Paso_Pattern_RS(Paso_SparseMatrix* Line 147  void Paso_Pattern_RS(Paso_SparseMatrix*
147      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");
148      return;      return;
149    }    }
150        /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
     #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)  
151      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
152        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
153          min_offdiagonal = DBL_MAX;          max_offdiagonal = DBL_MIN;
154          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) {
155              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
156                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);                  max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
157              }              }
158          }          }
159            
160          threshold = theta*min_offdiagonal;          threshold = theta*max_offdiagonal;
161          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) {
162              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
163              if(A->val[iptr]<=threshold) {              if((-A->val[iptr])>=threshold) {
                if(j!=i) {  
164                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
165                  }                  Paso_IndexList_insertIndex(&(index_list[j]),i);
166              }              }
167          }          }
168         }         }
169        }        }
170            
171      
172      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);
173            
174       /* clean up */       /* clean up */
# Line 198  void Paso_Pattern_RS(Paso_SparseMatrix* Line 178  void Paso_Pattern_RS(Paso_SparseMatrix*
178       }       }
179      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
180    
181      Paso_Pattern_mis(out,mis_marker);      /*Paso_Pattern_mis(out,mis_marker);*/
182        Paso_Pattern_greedy(out,mis_marker);
183        Paso_Pattern_free(out);
184  }  }
185    
186  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 207  void Paso_Pattern_Aggregiation(Paso_Spar
207      }      }
208            
209    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
210      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");
211      return;      return;
212    }    }
213    
# Line 234  void Paso_Pattern_Aggregiation(Paso_Spar Line 216  void Paso_Pattern_Aggregiation(Paso_Spar
216        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
217          diag = 0;          diag = 0;
218          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) {
219              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] == i){
220                  diag+=A->val[iptr];                  diag+=A->val[iptr];
221              }              }
222          }          }
# Line 246  void Paso_Pattern_Aggregiation(Paso_Spar Line 228  void Paso_Pattern_Aggregiation(Paso_Spar
228       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
229         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
230          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
         val=0.;  
231          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) {
232              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
233              val=A->val[iptr];              val=A->val[iptr];
234              if(j!= i) {                if((val*val)>=(eps_Aii*diags[j])) {
               if(val*val>=eps_Aii * diags[j]) {  
235                 Paso_IndexList_insertIndex(&(index_list[i]),j);                 Paso_IndexList_insertIndex(&(index_list[i]),j);
236                }                }
             }  
237          }          }
238         }         }
239       }       }
# Line 270  void Paso_Pattern_Aggregiation(Paso_Spar Line 249  void Paso_Pattern_Aggregiation(Paso_Spar
249      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
250      MEMFREE(diags);      MEMFREE(diags);
251            
252      Paso_Pattern_mis(out,mis_marker);      
253        /*Paso_Pattern_mis(out,mis_marker);*/
254        Paso_Pattern_greedy(out,mis_marker);
255        Paso_Pattern_free(out);
256    
257    }
258    
259    /* Greedy algorithm */
260    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
261    
262      dim_t i,j;
263      /*double sum;*/
264      index_t iptr;
265      bool_t passed=FALSE;
266      dim_t n=pattern->numOutput;
267    
268      if (pattern->type & PATTERN_FORMAT_SYM) {
269        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
270        return;
271      }
272      
273       #pragma omp parallel for private(i) schedule(static)
274       for (i=0;i<n;++i)
275            if(mis_marker[i]==IS_AVAILABLE)
276                        mis_marker[i]=IS_IN_C;
277    
278    
279        for (i=0;i<n;++i) {
280          if (mis_marker[i]==IS_IN_C) {
281            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
282                 j=pattern->index[iptr];
283                 mis_marker[j]=IS_IN_F;
284            }
285          }
286        }
287        
288        
289        
290        for (i=0;i<n;i++) {
291            if (mis_marker[i]==IS_IN_F) {
292               passed=TRUE;
293               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
294                  j=pattern->index[iptr];
295                    if (mis_marker[j]==IS_IN_F) {
296                        passed=TRUE;
297                    }
298                    else {
299                        passed=FALSE;
300                        break;
301                    }
302                  }
303               if (passed) mis_marker[i]=IS_IN_C;
304               }
305            }
306    
307         /* swap to TRUE/FALSE in mis_marker */
308         #pragma omp parallel for private(i) schedule(static)
309         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
310        
311  }  }
312    
313    
314    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
315    
316      dim_t i,j;
317      /*double sum;*/
318      index_t iptr;
319      index_t num_colors;
320      index_t* colorOf;
321      register index_t color;
322      bool_t passed=FALSE;
323      dim_t n=pattern->numOutput;
324    
325      
326      colorOf=MEMALLOC(n,index_t);
327    
328      if (pattern->type & PATTERN_FORMAT_SYM) {
329        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
330        return;
331      }
332      
333       Paso_Pattern_color(pattern,&num_colors,colorOf);
334      
335       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
336       #pragma omp parallel for private(i) schedule(static)
337       for (i=0;i<n;++i)
338            if(mis_marker[i]==IS_AVAILABLE)
339                        mis_marker[i]=IS_IN_F;
340    
341       #pragma omp barrier
342       for (color=0;color<num_colors;++color) {
343        #pragma omp parallel for schedule(static) private(i,iptr,j)
344        for (i=0;i<n;++i) {
345         if (colorOf[i]==color) {  
346          if (mis_marker[i]==IS_IN_F) {
347            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
348                 j=pattern->index[iptr];
349                 if (colorOf[j]<color)
350                  mis_marker[j]=IS_IN_C;
351            }
352          }
353         }
354        }
355       }
356        
357        
358       #pragma omp barrier
359       for (color=0;color<num_colors;++color) {
360       #pragma omp parallel for schedule(static) private(i,iptr,j)
361        for (i=0;i<n;i++) {
362          if (colorOf[i]==color) {  
363            if (mis_marker[i]==IS_IN_C) {
364               passed=TRUE;
365               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
366                  j=pattern->index[iptr];
367                   if (colorOf[j]<color && passed) {
368                    if (mis_marker[j]==IS_IN_C) {
369                        passed=TRUE;
370                    }
371                    else {
372                        passed=FALSE;
373                        /*break;*/
374                    }
375                  }
376               }
377               if (passed) mis_marker[i]=IS_IN_F;
378               }
379              
380            }
381        }
382       }
383    
384         /* swap to TRUE/FALSE in mis_marker */
385         #pragma omp parallel for private(i) schedule(static)
386         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
387        
388        MEMFREE(colorOf);
389    }
390    
391    /*For testing */
392    void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
393    
394      dim_t i,j=0,k;
395      double *theta;
396      index_t iptr;
397      dim_t n=A->numRows;
398      double rsum,diag=0;
399      index_t *AvADJ;
400      theta=MEMALLOC(n,double);
401      AvADJ=MEMALLOC(n,index_t);
402    
403    
404      
405    
406      if (A->pattern->type & PATTERN_FORMAT_SYM) {
407        Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
408        return;
409      }
410      
411    
412        #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
413        for (i=0;i<n;++i) {
414            rsum=0;
415            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
416                j=A->pattern->index[iptr];
417                if(j!=i) {
418                  rsum+=ABS(A->val[iptr]);    
419                }
420                else {
421                    diag=ABS(A->val[iptr]);
422                }
423            }
424            theta[i]=diag/rsum;
425            if(theta[i]>threshold) {
426                mis_marker[i]=IS_IN_F;
427            }
428        }
429        
430        while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
431             k=0;
432             #pragma omp parallel for private(i,j,k) schedule(static)
433             for (i=0;i<n;++i) {
434               if(mis_marker[i]==IS_AVAILABLE) {
435                    if(k==0) {
436                        j=i;
437                        k++;
438                    }
439                    if(theta[j]>theta[i]) {
440                        j=i;
441                    }
442                }
443             }
444             mis_marker[j]=IS_IN_C;
445            
446             for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
447                k=A->pattern->index[iptr];
448                if(mis_marker[k]==IS_AVAILABLE) {
449                   AvADJ[k]=1;
450                }
451                else {
452                    AvADJ[k]=-1;
453                }
454                
455             }
456                
457             #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
458            for (i=0;i<n;++i) {
459                if(AvADJ[i]) {
460                    rsum=0;
461                    for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
462                        k=A->pattern->index[iptr];
463                        if(k!=i && mis_marker[k]!=IS_IN_C ) {
464                          rsum+=ABS(A->val[iptr]);    
465                        }
466                        if(j==i) {
467                            diag=ABS(A->val[iptr]);
468                        }
469                    }
470                    theta[i]=diag/rsum;
471                    if(theta[i]>threshold) {
472                       mis_marker[i]=IS_IN_F;
473                    }
474                }
475            }
476            
477            
478        }
479    
480         /* swap to TRUE/FALSE in mis_marker */
481         #pragma omp parallel for private(i) schedule(static)
482         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
483        
484         MEMFREE(AvADJ);
485         MEMFREE(theta);
486    }
487    
488  #undef IS_AVAILABLE  #undef IS_AVAILABLE
489  #undef IS_IN_SET  #undef IS_IN_F
490  #undef IS_REMOVED  #undef IS_IN_C
491    
492    
493    
494    
495    
496    
497    
498    
499    
500    

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

  ViewVC Help
Powered by ViewVC 1.1.26