/[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 1980 by artak, Thu Nov 6 04:45:08 2008 UTC revision 2107 by artak, Fri Nov 28 04:39:07 2008 UTC
# Line 40  Line 40 
40    
41  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) {
42    
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
43    dim_t i,j;    dim_t i,j;
44    /*double threshold=0.05;*/    /*double threshold=0.05;*/
45    index_t iptr,*index,*where_p,diagptr;    index_t iptr,*index,*where_p,*diagptr;
46    bool_t fail;    bool_t passed=FALSE;
47    dim_t n=A->numRows;    dim_t n=A->numRows;
48      diagptr=MEMALLOC(n,index_t);
49    /*double sum;*/    /*double sum;*/
50    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
51      Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
52      return;      return;
53    }    }
54        
55       /* is there any vertex available ?*/     #pragma omp parallel for private(i) schedule(static)
56       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {     for (i=0;i<n;++i)
57            if(mis_marker[i]==IS_AVAILABLE)
58                        mis_marker[i]=IS_IN_SET;
59    
60              #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)      #pragma omp parallel for private(i,index,where_p) schedule(static)
61              for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
62                if (mis_marker[i]==IS_AVAILABLE) {           diagptr[i]=A->pattern->ptr[i];
63                   diagptr=A->pattern->ptr[i];           index=&(A->pattern->index[diagptr[i]]);
64                   index=&(A->pattern->index[diagptr]);           where_p=(index_t*)bsearch(&i,
65                   where_p=(index_t*)bsearch(&i,                                  index,
66                                          index,                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],
67                                          A->pattern->ptr[i + 1]-A->pattern->ptr[i],                                  sizeof(index_t),
68                                          sizeof(index_t),                                  Paso_comparIndex);
69                                          Paso_comparIndex);          if (where_p==NULL) {
70                  if (where_p==NULL) {              Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
71                      Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");          } else {
72                  } else {                  diagptr[i]+=(index_t)(where_p-index);
73                      diagptr+=(index_t)(where_p-index);          }
74                  }      }
75                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {      
76                       j=A->pattern->index[iptr]-index_offset;  
77                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {  
78                          mis_marker[j]=IS_REMOVED;      /*This loop cannot be parallelized, as order matters here.*/
79                         }      for (i=0;i<n;++i) {
80                   }        if (mis_marker[i]==IS_IN_SET) {
81            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
82                 j=A->pattern->index[iptr];
83                 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
84                    mis_marker[j]=IS_REMOVED;
85                 }
86            }
87          }
88        }
89        
90        
91        
92          /*This loop cannot be parallelized, as order matters here.*/
93        for (i=0;i<n;i++) {
94            if (mis_marker[i]==IS_REMOVED) {
95               for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
96                  j=A->pattern->index[iptr];
97                  if (mis_marker[j]==IS_IN_SET) {
98                    if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
99                        passed=TRUE;
100                  }                  }
101              }                  else {
102                                    passed=FALSE;
103              #pragma omp parallel for private(i,iptr) schedule(static)                      break;
             for (i=0;i<n;++i)  
                 if(mis_marker[i]==IS_AVAILABLE)  
                            mis_marker[i]=IS_IN_SET;  
               
               #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)  
               for (i=0;i<n;i++) {  
                if (mis_marker[i]==IS_REMOVED) {  
                  fail=FALSE;  
                  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);  
104                  }                  }
                  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_SET && (A->val[iptr]/A->val[diagptr])<-threshold){  
                          fail=TRUE;  
                          #ifndef _OPENMP    
                          break;  
                          #endif  
                      }  
                  }  
                  if(!fail) {  
                     mis_marker[i]=IS_IN_SET;  
                  }  
                }  
