/[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 2131 by artak, Thu Dec 4 06:09:50 2008 UTC revision 2381 by artak, Tue Apr 14 03:46:59 2009 UTC
# 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  /***************************************************************/  /***************************************************************/
# Line 49  void Paso_Pattern_coup(Paso_SparseMatrix Line 50  void Paso_Pattern_coup(Paso_SparseMatrix
50    diagptr=MEMALLOC(n,index_t);    diagptr=MEMALLOC(n,index_t);
51    
52    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
53      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");
54      return;      return;
55    }    }
56        
# Line 68  void Paso_Pattern_coup(Paso_SparseMatrix Line 69  void Paso_Pattern_coup(Paso_SparseMatrix
69                                  sizeof(index_t),                                  sizeof(index_t),
70                                  Paso_comparIndex);                                  Paso_comparIndex);
71          if (where_p==NULL) {          if (where_p==NULL) {
72              Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");              Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73          } else {          } else {
74                  diagptr[i]+=(index_t)(where_p-index);                  diagptr[i]+=(index_t)(where_p-index);
75          }          }
# 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) {
110                 mis_marker[i]=IS_IN_SET;                 mis_marker[i]=IS_IN_SET;
# Line 112  void Paso_Pattern_coup(Paso_SparseMatrix Line 113  void Paso_Pattern_coup(Paso_SparseMatrix
113          }          }
114      }      }
115    
116      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsing process.*/      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117      /* 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'. */
118       /*This loop cannot be parallelized, as order matters here.*/      #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
120          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_REMOVED) {
121             sum=0;             sum=0;
# Line 123  void Paso_Pattern_coup(Paso_SparseMatrix Line 124  void Paso_Pattern_coup(Paso_SparseMatrix
124               if (mis_marker[j]==IS_REMOVED)               if (mis_marker[j]==IS_REMOVED)
125                  sum+=A->val[iptr];                  sum+=A->val[iptr];
126             }             }
127             if(ABS(sum)<1.e-12)             if(ABS(sum)<1.e-25)
128               mis_marker[i]=IS_IN_SET;               mis_marker[i]=IS_IN_SET;
129          }          }
130      }      }
# Line 137  void Paso_Pattern_coup(Paso_SparseMatrix Line 138  void Paso_Pattern_coup(Paso_SparseMatrix
138  }  }
139    
140  /*  /*
141     * Ruge-Stueben strength of connection mask.
142   *   *
143   * Return a strength of connection mask using the classical   */
144   * 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)  
145  {  {
146    dim_t i;    dim_t i,n,j;
147    index_t iptr;    index_t iptr;
148    double threshold,max_offdiagonal;    double threshold,min_offdiagonal;
149    dim_t n=A->numRows;    
150      Paso_Pattern *out=NULL;
151      
152      Paso_IndexList* index_list=NULL;
153    
154     index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
155       if (! Paso_checkPtr(index_list)) {
156            #pragma omp parallel for private(i) schedule(static)
157            for(i=0;i<A->pattern->numOutput;++i) {
158                 index_list[i].extension=NULL;
159                 index_list[i].n=0;
160            }
161        }
162      
163      
164      n=A->numRows;
165    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
166      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");
167      return;      return;
168    }    }
169    
170      #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold) schedule(static)      #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
171      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
172        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
173          /*min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];*/          min_offdiagonal = DBL_MAX;
         max_offdiagonal = -1.e+15;  
