/[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 1913 by phornby, Thu Oct 23 08:27:33 2008 UTC revision 1954 by artak, Fri Oct 31 03:22:34 2008 UTC
# Line 28  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
 /*  
 #include "mpi_C.h"  
 #include "Paso.h"  
31  #include "PasoUtil.h"  #include "PasoUtil.h"
 #include "Pattern.h"  
 */  
32  #include "Pattern_coupling.h"  #include "Pattern_coupling.h"
33    
34    
# Line 43  Line 38 
38  #define IS_IN_MIS_NOW -2  #define IS_IN_MIS_NOW -2
39  #define IS_IN_MIS -3  #define IS_IN_MIS -3
40  #define IS_CONNECTED_TO_MIS -4  #define IS_CONNECTED_TO_MIS -4
41    #define ZERO 1.e-10
42    
43    
44    void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {  
45    
46    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
47    dim_t i,j;    dim_t i,j;
48    double threshold=0.05;    /*double threshold=0.05;*/
49    index_t iptr,*index,*where_p,diagptr;    index_t iptr,*index,*where_p,diagptr;
50    bool_t flag;    bool_t fail;
51    dim_t n=A->pattern->numOutput;    dim_t n=A->numRows;
52      double sum;
53    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
54      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");
55      return;      return;
56    }    }
57        
58       /* is there any vertex available ?*/       /* is there any vertex available ?*/
59       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
60    
61             #pragma omp parallel for private(i,iptr,flag) schedule(static)              #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
62             for (i=0;i<n;++i) {              for (i=0;i<n;++i) {
63                if (mis_marker[i]==IS_AVAILABLE) {                if (mis_marker[i]==IS_AVAILABLE) {
                  flag=IS_IN_MIS;  
64                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
65                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
66                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 80  void Paso_Pattern_coup(Paso_SparseMatrix Line 75  void Paso_Pattern_coup(Paso_SparseMatrix
75                  }                  }
76                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
77                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
78                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
79                          flag=IS_AVAILABLE;                          mis_marker[j]=IS_CONNECTED_TO_MIS;
80                          break;                         }
                      }  
81                   }                   }
                  mis_marker[i]=flag;  
82                  }                  }
83              }              }
84                          
85                #pragma omp parallel for private(i,iptr) schedule(static)              #pragma omp parallel for private(i,sum,iptr) schedule(static)
86                for (i=0;i<n;++i)
87                    if(mis_marker[i]==IS_AVAILABLE)
88                               mis_marker[i]=IS_IN_MIS;
89                
90                  #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
91                for (i=0;i<n;i++) {                for (i=0;i<n;i++) {
92                 if (mis_marker[i]==IS_AVAILABLE) {                 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {
93                     fail=FALSE;
94                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
95                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
96                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 106  void Paso_Pattern_coup(Paso_SparseMatrix Line 105  void Paso_Pattern_coup(Paso_SparseMatrix
105                  }                  }
106                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {                   for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
107                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
108                       if (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){                       if (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){
109                           mis_marker[i]=IS_IN_MIS;                           fail=TRUE;
110                       }                           #ifndef _OPENMP  
111                       else {                           break;
112                           mis_marker[i]=IS_CONNECTED_TO_MIS;                           #endif
113                       }                       }
114                   }                   }
115                     if(!fail) {
116                        mis_marker[i]=IS_IN_MIS;
117                     }
118                 }                 }
119                }                }
120                  
121                #pragma omp parallel for private(i,sum,iptr) schedule(static)
122                for (i=0;i<n;++i)
123                    if(mis_marker[i]==IS_IN_MIS)
124                        {
125                            sum=0.;
126                            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
127                                sum+=A->val[iptr];
128                            }
129                            if(ABS(sum)<ZERO)
130                               mis_marker[i]=IS_CONNECTED_TO_MIS;
131                        }
132                
133          }          }
134       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
135       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
# Line 137  void Paso_Pattern_coup(Paso_SparseMatrix Line 152  void Paso_Pattern_coup(Paso_SparseMatrix
152  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)
153  {  {
154    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
155    dim_t i,j;    dim_t i;
156    index_t iptr;    index_t iptr;
157    double threshold,min_offdiagonal;    double threshold,min_offdiagonal;
158    bool_t flag;    dim_t n=A->numRows;
   dim_t n=A->pattern->numOutput;  
159    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
160      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");
161      return;      return;
# Line 154  void Paso_Pattern_RS(Paso_SparseMatrix* Line 168  void Paso_Pattern_RS(Paso_SparseMatrix*
168        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
169          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
170          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
171              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr-index_offset] != i){
172                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
173              }              }
174          }          }
175    
176          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
177          mis_marker[i]=IS_CONNECTED_TO_MIS;          mis_marker[i]=IS_IN_MIS;
178          #pragma omp parallel for private(iptr) schedule(static)          #pragma omp parallel for private(iptr) schedule(static)
179          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {          for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
180              if(-1.*(A->val[iptr-index_offset]) <= threshold){              if(A->val[iptr-index_offset] < threshold){
181                      mis_marker[i]=IS_IN_MIS;                      mis_marker[i]=IS_CONNECTED_TO_MIS;
182                        #ifndef _OPENMP    
183                         break;
184                        #endif
185              }              }
186          }          }
187        }        }
# Line 173  void Paso_Pattern_RS(Paso_SparseMatrix* Line 190  void Paso_Pattern_RS(Paso_SparseMatrix*
190       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
191       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_MIS);
192  }  }
193    
194    
195    
196    
197    void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
198    {
199      index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
200      dim_t i,j;
201      index_t iptr;
202      double diag,eps_Aii,val;
203      dim_t n=A->numRows;
204      double* diags=MEMALLOC(n,double);
205      
206      if (A->pattern->type & PATTERN_FORMAT_SYM) {
207        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
208        return;
209      }
210        
211       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
212        #pragma omp parallel for private(i,iptr,diag) schedule(static)
213          for (i=0;i<n;++i) {
214            diag = 0;
215            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
216                if(A->pattern->index[iptr-index_offset] != i){
217                    diag+=A->val[iptr-index_offset];
218                }
219            }
220            diags[i]=ABS(diag);
221          }
222        
223        #pragma omp parallel for private(i,iptr,diag) schedule(static)
224          for (i=0;i<n;++i) {
225            eps_Aii = theta*theta*diags[i];
226            mis_marker[i]=IS_CONNECTED_TO_MIS;
227    
228            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
229                j=A->pattern->index[iptr-index_offset];
230                val=A->val[iptr-index_offset];
231                if(j!= i && val*val>=eps_Aii * diags[j]){
232                    mis_marker[i]=IS_IN_MIS;
233                }
234            }
235          }
236       }
237        /* swap to TRUE/FALSE in mis_marker */
238         #pragma omp parallel for private(i) schedule(static)
239         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
240        
241         MEMFREE(diags);
242    }
243    
244    
245    
246    
247    
248    
249    
250    
251    
252    
253    
254    
255    
256  #undef IS_AVAILABLE  #undef IS_AVAILABLE
257  #undef IS_IN_MIS_NOW  #undef IS_IN_MIS_NOW
258  #undef IS_IN_MIS  #undef IS_IN_MIS
259  #undef IS_CONNECTED_TO_MIS  #undef IS_CONNECTED_TO_MIS
260    #undef ZERO

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

  ViewVC Help
Powered by ViewVC 1.1.26