/[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 1882 by artak, Tue Oct 14 04:39:02 2008 UTC revision 1980 by artak, Thu Nov 6 04:45:08 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) {  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
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 naib,iptr;    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(naib,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) {
61                   flag=IS_AVAILABLE;                   diagptr=A->pattern->ptr[i];
62                     index=&(A->pattern->index[diagptr]);
63                     where_p=(index_t*)bsearch(&i,
64                                            index,
65                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
66                                            sizeof(index_t),
67                                            Paso_comparIndex);
68                    if (where_p==NULL) {
69                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
70                    } else {
71                        diagptr+=(index_t)(where_p-index);
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                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
75                       if (naib!=i && A->val[naib]<threshold*A->val[i]) {                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
76                          flag=IS_IN_MIS;                          mis_marker[j]=IS_REMOVED;
77                          break;                         }
                      }  
78                   }                   }
79                   mis_marker[i]=flag;                  }
80                }              }
81             }              
82                #pragma omp parallel for private(i,iptr) schedule(static)
83             #pragma omp parallel for private(naib,i,iptr) schedule(static)              for (i=0;i<n;++i)
84             for (i=0;i<n;i++) {                  if(mis_marker[i]==IS_AVAILABLE)
85                if (mis_marker[i]==IS_AVAILABLE) {                             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++) {
89                   if (mis_marker[i]==IS_REMOVED) {
90                     fail=FALSE;
91                     diagptr=A->pattern->ptr[i];
92                     index=&(A->pattern->index[diagptr]);
93                     where_p=(index_t*)bsearch(&i,
94                                            index,
95                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
96                                            sizeof(index_t),
97                                            Paso_comparIndex);
98                    if (where_p==NULL) {
99                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
100                    } else {
101                        diagptr+=(index_t)(where_p-index);
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                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
105                        if (naib!=i && mis_marker[naib]==IS_IN_MIS && A->val[iptr]/A->val[i]>=-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                       else                           #ifndef _OPENMP  
108                           mis_marker[i]=IS_CONNECTED_TO_MIS;                           break;
109                             #endif
110                         }
111                     }
112                     if(!fail) {
113                        mis_marker[i]=IS_IN_SET;
114                   }                   }
115                   }
116                }                }
117             }                
118       }          }
119       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
120       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
121       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);
122  }  }
123    
124    /*
125     *
126     * Return a strength of connection mask using the classical
127     * strength of connection measure by Ruge and Stuben.
128     *
129     * Specifically, an off-diagonal entry A[i.j] is a strong
130     * connection if:
131     *  
132     *      -A[i,j] >= theta * max( -A[i,k] )   where k != i
133     *
134     * Otherwise, the connection is weak.
135     *  
136     */  
137    void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
138    {
139      index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
140      dim_t i;
141      index_t iptr;
142      double threshold,min_offdiagonal;
143      dim_t n=A->numRows;
144      if (A->pattern->type & PATTERN_FORMAT_SYM) {
145        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
146        return;
147      }
148    
149    /* is there any vertex available ?*/
150         if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
151    
152         #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
153          for (i=0;i<n;++i) {
154            min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
155            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
156                if(A->pattern->index[iptr-index_offset] != i){
157                    min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
158                }
159            }
160    
161            threshold = theta*min_offdiagonal;
162            mis_marker[i]=IS_IN_SET;
163            #pragma omp parallel for private(iptr) schedule(static)
164            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
165                if(A->val[iptr-index_offset] < threshold){
166                        mis_marker[i]=IS_REMOVED;
167                        #ifndef _OPENMP    
168                         break;
169                        #endif
170                }
171            }
172          }
173        }
174         /* swap to TRUE/FALSE in mis_marker */
175         #pragma omp parallel for private(i) schedule(static)
176         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
177    }
178    
179    
180    
181    
182    void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
183    {
184      index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
185      dim_t i,j;
186      index_t iptr;
187      double diag,eps_Aii,val;
188      dim_t n=A->numRows;
189      double* diags=MEMALLOC(n,double);
190      
191      if (A->pattern->type & PATTERN_FORMAT_SYM) {
192        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
193        return;
194      }
195        
196       if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
197        #pragma omp parallel for private(i,iptr,diag) schedule(static)
198          for (i=0;i<n;++i) {
199            diag = 0;
200            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
201                if(A->pattern->index[iptr-index_offset] != i){
202                    diag+=A->val[iptr-index_offset];
203                }
204            }
205            diags[i]=ABS(diag);
206          }
207        
208        #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)
209          for (i=0;i<n;++i) {
210            eps_Aii = theta*theta*diags[i];
211            mis_marker[i]=IS_REMOVED;
212    
213            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
214                j=A->pattern->index[iptr-index_offset];
215                val=A->val[iptr-index_offset];
216                if(j!= i && val*val>=eps_Aii * diags[j]){
217                    mis_marker[i]=IS_IN_SET;
218                }
219            }
220          }
221       }
222        /* swap to TRUE/FALSE in mis_marker */
223         #pragma omp parallel for private(i) schedule(static)
224         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
225        
226         MEMFREE(diags);
227    }
228    
229    
230  #undef IS_AVAILABLE  #undef IS_AVAILABLE
231  #undef IS_IN_MIS_NOW  #undef IS_IN_SET
232  #undef IS_IN_MIS  #undef IS_REMOVED
 #undef IS_CONNECTED_TO_MIS  

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

  ViewVC Help
Powered by ViewVC 1.1.26