/[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 1914 by phornby, Thu Oct 23 08:31:22 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"  
 #include "Pattern.h"  
 */  
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_MIS_NOW -2  #define IS_IN_SET -3
40  #define IS_IN_MIS -3  #define IS_REMOVED -4
 #define IS_CONNECTED_TO_MIS -4  
   
41    
42    void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43    
 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {  
   
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
44    dim_t i,j;    dim_t i,j;
45    double threshold=0.05;    /*double threshold=0.05;*/
46    index_t iptr,*index,*where_p,diagptr;    double sum;
47    bool_t flag;    index_t iptr,*index,*where_p,*diagptr;
48    dim_t n=A->pattern->numOutput;    bool_t passed=FALSE;
49      dim_t n=A->numRows;
50      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(i,index,where_p) schedule(static)
63        for (i=0;i<n;++i) {
64             diagptr[i]=A->pattern->ptr[i];
65             index=&(A->pattern->index[diagptr[i]]);
66             where_p=(index_t*)bsearch(&i,
67                                    index,
68                                    A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69                                    sizeof(index_t),
70                                    Paso_comparIndex);
71            if (where_p==NULL) {
72                Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73            } else {
74                    diagptr[i]+=(index_t)(where_p-index);
75            }
76        }
77        
78    
79             #pragma omp parallel for private(i,iptr,flag) 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                   flag=IS_IN_MIS;        if (mis_marker[i]==IS_IN_SET) {
83                   diagptr=A->pattern->ptr[i];          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
84                   index=&(A->pattern->index[diagptr]);               j=A->pattern->index[iptr];
85                   where_p=(index_t*)bsearch(&i,               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
86                                          index,                  mis_marker[j]=IS_REMOVED;
87                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],               }
88                                          sizeof(index_t),          }
89                                          Paso_comparIndex);        }
90                  if (where_p==NULL) {      }
91                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");      
92                  } else {      
93                      diagptr+=(index_t)(where_p-index);      
94                  }        /*This loop cannot be parallelized, as order matters here.*/
95                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {      for (i=0;i<n;i++) {
96                       j=A->pattern->index[iptr]-index_offset;          if (mis_marker[i]==IS_REMOVED) {
97                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98                          flag=IS_AVAILABLE;                j=A->pattern->index[iptr];
99                          break;                if (mis_marker[j]==IS_IN_SET) {
100                       }                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101                   }                      passed=TRUE;
                  mis_marker[i]=flag;  
102                  }                  }
103              }                  else {
104                                  passed=FALSE;
105                #pragma omp parallel for private(i,iptr) schedule(static)                      break;
               for (i=0;i<n;i++) {  
                if (mis_marker[i]==IS_AVAILABLE) {  
                  diagptr=A->pattern->ptr[i];  
                  index=&(A->pattern->index[diagptr]);  
                  where_p=(index_t*)bsearch(&i,  
                                         index,  
                                         A->pattern->ptr[i + 1]-A->pattern->ptr[i],  
                                         sizeof(index_t),  
                                         Paso_comparIndex);  
                 if (where_p==NULL) {  
                     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");  
                 } else {  
                     diagptr+=(index_t)(where_p-index);  
106                  }                  }
                  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 (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){  
                          mis_marker[i]=IS_IN_MIS;  
                      }  
                      else {  
                          mis_marker[i]=IS_CONNECTED_TO_MIS;  
                      }  
                  }  
                }  
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   * 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    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    dim_t i,n,j;
   dim_t i,j;  
147    index_t iptr;    index_t iptr;
148    double threshold,min_offdiagonal;    double threshold,min_offdiagonal;
149    bool_t flag;    
150    dim_t n=A->pattern->numOutput;    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  /* is there any vertex available ?*/      #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
171       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {      for (i=0;i<n;++i) {
172          if(mis_marker[i]==IS_AVAILABLE) {
173       #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)          min_offdiagonal = DBL_MAX;
174        for (i=0;i<n;++i) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
         min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];  
         for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {  
175              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
176                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
177              }              }
178          }          }
179    
180          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
181          mis_marker[i]=IS_CONNECTED_TO_MIS;          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
182          #pragma omp parallel for private(iptr) schedule(static)              j=A->pattern->index[iptr];
183          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {              if(A->val[iptr]<=threshold) {
184              if(-1.*(A->val[iptr-index_offset]) <= threshold){                 if(j!=i) {
185                      mis_marker[i]=IS_IN_MIS;                  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       /* swap to TRUE/FALSE in mis_marker */      
227       #pragma omp parallel for private(i) schedule(static)    if (A->pattern->type & PATTERN_FORMAT_SYM) {
228       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);      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.1914  
changed lines
  Added in v.2381

  ViewVC Help
Powered by ViewVC 1.1.26