/[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 2380 by caltinay, Mon Mar 16 02:42:13 2009 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 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      #pragma omp parallel for private(i,iptr,j,sum) schedule(static)      #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
# 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 136  void Paso_Pattern_coup(Paso_SparseMatrix Line 137  void Paso_Pattern_coup(Paso_SparseMatrix
137       MEMFREE(diagptr);       MEMFREE(diagptr);
138  }  }
139    
140    /*
141     * Ruge-Stueben strength of connection mask.
142     *
143     */
144    void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
145    {
146      dim_t i,n,j;
147      index_t iptr;
148      double threshold,min_offdiagonal;
149      
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) {
166        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
167        return;
168      }
169    
170        #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
171        for (i=0;i<n;++i) {
172          if(mis_marker[i]==IS_AVAILABLE) {
173            min_offdiagonal = DBL_MAX;
174            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
175                if(A->pattern->index[iptr] != i){
176                    min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
177                }
178            }
179    
180            threshold = theta*min_offdiagonal;
181            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
182                j=A->pattern->index[iptr];
183                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_mis(out,mis_marker);
202    }
203    
204    void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
205    {
206      dim_t i,j,n;
207      index_t iptr;
208      double diag,eps_Aii,val;
209      double* diags;
210    
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) {
228        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
229        return;
230      }
231    
232    
233        #pragma omp parallel for private(i,iptr,diag) schedule(static)
234          for (i=0;i<n;++i) {
235            diag = 0;
236            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
237                if(A->pattern->index[iptr] != i){
238                    diag+=A->val[iptr];
239                }
240            }
241            diags[i]=ABS(diag);
242          }
243    
244    
245        #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
246         for (i=0;i<n;++i) {
247           if (mis_marker[i]==IS_AVAILABLE) {
248            eps_Aii = theta*theta*diags[i];
249            val=0.;
250            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
251                j=A->pattern->index[iptr];
252                val=A->val[iptr];
253                if(j!= i) {
254                  if(val*val>=eps_Aii * diags[j]) {
255                   Paso_IndexList_insertIndex(&(index_list[i]),j);
256                  }
257                }
258            }
259           }
260         }
261    
262        out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
263        
264         /* clean up */
265        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
Legend:
Removed from v.2380  
changed lines
  Added in v.2381
  ViewVC Help
Powered by ViewVC 1.1.26