174          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) {
175              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
176                  max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
177              }              }
178          }          }
179    
180          threshold = theta*max_offdiagonal;          threshold = theta*min_offdiagonal;
181          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) {
182              if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){              j=A->pattern->index[iptr];
183                      A->val[iptr]=0.;              if(A->val[iptr]<=threshold) {
184                   if(j!=i) {
185                    Paso_IndexList_insertIndex(&(index_list[i]),j);
186                    }
187              }              }
188          }          }
189         }         }
190        }        }
191        
192        out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
193        
194         /* clean up */
195       if (index_list!=NULL) {
196            #pragma omp parallel for private(i) schedule(static)
197            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
198         }
199        TMPMEMFREE(index_list);
200    
201      Paso_Pattern_coup(A,mis_marker,0.05);      Paso_Pattern_mis(out,mis_marker);
202  }  }
203    
204    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)  
205  {  {
206    dim_t i,j;    dim_t i,j,n;
207    index_t iptr;    index_t iptr;
208    double diag,eps_Aii,val;    double diag,eps_Aii,val;
209    dim_t n=A->numRows;    double* diags;
210    double* diags=MEMALLOC(n,double);  
211      
212      Paso_Pattern *out=NULL;
213      Paso_IndexList* index_list=NULL;
214    
215      n=A->numRows;  
216      diags=MEMALLOC(n,double);
217    
218      index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
219       if (! Paso_checkPtr(index_list)) {
220            #pragma omp parallel for private(i) schedule(static)
221            for(i=0;i<A->pattern->numOutput;++i) {
222                 index_list[i].extension=NULL;
223                 index_list[i].n=0;
224            }
225        }
226        
227    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
228      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");
229      return;      return;
230    }    }
       
231    
232      #pragma omp parallel for private(i,iptr,diag) schedule(static)  
233        #pragma omp parallel for private(i,iptr,diag) schedule(static)
234        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
235          diag = 0;          diag = 0;
236          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) {
# Line 210  void Paso_Pattern_Aggregiation(Paso_Spar Line 240  void Paso_Pattern_Aggregiation(Paso_Spar
240          }          }
241          diags[i]=ABS(diag);          diags[i]=ABS(diag);
242        }        }
243        
244      
245      #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)
246       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
247         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
248          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
249          val=0.;          val=0.;
250          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) {
251              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
252              val=A->val[iptr];              val=A->val[iptr];
253              if(j!= i) {              if(j!= i) {
254                if(val*val>=eps_Aii * diags[j]) {                if(val*val>=eps_Aii * diags[j]) {
255                 A->val[iptr]=val;                 Paso_IndexList_insertIndex(&(index_list[i]),j);
               }  
               else {  
                A->val[iptr]=0.;  
256                }                }
257              }              }
258          }          }
259         }         }
260       }       }
261    
262        out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
263            
264      Paso_Pattern_coup(A,mis_marker,0.05);       /* clean up */
265       MEMFREE(diags);      if (index_list!=NULL) {
266            #pragma omp parallel for private(i) schedule(static)
267            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
268         }
269    
270        TMPMEMFREE(index_list);
271        MEMFREE(diags);
272        
273        Paso_Pattern_mis(out,mis_marker);
274    
275    
276  }  }
277    
278    
279    
280    
281    
282  #undef IS_AVAILABLE  #undef IS_AVAILABLE
283  #undef IS_IN_SET  #undef IS_IN_SET
284  #undef IS_REMOVED  #undef IS_REMOVED
   
 void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {  
   index_t out=0, *mis_marker=NULL,i;  
   dim_t n=A->numRows;  
   mis_marker=TMPMEMALLOC(n,index_t);  
   if ( !Paso_checkPtr(mis_marker) ) {  
     /* get coloring */  
     #pragma omp parallel for private(i) schedule(static)  
     for (i = 0; i < n; ++i) {  
         colorOf[i]=-1;  
         mis_marker[i]=-1;  
     }  
   
     while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {  
        Paso_Pattern_coup(A,mis_marker,0.05);  
   
        #pragma omp parallel for private(i) schedule(static)  
        for (i = 0; i < n; ++i) {  
         if (mis_marker[i]) colorOf[i]=out;  
         mis_marker[i]=colorOf[i];  
        }  
        ++out;  
     }  
     TMPMEMFREE(mis_marker);  
     *num_colors=out;  
   }  
 }  

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

  ViewVC Help
Powered by ViewVC 1.1.26