/[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 1954 by artak, Fri Oct 31 03:22:34 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  /***************************************************************/  /***************************************************************/
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  
 #define ZERO 1.e-10  
   
41    
42  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43    
   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;*/
   index_t iptr,*index,*where_p,diagptr;  
   bool_t fail;  
   dim_t n=A->numRows;  
46    double sum;    double sum;
47      index_t iptr,*index,*where_p,*diagptr;
48      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       if (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    
80              #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)      /*This loop cannot be parallelized, as order matters here.*/
81              for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
82                if (mis_marker[i]==IS_AVAILABLE) {        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 (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                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {                  else {
104                       j=A->pattern->index[iptr]-index_offset;                      passed=FALSE;
105                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {                      break;
                         mis_marker[j]=IS_CONNECTED_TO_MIS;  
                        }  
                  }  
106                  }                  }
107                  }
108              }              }
109                            if (passed) {
110              #pragma omp parallel for private(i,sum,iptr) schedule(static)                 mis_marker[i]=IS_IN_SET;
111              for (i=0;i<n;++i)                 passed=FALSE;
112                  if(mis_marker[i]==IS_AVAILABLE)              }
113                             mis_marker[i]=IS_IN_MIS;          }
114                    }
115                #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)  
116                for (i=0;i<n;i++) {      /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117                 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {      /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
118                   fail=FALSE;      #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119                   diagptr=A->pattern->ptr[i];      for (i=0;i<n;i++) {
120                   index=&(A->pattern->index[diagptr]);          if (mis_marker[i]==IS_REMOVED) {
121                   where_p=(index_t*)bsearch(&i,             sum=0;
122                                          index,             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
123                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],               j=A->pattern->index[iptr];
124                                          sizeof(index_t),               if (mis_marker[j]==IS_REMOVED)
125                                          Paso_comparIndex);                  sum+=A->val[iptr];
126                  if (where_p==NULL) {             }
127                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");             if(ABS(sum)<1.e-25)
128                  } else {               mis_marker[i]=IS_IN_SET;
                     diagptr+=(index_t)(where_p-index);  
                 }  
                  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 (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){  
                          fail=TRUE;  
                          #ifndef _OPENMP    
                          break;  
                          #endif  
                      }  
                  }  
                  if(!fail) {  
                     mis_marker[i]=IS_IN_MIS;  
                  }  
                }  
               }  
                 
             #pragma omp parallel for private(i,sum,iptr) schedule(static)  
             for (i=0;i<n;++i)  
                 if(mis_marker[i]==IS_IN_MIS)  
                     {  
                         sum=0.;  
                         for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {  
                             sum+=A->val[iptr];  
                         }  
                         if(ABS(sum)<ZERO)  
                            mis_marker[i]=IS_CONNECTED_TO_MIS;  
                     }  
               
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;  
147    index_t iptr;    index_t iptr;
148    double threshold,min_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  /* 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) {
175          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];              if(A->pattern->index[iptr] != i){
176          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
             if(A->pattern->index[iptr-index_offset] != i){  
                 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);  
177              }              }
178          }          }
179    
180          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
181          mis_marker[i]=IS_IN_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(A->val[iptr-index_offset] < threshold){                 if(j!=i) {
185                      mis_marker[i]=IS_CONNECTED_TO_MIS;                  Paso_IndexList_insertIndex(&(index_list[i]),j);
186                      #ifndef _OPENMP                      }
                      break;  
                     #endif  
187              }              }
188          }          }
189           }
190        }        }
191      }      
192       /* swap to TRUE/FALSE in mis_marker */      out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
193       #pragma omp parallel for private(i) schedule(static)      
194       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);       /* 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)  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
205  {  {
206    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    dim_t i,j,n;
   dim_t i,j;  
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     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
233      #pragma omp parallel for private(i,iptr,diag) schedule(static)      #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]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
237              if(A->pattern->index[iptr-index_offset] != i){              if(A->pattern->index[iptr] != i){
238                  diag+=A->val[iptr-index_offset];                  diag+=A->val[iptr];
239              }              }
240          }          }
241          diags[i]=ABS(diag);          diags[i]=ABS(diag);
242        }        }
       
     #pragma omp parallel for private(i,iptr,diag) schedule(static)  
       for (i=0;i<n;++i) {  
         eps_Aii = theta*theta*diags[i];  
         mis_marker[i]=IS_CONNECTED_TO_MIS;  
243    
244          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {  
245              j=A->pattern->index[iptr-index_offset];      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
246              val=A->val[iptr-index_offset];       for (i=0;i<n;++i) {
247              if(j!= i && val*val>=eps_Aii * diags[j]){         if (mis_marker[i]==IS_AVAILABLE) {
248                  mis_marker[i]=IS_IN_MIS;          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     }       }
     /* swap to TRUE/FALSE in mis_marker */  
      #pragma omp parallel for private(i) schedule(static)  
      for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);  
       
      MEMFREE(diags);  
 }  
   
   
   
   
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  
 #undef ZERO  

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

  ViewVC Help
Powered by ViewVC 1.1.26