/[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 1974 by artak, Fri Oct 31 03:22:34 2008 UTC revision 1975 by artak, Thu Nov 6 03:07:02 2008 UTC
# Line 35  Line 35 
35  /***************************************************************/  /***************************************************************/
36    
37  #define IS_AVAILABLE -1  #define IS_AVAILABLE -1
38  #define IS_IN_MIS_NOW -2  #define IS_IN_SET -3
39  #define IS_IN_MIS -3  #define IS_REMOVED -4
 #define IS_CONNECTED_TO_MIS -4  
 #define ZERO 1.e-10  
   
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    
# Line 49  void Paso_Pattern_coup(Paso_SparseMatrix Line 46  void Paso_Pattern_coup(Paso_SparseMatrix
46    index_t iptr,*index,*where_p,diagptr;    index_t iptr,*index,*where_p,diagptr;
47    bool_t fail;    bool_t fail;
48    dim_t n=A->numRows;    dim_t n=A->numRows;
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;
# Line 76  void Paso_Pattern_coup(Paso_SparseMatrix Line 73  void Paso_Pattern_coup(Paso_SparseMatrix
73                   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) {
74                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
75                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
76                          mis_marker[j]=IS_CONNECTED_TO_MIS;                          mis_marker[j]=IS_REMOVED;
77                         }                         }
78                   }                   }
79                  }                  }
# Line 85  void Paso_Pattern_coup(Paso_SparseMatrix Line 82  void Paso_Pattern_coup(Paso_SparseMatrix
82              #pragma omp parallel for private(i,sum,iptr) schedule(static)              #pragma omp parallel for private(i,sum,iptr) schedule(static)
83              for (i=0;i<n;++i)              for (i=0;i<n;++i)
84                  if(mis_marker[i]==IS_AVAILABLE)                  if(mis_marker[i]==IS_AVAILABLE)
85                             mis_marker[i]=IS_IN_MIS;                             mis_marker[i]=IS_IN_SET;
86                            
87                #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)                #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
88                for (i=0;i<n;i++) {                for (i=0;i<n;i++) {
89                 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {                 if (mis_marker[i]==IS_REMOVED) {
90                   fail=FALSE;                   fail=FALSE;
91                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
92                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
# Line 105  void Paso_Pattern_coup(Paso_SparseMatrix Line 102  void Paso_Pattern_coup(Paso_SparseMatrix
102                  }                  }
103                   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) {
104                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
105                       if (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){                       if (mis_marker[j]==IS_IN_SET && (A->val[iptr]/A->val[diagptr])<-threshold){
106                           fail=TRUE;                           fail=TRUE;
107                           #ifndef _OPENMP                             #ifndef _OPENMP  
108                           break;                           break;
# Line 113  void Paso_Pattern_coup(Paso_SparseMatrix Line 110  void Paso_Pattern_coup(Paso_SparseMatrix
110                       }                       }
111                   }                   }
112                   if(!fail) {                   if(!fail) {
113                      mis_marker[i]=IS_IN_MIS;                      mis_marker[i]=IS_IN_SET;
114                   }                   }
115                 }                 }
116                }                }
117                                
118              #pragma omp parallel for private(i,sum,iptr) schedule(static)           /*   #pragma omp parallel for private(i,sum,iptr) schedule(static)
119              for (i=0;i<n;++i)              for (i=0;i<n;++i)
120                  if(mis_marker[i]==IS_IN_MIS)                  if(mis_marker[i]==IS_IN_SET)
121                      {                      {
122                          sum=0.;                          sum=0.;
123                          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) {
124                              sum+=A->val[iptr];                              sum+=A->val[iptr];
125                          }                          }
126                          if(ABS(sum)<ZERO)                          if(ABS(sum)<ZERO)
127                             mis_marker[i]=IS_CONNECTED_TO_MIS;                             mis_marker[i]=IS_REMOVED;
128                      }                      }
129                         */  
130          }          }
131       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
132       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
133       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);
134  }  }
135    
136  /*  /*
# Line 174  void Paso_Pattern_RS(Paso_SparseMatrix* Line 171  void Paso_Pattern_RS(Paso_SparseMatrix*
171          }          }
172    
173          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
174          mis_marker[i]=IS_IN_MIS;          mis_marker[i]=IS_IN_SET;
175          #pragma omp parallel for private(iptr) schedule(static)          #pragma omp parallel for private(iptr) schedule(static)
176          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) {
177              if(A->val[iptr-index_offset] < threshold){              if(A->val[iptr-index_offset] < threshold){
178                      mis_marker[i]=IS_CONNECTED_TO_MIS;                      mis_marker[i]=IS_REMOVED;
179                      #ifndef _OPENMP                          #ifndef _OPENMP    
180                       break;                       break;
181                      #endif                      #endif
# Line 188  void Paso_Pattern_RS(Paso_SparseMatrix* Line 185  void Paso_Pattern_RS(Paso_SparseMatrix*
185      }      }
186       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
187       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
188       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);
189  }  }
190    
191    
# Line 220  void Paso_Pattern_Aggregiation(Paso_Spar Line 217  void Paso_Pattern_Aggregiation(Paso_Spar
217          diags[i]=ABS(diag);          diags[i]=ABS(diag);
218        }        }
219            
220      #pragma omp parallel for private(i,iptr,diag) schedule(static)      #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)
221        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
222          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
223          mis_marker[i]=IS_CONNECTED_TO_MIS;          mis_marker[i]=IS_REMOVED;
224    
225          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) {
226              j=A->pattern->index[iptr-index_offset];              j=A->pattern->index[iptr-index_offset];
227              val=A->val[iptr-index_offset];              val=A->val[iptr-index_offset];
228              if(j!= i && val*val>=eps_Aii * diags[j]){              if(j!= i && val*val>=eps_Aii * diags[j]){
229                  mis_marker[i]=IS_IN_MIS;                  mis_marker[i]=IS_IN_SET;
230              }              }
231          }          }
232        }        }
233     }     }
234      /* swap to TRUE/FALSE in mis_marker */      /* swap to TRUE/FALSE in mis_marker */
235       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
236       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);
237            
238       MEMFREE(diags);       MEMFREE(diags);
239  }  }
240    
241    
   
   
   
   
   
   
   
   
   
   
   
242  #undef IS_AVAILABLE  #undef IS_AVAILABLE
243  #undef IS_IN_MIS_NOW  #undef IS_IN_SET
244  #undef IS_IN_MIS  #undef IS_REMOVED
 #undef IS_CONNECTED_TO_MIS  
 #undef ZERO  

Legend:
Removed from v.1974  
changed lines
  Added in v.1975

  ViewVC Help
Powered by ViewVC 1.1.26