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