/[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 1872 by jfenwick, Mon Oct 13 00:18:55 2008 UTC revision 2381 by artak, Tue Apr 14 03:46:59 2009 UTC
# Line 28  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
 #include "mpi_C.h"  
 #include "Paso.h"  
31  #include "PasoUtil.h"  #include "PasoUtil.h"
32  #include "Pattern.h"  #include "Pattern_coupling.h"
33  #include "Solver.h"  #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
40  #define IS_IN_MIS -3  #define IS_REMOVED -4
41  #define IS_CONNECTED_TO_MIS -4  
42    void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43  void Paso_Pattern_cop(Paso_SparseMatrix* A, index_t* mis_marker) {  
44      dim_t i,j;
45    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    /*double threshold=0.05;*/
46    dim_t i;    double sum;
47    double threshold=0.05;    index_t iptr,*index,*where_p,*diagptr;
48    index_t naib,iptr;    bool_t passed=FALSE;
49    bool_t flag;    dim_t n=A->numRows;
50    dim_t n=A->pattern->numOutput;    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        
57       /* is there any vertex available ?*/     #pragma omp parallel for private(i) schedule(static)
58       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {     for (i=0;i<n;++i)
59            if(mis_marker[i]==IS_AVAILABLE)
60                        mis_marker[i]=IS_IN_SET;
61    
62             #pragma omp parallel for private(naib,i,iptr,flag) schedule(static)      #pragma omp parallel for private(i,index,where_p) schedule(static)
63             for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
64                if (mis_marker[i]==IS_AVAILABLE) {           diagptr[i]=A->pattern->ptr[i];
65                   flag=IS_IN_MIS_NOW;           index=&(A->pattern->index[diagptr[i]]);
66                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {           where_p=(index_t*)bsearch(&i,
67                       naib=A->pattern->index[iptr]-index_offset;                                  index,
68                       if (naib!=i && A->val[naib]<threshold*A->val[i]) {                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69                          flag=IS_AVAILABLE;                                  sizeof(index_t),
70                          break;                                  Paso_comparIndex);
71                       }          if (where_p==NULL) {
72                   }              Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73                   mis_marker[i]=flag;          } else {
74                }                  diagptr[i]+=(index_t)(where_p-index);
75             }          }
76        }
77        
78    
79             #pragma omp parallel for private(naib,i,iptr) schedule(static)  
80             for (i=0;i<n;i++) {      /*This loop cannot be parallelized, as order matters here.*/
81                if (mis_marker[i]==IS_AVAILABLE) {      for (i=0;i<n;++i) {
82                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {        if (mis_marker[i]==IS_IN_SET) {
83                       naib=A->pattern->index[iptr]-index_offset;          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
84                       if (naib!=i && mis_marker[i]==IS_IN_MIS_NOW && A->val[naib]/A->val[i]>=-threshold)               j=A->pattern->index[iptr];
85                           mis_marker[naib]=IS_IN_MIS;               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
86                   }                  mis_marker[j]=IS_REMOVED;
87                   mis_marker[i]=IS_CONNECTED_TO_MIS;               }
88                }          }
89          }
90        }
91        
92        
93        
94          /*This loop cannot be parallelized, as order matters here.*/
95        for (i=0;i<n;i++) {
96            if (mis_marker[i]==IS_REMOVED) {
97               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                        break;
106                    }
107                  }
108                }
109                if (passed) {
110                   mis_marker[i]=IS_IN_SET;
111                   passed=FALSE;
112                }
113            }
114        }
115    
116        /* 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'. */
118        #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119        for (i=0;i<n;i++) {
120            if (mis_marker[i]==IS_REMOVED) {
121               sum=0;
122               for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
123                 j=A->pattern->index[iptr];
124                 if (mis_marker[j]==IS_REMOVED)
125                    sum+=A->val[iptr];
126             }             }
127       }             if(ABS(sum)<1.e-25)
128                 mis_marker[i]=IS_IN_SET;
129            }
130        }
131    
132    
133       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
134       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
135       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);
136        
137         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_MIS_NOW  #undef IS_IN_SET
284  #undef IS_IN_MIS  #undef IS_REMOVED
 #undef IS_CONNECTED_TO_MIS  

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

  ViewVC Help
Powered by ViewVC 1.1.26