/[escript]/trunk/paso/src/Pattern_coupling.c
ViewVC logotype

Contents of /trunk/paso/src/Pattern_coupling.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1954 - (show annotations)
Fri Oct 31 03:22:34 2008 UTC (10 years, 9 months ago) by artak
File MIME type: text/plain
File size: 8987 byte(s)
condition is added to check wheter row sum is ZERO and some new coarsening algorithm is also added
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**********************************************************************/
16
17 /* Paso: Pattern: Paso_Pattern_coupling
18
19 searches for a maximal independent set MIS in the matrix pattern
20 vertices in the maximal independent set are marked in mis_marker
21 nodes to be considered are marked by -1 on the input in mis_marker
22
23 */
24 /**********************************************************************/
25
26 /* Copyrights by ACcESS Australia 2003,2004,2005 */
27 /* Author: artak@uq.edu.au */
28
29 /**************************************************************/
30
31 #include "PasoUtil.h"
32 #include "Pattern_coupling.h"
33
34
35 /***************************************************************/
36
37 #define IS_AVAILABLE -1
38 #define IS_IN_MIS_NOW -2
39 #define IS_IN_MIS -3
40 #define IS_CONNECTED_TO_MIS -4
41 #define ZERO 1.e-10
42
43
44 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);
47 dim_t i,j;
48 /*double threshold=0.05;*/
49 index_t iptr,*index,*where_p,diagptr;
50 bool_t fail;
51 dim_t n=A->numRows;
52 double sum;
53 if (A->pattern->type & PATTERN_FORMAT_SYM) {
54 Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
55 return;
56 }
57
58 /* is there any vertex available ?*/
59 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
60
61 #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
62 for (i=0;i<n;++i) {
63 if (mis_marker[i]==IS_AVAILABLE) {
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) {
77 j=A->pattern->index[iptr]-index_offset;
78 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
79 mis_marker[j]=IS_CONNECTED_TO_MIS;
80 }
81 }
82 }
83 }
84
85 #pragma omp parallel for private(i,sum,iptr) schedule(static)
86 for (i=0;i<n;++i)
87 if(mis_marker[i]==IS_AVAILABLE)
88 mis_marker[i]=IS_IN_MIS;
89
90 #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
91 for (i=0;i<n;i++) {
92 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {
93 fail=FALSE;
94 diagptr=A->pattern->ptr[i];
95 index=&(A->pattern->index[diagptr]);
96 where_p=(index_t*)bsearch(&i,
97 index,
98 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
99 sizeof(index_t),
100 Paso_comparIndex);
101 if (where_p==NULL) {
102 Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
103 } else {
104 diagptr+=(index_t)(where_p-index);
105 }
106 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
107 j=A->pattern->index[iptr]-index_offset;
108 if (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){
109 fail=TRUE;
110 #ifndef _OPENMP
111 break;
112 #endif
113 }
114 }
115 if(!fail) {
116 mis_marker[i]=IS_IN_MIS;
117 }
118 }
119 }
120
121 #pragma omp parallel for private(i,sum,iptr) schedule(static)
122 for (i=0;i<n;++i)
123 if(mis_marker[i]==IS_IN_MIS)
124 {
125 sum=0.;
126 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
127 sum+=A->val[iptr];
128 }
129 if(ABS(sum)<ZERO)
130 mis_marker[i]=IS_CONNECTED_TO_MIS;
131 }
132
133 }
134 /* swap to TRUE/FALSE in mis_marker */
135 #pragma omp parallel for private(i) schedule(static)
136 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
137 }
138
139 /*
140 *
141 * Return a strength of connection mask using the classical
142 * strength of connection measure by Ruge and Stuben.
143 *
144 * Specifically, an off-diagonal entry A[i.j] is a strong
145 * connection if:
146 *
147 * -A[i,j] >= theta * max( -A[i,k] ) where k != i
148 *
149 * Otherwise, the connection is weak.
150 *
151 */
152 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
153 {
154 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
155 dim_t i;
156 index_t iptr;
157 double threshold,min_offdiagonal;
158 dim_t n=A->numRows;
159 if (A->pattern->type & PATTERN_FORMAT_SYM) {
160 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
161 return;
162 }
163
164 /* is there any vertex available ?*/
165 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
166
167 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
168 for (i=0;i<n;++i) {
169 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
170 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
171 if(A->pattern->index[iptr-index_offset] != i){
172 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
173 }
174 }
175
176 threshold = theta*min_offdiagonal;
177 mis_marker[i]=IS_IN_MIS;
178 #pragma omp parallel for private(iptr) schedule(static)
179 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
180 if(A->val[iptr-index_offset] < threshold){
181 mis_marker[i]=IS_CONNECTED_TO_MIS;
182 #ifndef _OPENMP
183 break;
184 #endif
185 }
186 }
187 }
188 }
189 /* swap to TRUE/FALSE in mis_marker */
190 #pragma omp parallel for private(i) schedule(static)
191 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
192 }
193
194
195
196
197 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
198 {
199 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
200 dim_t i,j;
201 index_t iptr;
202 double diag,eps_Aii,val;
203 dim_t n=A->numRows;
204 double* diags=MEMALLOC(n,double);
205
206 if (A->pattern->type & PATTERN_FORMAT_SYM) {
207 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
208 return;
209 }
210
211 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
212 #pragma omp parallel for private(i,iptr,diag) schedule(static)
213 for (i=0;i<n;++i) {
214 diag = 0;
215 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
216 if(A->pattern->index[iptr-index_offset] != i){
217 diag+=A->val[iptr-index_offset];
218 }
219 }
220 diags[i]=ABS(diag);
221 }
222
223 #pragma omp parallel for private(i,iptr,diag) schedule(static)
224 for (i=0;i<n;++i) {
225 eps_Aii = theta*theta*diags[i];
226 mis_marker[i]=IS_CONNECTED_TO_MIS;
227
228 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
229 j=A->pattern->index[iptr-index_offset];
230 val=A->val[iptr-index_offset];
231 if(j!= i && val*val>=eps_Aii * diags[j]){
232 mis_marker[i]=IS_IN_MIS;
233 }
234 }
235 }
236 }
237 /* swap to TRUE/FALSE in mis_marker */
238 #pragma omp parallel for private(i) schedule(static)
239 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
240
241 MEMFREE(diags);
242 }
243
244
245
246
247
248
249
250
251
252
253
254
255
256 #undef IS_AVAILABLE
257 #undef IS_IN_MIS_NOW
258 #undef IS_IN_MIS
259 #undef IS_CONNECTED_TO_MIS
260 #undef ZERO

  ViewVC Help
Powered by ViewVC 1.1.26