105                }                }
106                              }
107                if (passed) {
108                   mis_marker[i]=IS_IN_SET;
109                   passed=FALSE;
110                }
111          }          }
112        }
113          
114    
115       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
116       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
117       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
118        
119         MEMFREE(diagptr);
120  }  }
121    
122  /*  /*
# Line 136  void Paso_Pattern_coup(Paso_SparseMatrix Line 134  void Paso_Pattern_coup(Paso_SparseMatrix
134   */     */  
135  void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
136  {  {
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
137    dim_t i;    dim_t i;
138    index_t iptr;    index_t iptr;
139    double threshold,min_offdiagonal;    double threshold,max_offdiagonal;
140    dim_t n=A->numRows;    dim_t n=A->numRows;
141    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
142      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");
143      return;      return;
144    }    }
145    
146  /* is there any vertex available ?*/      #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold) schedule(static)
147       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {      for (i=0;i<n;++i) {
148          if(mis_marker[i]==IS_AVAILABLE) {
149       #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)          /*min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];*/
150        for (i=0;i<n;++i) {          max_offdiagonal = -1.e+15;
151          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
152          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {              if(A->pattern->index[iptr] != i){
153              if(A->pattern->index[iptr-index_offset] != i){                  max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));
                 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);  
154              }              }
155          }          }
156    
157          threshold = theta*min_offdiagonal;          threshold = theta*max_offdiagonal;
158          mis_marker[i]=IS_IN_SET;          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
159          #pragma omp parallel for private(iptr) schedule(static)              if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){
160          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {                      A->val[iptr]=0.;
             if(A->val[iptr-index_offset] < threshold){  
                     mis_marker[i]=IS_REMOVED;  
                     #ifndef _OPENMP      
                      break;  
                     #endif  
161              }              }
162          }          }
163           }
164        }        }
165      }  
166       /* swap to TRUE/FALSE in mis_marker */      Paso_Pattern_coup(A,mis_marker,0.05);
      #pragma omp parallel for private(i) schedule(static)  
      for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);  
167  }  }
168    
169    
# Line 181  void Paso_Pattern_RS(Paso_SparseMatrix* Line 171  void Paso_Pattern_RS(Paso_SparseMatrix*
171    
172  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)
173  {  {
   index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);  
174    dim_t i,j;    dim_t i,j;
175    index_t iptr;    index_t iptr;
176    double diag,eps_Aii,val;    double diag,eps_Aii,val;
# Line 193  void Paso_Pattern_Aggregiation(Paso_Spar Line 182  void Paso_Pattern_Aggregiation(Paso_Spar
182      return;      return;
183    }    }
184            
185     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
186      #pragma omp parallel for private(i,iptr,diag) schedule(static)      #pragma omp parallel for private(i,iptr,diag) schedule(static)
187        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
188          diag = 0;          diag = 0;
189          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) {
190              if(A->pattern->index[iptr-index_offset] != i){              if(A->pattern->index[iptr] != i){
191                  diag+=A->val[iptr-index_offset];                  diag+=A->val[iptr];
192              }              }
193          }          }
194          diags[i]=ABS(diag);          diags[i]=ABS(diag);
195        }        }
196            
197      #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)    
198        for (i=0;i<n;++i) {      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
199         for (i=0;i<n;++i) {
200           if (mis_marker[i]==IS_AVAILABLE) {
201          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
202          mis_marker[i]=IS_REMOVED;          val=0.;
203            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
204          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {              j=A->pattern->index[iptr];
205              j=A->pattern->index[iptr-index_offset];              val=A->val[iptr];
206              val=A->val[iptr-index_offset];              if(j!= i) {
207              if(j!= i && val*val>=eps_Aii * diags[j]){                if(val*val>=eps_Aii * diags[j]) {
208                  mis_marker[i]=IS_IN_SET;                 A->val[iptr]=val;
209                  }
210                  else {
211                   A->val[iptr]=0.;
212                  }
213              }              }
214          }          }
215        }         }
216     }       }
217      /* swap to TRUE/FALSE in mis_marker */      
218       #pragma omp parallel for private(i) schedule(static)      Paso_Pattern_coup(A,mis_marker,0.05);
      for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);  
       
219       MEMFREE(diags);       MEMFREE(diags);
220  }  }
221    
# Line 230  void Paso_Pattern_Aggregiation(Paso_Spar Line 223  void Paso_Pattern_Aggregiation(Paso_Spar
223  #undef IS_AVAILABLE  #undef IS_AVAILABLE
224  #undef IS_IN_SET  #undef IS_IN_SET
225  #undef IS_REMOVED  #undef IS_REMOVED
226    
227    void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {
228      index_t out=0, *mis_marker=NULL,i;
229      dim_t n=A->numRows;
230      mis_marker=TMPMEMALLOC(n,index_t);
231      if ( !Paso_checkPtr(mis_marker) ) {
232        /* get coloring */
233        #pragma omp parallel for private(i) schedule(static)
234        for (i = 0; i < n; ++i) {
235            colorOf[i]=-1;
236            mis_marker[i]=-1;
237        }
238    
239        while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {
240           Paso_Pattern_coup(A,mis_marker,0.05);
241    
242           #pragma omp parallel for private(i) schedule(static)
243           for (i = 0; i < n; ++i) {
244            if (mis_marker[i]) colorOf[i]=out;
245            mis_marker[i]=colorOf[i];
246           }
247           ++out;
248        }
249        TMPMEMFREE(mis_marker);
250        *num_colors=out;
251      }
252    }
253    
254    

Legend:
Removed from v.1980  
changed lines
  Added in v.2107

  ViewVC Help
Powered by ViewVC 1.1.26