/[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 1915 by phornby, Thu Oct 23 08:40:09 2008 UTC revision 1936 by artak, Mon Oct 27 23:30:57 2008 UTC
# Line 41  Line 41 
41    
42    
43    
44  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
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->pattern->numOutput;
52    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
53      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");
# Line 57  void Paso_Pattern_coup(Paso_SparseMatrix Line 57  void Paso_Pattern_coup(Paso_SparseMatrix
57       /* is there any vertex available ?*/       /* is there any vertex available ?*/
58       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
59    
60             #pragma omp parallel for private(i,iptr,flag) schedule(static)              #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
61             for (i=0;i<n;++i) {              for (i=0;i<n;++i) {
62                if (mis_marker[i]==IS_AVAILABLE) {                if (mis_marker[i]==IS_AVAILABLE) {
                  flag=IS_IN_MIS;  
63                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
64                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
65                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 75  void Paso_Pattern_coup(Paso_SparseMatrix Line 74  void Paso_Pattern_coup(Paso_SparseMatrix
74                  }                  }
75                   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) {
76                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
77                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {                       if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
78                          flag=IS_AVAILABLE;                          mis_marker[j]=IS_CONNECTED_TO_MIS;
                         break;  
79                       }                       }
80                   }                   }
                  mis_marker[i]=flag;  
81                  }                  }
82              }              }
83                
84                #pragma omp parallel for private(i) schedule(static)
85                for (i=0;i<n;++i)
86                    if(mis_marker[i]==IS_AVAILABLE)
87                        mis_marker[i]=IS_IN_MIS;
88                        
89                #pragma omp parallel for private(i,iptr) schedule(static)                #pragma omp parallel for private(i,iptr,fail,index,where_p,diagptr) schedule(static)
90                for (i=0;i<n;i++) {                for (i=0;i<n;i++) {
91                 if (mis_marker[i]==IS_AVAILABLE) {                 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {
92                     fail=FALSE;
93                   diagptr=A->pattern->ptr[i];                   diagptr=A->pattern->ptr[i];
94                   index=&(A->pattern->index[diagptr]);                   index=&(A->pattern->index[diagptr]);
95                   where_p=(index_t*)bsearch(&i,                   where_p=(index_t*)bsearch(&i,
# Line 101  void Paso_Pattern_coup(Paso_SparseMatrix Line 104  void Paso_Pattern_coup(Paso_SparseMatrix
104                  }                  }
105                   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) {
106                       j=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
107                       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){
108                           mis_marker[i]=IS_IN_MIS;                           fail=TRUE;
109                       }                           break;
                      else {  
                          mis_marker[i]=IS_CONNECTED_TO_MIS;  
110                       }                       }
111                   }                   }
112                     if(!fail)
113                        mis_marker[i]=IS_IN_MIS;
114                 }                 }
115                }                }
116          }          }
# Line 132  void Paso_Pattern_coup(Paso_SparseMatrix Line 135  void Paso_Pattern_coup(Paso_SparseMatrix
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  {  {
137    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
138    dim_t i,j;    dim_t i;
139    index_t iptr;    index_t iptr;
140    double threshold,min_offdiagonal;    double threshold,min_offdiagonal;
   bool_t flag;  
141    dim_t n=A->pattern->numOutput;    dim_t n=A->pattern->numOutput;
142    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
143      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");
# Line 155  void Paso_Pattern_RS(Paso_SparseMatrix* Line 157  void Paso_Pattern_RS(Paso_SparseMatrix*
157          }          }
158    
159          threshold = theta*min_offdiagonal;          threshold = theta*min_offdiagonal;
160          mis_marker[i]=IS_CONNECTED_TO_MIS;          mis_marker[i]=IS_IN_MIS;
161          #pragma omp parallel for private(iptr) schedule(static)          #pragma omp parallel for private(iptr) schedule(static)
162          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) {
163              if(-1.*(A->val[iptr-index_offset]) <= threshold){              if(A->val[iptr-index_offset] < threshold){
164                      mis_marker[i]=IS_IN_MIS;                      mis_marker[i]=IS_CONNECTED_TO_MIS;
165                        #ifndef _OPENMP    
166                         break;
167                        #endif
168              }              }
169          }          }
170        }        }

Legend:
Removed from v.1915  
changed lines
  Added in v.1936

  ViewVC Help
Powered by ViewVC 1.1.26