/[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 1881 by artak, Tue Oct 14 04:34:09 2008 UTC revision 1907 by phornby, Wed Oct 22 11:41:18 2008 UTC
# Line 42  Line 42 
42  #define IS_IN_MIS -3  #define IS_IN_MIS -3
43  #define IS_CONNECTED_TO_MIS -4  #define IS_CONNECTED_TO_MIS -4
44    
45  void Paso_Pattern_cop(Paso_SparseMatrix* A, index_t* mis_marker) {  
46    
47    void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
48    
49    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
50    dim_t i,j;    dim_t i,j;
51    double threshold=0.05;    double threshold=0.05;
52    index_t naib,iptr;    index_t iptr,*index,*where_p,diagptr;
53    bool_t flag;    bool_t flag;
54    dim_t n=A->pattern->numOutput;    dim_t n=A->pattern->numOutput;
55    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
# Line 58  void Paso_Pattern_cop(Paso_SparseMatrix* Line 60  void Paso_Pattern_cop(Paso_SparseMatrix*
60       /* is there any vertex available ?*/       /* is there any vertex available ?*/
61       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
62    
63             #pragma omp parallel for private(naib,i,iptr,flag) schedule(static)             #pragma omp parallel for private(i,iptr,flag) schedule(static)
64             for (i=0;i<n;++i) {             for (i=0;i<n;++i) {
65                if (mis_marker[i]==IS_AVAILABLE) {                if (mis_marker[i]==IS_AVAILABLE) {
66                   flag=IS_AVAILABLE;                   flag=IS_IN_MIS;
67                     diagptr=A->pattern->ptr[i];
68                     index=&(A->pattern->index[diagptr]);
69                     where_p=(index_t*)bsearch(&i,
70                                            index,
71                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
72                                            sizeof(index_t),
73                                            Paso_comparIndex);
74                    if (where_p==NULL) {
75                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
76                    } else {
77                        diagptr+=(index_t)(where_p-index);
78                    }
79                   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) {
80                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
81                       if (naib!=i && A->val[naib]<threshold*A->val[i]) {                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {
82                          flag=IS_IN_MIS;                          flag=IS_AVAILABLE;
83                          break;                          break;
84                       }                       }
85                   }                   }
86                   mis_marker[i]=flag;                   mis_marker[i]=flag;
87                }                  }
88             }              }
89              
90             #pragma omp parallel for private(naib,i,iptr) schedule(static)                #pragma omp parallel for private(i,iptr) 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_AVAILABLE) {
93                     diagptr=A->pattern->ptr[i];
94                     index=&(A->pattern->index[diagptr]);
95                     where_p=(index_t*)bsearch(&i,
96                                            index,
97                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
98                                            sizeof(index_t),
99                                            Paso_comparIndex);
100                    if (where_p==NULL) {
101                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
102                    } else {
103                        diagptr+=(index_t)(where_p-index);
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                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
107                        if (naib!=i && mis_marker[naib]==IS_IN_MIS && A->val[iptr]/A->val[i]>=-threshold)                       if (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){
108                           mis_marker[i]=IS_IN_MIS;                           mis_marker[i]=IS_IN_MIS;
109                       else                       }
110                         else {
111                           mis_marker[i]=IS_CONNECTED_TO_MIS;                           mis_marker[i]=IS_CONNECTED_TO_MIS;
112                         }
113                   }                   }
114                }                 }
115             }                }
116       }          }
117         /* swap to TRUE/FALSE in mis_marker */
118         #pragma omp parallel for private(i) schedule(static)
119         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
120    }
121    
122    /*
123     *
124     * Return a strength of connection mask using the classical
125     * strength of connection measure by Ruge and Stuben.
126     *
127     * Specifically, an off-diagonal entry A[i.j] is a strong
128     * connection if:
129     *  
130     *      -A[i,j] >= theta * max( -A[i,k] )   where k != i
131     *
132     * Otherwise, the connection is weak.
133     *  
134     */  
135    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);
138      dim_t i,j;
139      index_t iptr;
140      double threshold,min_offdiagonal;
141      bool_t flag;
142      dim_t n=A->pattern->numOutput;
143      if (A->pattern->type & PATTERN_FORMAT_SYM) {
144        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
145        return;
146      }
147    
148    /* is there any vertex available ?*/
149         if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
150    
151         #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
152          for (i=0;i<n;++i) {
153            min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
154            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
155                if(A->pattern->index[iptr] != i){
156                    min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
157                }
158            }
159    
160            threshold = theta*min_offdiagonal;
161            mis_marker[i]=IS_CONNECTED_TO_MIS;
162            #pragma omp parallel for private(iptr) schedule(static)
163            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
164                if(-1.*(A->val[iptr-index_offset]) <= threshold){
165                        mis_marker[i]=IS_IN_MIS;
166                }
167            }
168          }
169        }
170       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
171       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
172       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);
Legend:
Removed from v.1881  
changed lines
  Added in v.1907
  ViewVC Help
Powered by ViewVC 1.1.26