/[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 1883 by artak, Tue Oct 14 04:46:27 2008 UTC revision 1913 by phornby, Thu Oct 23 08:27:33 2008 UTC
# Line 28  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
31    /*
32  #include "mpi_C.h"  #include "mpi_C.h"
33  #include "Paso.h"  #include "Paso.h"
34  #include "PasoUtil.h"  #include "PasoUtil.h"
35  #include "Pattern.h"  #include "Pattern.h"
36  #include "Solver.h"  */
37    #include "Pattern_coupling.h"
38    
39    
40  /***************************************************************/  /***************************************************************/
# Line 42  Line 44 
44  #define IS_IN_MIS -3  #define IS_IN_MIS -3
45  #define IS_CONNECTED_TO_MIS -4  #define IS_CONNECTED_TO_MIS -4
46    
47    
48    
49  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
50    
51    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
52    dim_t i,j;    dim_t i,j;
53    double threshold=0.05;    double threshold=0.05;
54    index_t naib,iptr;    index_t iptr,*index,*where_p,diagptr;
55    bool_t flag;    bool_t flag;
56    dim_t n=A->pattern->numOutput;    dim_t n=A->pattern->numOutput;
57    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
# Line 58  void Paso_Pattern_coup(Paso_SparseMatrix Line 62  void Paso_Pattern_coup(Paso_SparseMatrix
62       /* is there any vertex available ?*/       /* is there any vertex available ?*/
63       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {       while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
64    
65             #pragma omp parallel for private(naib,i,iptr,flag) schedule(static)             #pragma omp parallel for private(i,iptr,flag) schedule(static)
66             for (i=0;i<n;++i) {             for (i=0;i<n;++i) {
67                if (mis_marker[i]==IS_AVAILABLE) {                if (mis_marker[i]==IS_AVAILABLE) {
68                   flag=IS_AVAILABLE;                   flag=IS_IN_MIS;
69                     diagptr=A->pattern->ptr[i];
70                     index=&(A->pattern->index[diagptr]);
71                     where_p=(index_t*)bsearch(&i,
72                                            index,
73                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
74                                            sizeof(index_t),
75                                            Paso_comparIndex);
76                    if (where_p==NULL) {
77                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
78                    } else {
79                        diagptr+=(index_t)(where_p-index);
80                    }
81                   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) {
82                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
83                       if (naib!=i && A->val[iptr]>=threshold*A->val[i]) {                       if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {
84                          flag=IS_IN_MIS;                          flag=IS_AVAILABLE;
85                          break;                          break;
86                       }                       }
87                   }                   }
88                   mis_marker[i]=flag;                   mis_marker[i]=flag;
89                }                  }
90             }              }
91              
92             #pragma omp parallel for private(naib,i,iptr) schedule(static)                #pragma omp parallel for private(i,iptr) schedule(static)
93             for (i=0;i<n;i++) {                for (i=0;i<n;i++) {
94                if (mis_marker[i]==IS_AVAILABLE) {                 if (mis_marker[i]==IS_AVAILABLE) {
95                     diagptr=A->pattern->ptr[i];
96                     index=&(A->pattern->index[diagptr]);
97                     where_p=(index_t*)bsearch(&i,
98                                            index,
99                                            A->pattern->ptr[i + 1]-A->pattern->ptr[i],
100                                            sizeof(index_t),
101                                            Paso_comparIndex);
102                    if (where_p==NULL) {
103                        Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
104                    } else {
105                        diagptr+=(index_t)(where_p-index);
106                    }
107                   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) {
108                       naib=A->pattern->index[iptr]-index_offset;                       j=A->pattern->index[iptr]-index_offset;
109                       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){
110                           mis_marker[i]=IS_IN_MIS;                           mis_marker[i]=IS_IN_MIS;
111                       }                       }
112                       else {                       else {
113                           mis_marker[i]=IS_CONNECTED_TO_MIS;                           mis_marker[i]=IS_CONNECTED_TO_MIS;
114                       }                       }
115                   }                   }
116                }                 }
117             }                }
118       }          }
119         /* swap to TRUE/FALSE in mis_marker */
120         #pragma omp parallel for private(i) schedule(static)
121         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
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,j;
141      index_t iptr;
142      double threshold,min_offdiagonal;
143      bool_t flag;
144      dim_t n=A->pattern->numOutput;
145      if (A->pattern->type & PATTERN_FORMAT_SYM) {
146        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
147        return;
148      }
149    
150    /* is there any vertex available ?*/
151         if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
152    
153         #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
154          for (i=0;i<n;++i) {
155            min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
156            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
157                if(A->pattern->index[iptr] != i){
158                    min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
159                }
160            }
161    
162            threshold = theta*min_offdiagonal;
163            mis_marker[i]=IS_CONNECTED_TO_MIS;
164            #pragma omp parallel for private(iptr) schedule(static)
165            for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
166                if(-1.*(A->val[iptr-index_offset]) <= threshold){
167                        mis_marker[i]=IS_IN_MIS;
168                }
169            }
170          }
171        }
172       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
173       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
174       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.1883  
changed lines
  Added in v.1913
  ViewVC Help
Powered by ViewVC 1.1.26