/[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 1907 by phornby, Wed Oct 22 11:41:18 2008 UTC revision 1975 by artak, Thu Nov 6 03:07:02 2008 UTC
# Line 28  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
 #include "mpi_C.h"  
 #include "Paso.h"  
31  #include "PasoUtil.h"  #include "PasoUtil.h"
32  #include "Pattern.h"  #include "Pattern_coupling.h"
 #include "Solver.h"  
33    
34    
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  
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) {  
42    
43    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    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;    index_t iptr,*index,*where_p,diagptr;
47    bool_t flag;    bool_t fail;
48    dim_t n=A->pattern->numOutput;    dim_t n=A->numRows;
49      /*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 ?*/       /* is there any vertex available ?*/
56       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
57    
58             #pragma omp parallel for private(i,iptr,flag) schedule(static)              #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
59             for (i=0;i<n;++i) {              for (i=0;i<n;++i) {
60                if (mis_marker[i]==IS_AVAILABLE) {                if (mis_marker[i]==IS_AVAILABLE) {
                  flag=IS_IN_MIS;  
61                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
62                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
63                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 78  void Paso_Pattern_coup(Paso_SparseMatrix Line 72  void Paso_Pattern_coup(Paso_SparseMatrix
72                  }                  }
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 && A->val[iptr]>=threshold*A->val[diagptr]) {                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
76                          flag=IS_AVAILABLE;                          mis_marker[j]=IS_REMOVED;
77                          break;                         }
                      }  
78                   }                   }
                  mis_marker[i]=flag;  
79                  }                  }
80              }              }
81                          
82                #pragma omp parallel for private(i,iptr) schedule(static)              #pragma omp parallel for private(i,sum,iptr) schedule(static)
83                for (i=0;i<n;++i)
84                    if(mis_marker[i]==IS_AVAILABLE)
85                               mis_marker[i]=IS_IN_SET;
86                
87                  #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_AVAILABLE) {                 if (mis_marker[i]==IS_REMOVED) {
90                     fail=FALSE;
91                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
92                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
93                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 104  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 (j!=i && 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                           mis_marker[i]=IS_IN_MIS;                           fail=TRUE;
107                       }                           #ifndef _OPENMP  
108                       else {                           break;
109                           mis_marker[i]=IS_CONNECTED_TO_MIS;                           #endif
110                       }                       }
111                   }                   }
112                     if(!fail) {
113                        mis_marker[i]=IS_IN_SET;
114                     }
115                 }                 }
116                }                }
117                  
118             /*   #pragma omp parallel for private(i,sum,iptr) schedule(static)
119                for (i=0;i<n;++i)
120                    if(mis_marker[i]==IS_IN_SET)
121                        {
122                            sum=0.;
123                            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
124                                sum+=A->val[iptr];
125                            }
126                            if(ABS(sum)<ZERO)
127                               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 135  void Paso_Pattern_coup(Paso_SparseMatrix Line 149  void Paso_Pattern_coup(Paso_SparseMatrix
149  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)
150  {  {
151    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
152    dim_t i,j;    dim_t i;
153    index_t iptr;    index_t iptr;
154    double threshold,min_offdiagonal;    double threshold,min_offdiagonal;
155    bool_t flag;    dim_t n=A->numRows;
   dim_t n=A->pattern->numOutput;  
156    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
157      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");
158      return;      return;
# Line 152  void Paso_Pattern_RS(Paso_SparseMatrix* Line 165  void Paso_Pattern_RS(Paso_SparseMatrix*
165        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
166          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];          min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
167          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) {
168              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr-index_offset] != i){
169                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
170              }              }
171          }          }
172    
173          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
174          mis_marker[i]=IS_CONNECTED_TO_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(-1.*(A->val[iptr-index_offset]) <= threshold){              if(A->val[iptr-index_offset] < threshold){
178                      mis_marker[i]=IS_IN_MIS;                      mis_marker[i]=IS_REMOVED;
179                        #ifndef _OPENMP    
180                         break;
181                        #endif
182              }              }
183          }          }
184        }        }
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    
192    
193    
194    void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
195    {
196      index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
197      dim_t i,j;
198      index_t iptr;
199      double diag,eps_Aii,val;
200      dim_t n=A->numRows;
201      double* diags=MEMALLOC(n,double);
202      
203      if (A->pattern->type & PATTERN_FORMAT_SYM) {
204        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
205        return;
206      }
207        
208       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
209        #pragma omp parallel for private(i,iptr,diag) schedule(static)
210          for (i=0;i<n;++i) {
211            diag = 0;
212            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
213                if(A->pattern->index[iptr-index_offset] != i){
214                    diag+=A->val[iptr-index_offset];
215                }
216            }
217            diags[i]=ABS(diag);
218          }
219        
220        #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)
221          for (i=0;i<n;++i) {
222            eps_Aii = theta*theta*diags[i];
223            mis_marker[i]=IS_REMOVED;
224    
225            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];
227                val=A->val[iptr-index_offset];
228                if(j!= i && val*val>=eps_Aii * diags[j]){
229                    mis_marker[i]=IS_IN_SET;
230                }
231            }
232          }
233       }
234        /* swap to TRUE/FALSE in mis_marker */
235         #pragma omp parallel for private(i) schedule(static)
236         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
237        
238         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  

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

  ViewVC Help
Powered by ViewVC 1.1.26