/[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 2475 - (show annotations)
Wed Jun 17 01:48:46 2009 UTC (10 years, 2 months ago) by artak
File MIME type: text/plain
File size: 12145 byte(s)
AMG now takes into account new SolverOptions class, namely one can select coarsening algorithms from python. Not all options are included into AMG, this will be implemented later on.
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 #include <limits.h>
34
35
36 /***************************************************************/
37
38 #define IS_AVAILABLE -1
39 #define IS_IN_SET -3 /* Week connection */
40 #define IS_REMOVED -4 /* strong */
41
42 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43
44 dim_t i,j;
45 /*double sum;*/
46 index_t iptr,*index,*where_p,*diagptr;
47 bool_t passed=FALSE;
48 dim_t n=A->numRows;
49 diagptr=MEMALLOC(n,index_t);
50
51 if (A->pattern->type & PATTERN_FORMAT_SYM) {
52 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
53 return;
54 }
55
56 #pragma omp parallel for private(i) schedule(static)
57 for (i=0;i<n;++i)
58 if(mis_marker[i]==IS_AVAILABLE)
59 mis_marker[i]=IS_IN_SET;
60
61 #pragma omp parallel for private(i,index,where_p) schedule(static)
62 for (i=0;i<n;++i) {
63 diagptr[i]=A->pattern->ptr[i];
64 index=&(A->pattern->index[diagptr[i]]);
65 where_p=(index_t*)bsearch(&i,
66 index,
67 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
68 sizeof(index_t),
69 Paso_comparIndex);
70 if (where_p==NULL) {
71 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72 } else {
73 diagptr[i]+=(index_t)(where_p-index);
74 }
75 }
76
77
78
79 /*This loop cannot be parallelized, as order matters here.*/
80 for (i=0;i<n;++i) {
81 if (mis_marker[i]==IS_IN_SET) {
82 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
83 j=A->pattern->index[iptr];
84 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85 mis_marker[j]=IS_REMOVED;
86 }
87 }
88 }
89 }
90
91
92
93 /*This loop cannot be parallelized, as order matters here.*/
94 for (i=0;i<n;i++) {
95 if (mis_marker[i]==IS_REMOVED) {
96 passed=TRUE;
97 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98 j=A->pattern->index[iptr];
99 if (mis_marker[j]==IS_IN_SET) {
100 if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101 passed=TRUE;
102 }
103 else {
104 passed=FALSE;
105 break;
106 }
107 }
108 }
109 if (passed) mis_marker[i]=IS_IN_SET;
110 }
111 }
112 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114 /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115 for (i=0;i<n;i++) {
116 if (mis_marker[i]==IS_REMOVED) {
117 sum=0;
118 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
119 j=A->pattern->index[iptr];
120 if (mis_marker[j]==IS_REMOVED)
121 sum+=A->val[iptr];
122 }
123 if(ABS(sum)<1.e-25)
124 mis_marker[i]=IS_IN_SET;
125 }
126 }
127 */
128
129 /* swap to TRUE/FALSE in mis_marker */
130 #pragma omp parallel for private(i) schedule(static)
131 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
132
133 MEMFREE(diagptr);
134 }
135
136 /*
137 * Ruge-Stueben strength of connection mask.
138 *
139 */
140 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
141 {
142 dim_t i,n,j;
143 index_t iptr;
144 double threshold,min_offdiagonal;
145
146 Paso_Pattern *out=NULL;
147
148 Paso_IndexList* index_list=NULL;
149
150 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
151 if (! Paso_checkPtr(index_list)) {
152 #pragma omp parallel for private(i) schedule(static)
153 for(i=0;i<A->pattern->numOutput;++i) {
154 index_list[i].extension=NULL;
155 index_list[i].n=0;
156 }
157 }
158
159
160 n=A->numRows;
161 if (A->pattern->type & PATTERN_FORMAT_SYM) {
162 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
163 return;
164 }
165 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold,j) schedule(static)
166 for (i=0;i<n;++i) {
167 if(mis_marker[i]==IS_AVAILABLE) {
168 min_offdiagonal = DBL_MAX;
169 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170 if(A->pattern->index[iptr] != i){
171 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
172 }
173 }
174
175 threshold = theta*min_offdiagonal;
176 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
177 j=A->pattern->index[iptr];
178 if(A->val[iptr]<=threshold) {
179 if(j!=i) {
180 Paso_IndexList_insertIndex(&(index_list[i]),j);
181 Paso_IndexList_insertIndex(&(index_list[j]),i);
182 }
183 }
184 }
185 }
186 }
187
188
189 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
190
191 /* clean up */
192 if (index_list!=NULL) {
193 #pragma omp parallel for private(i) schedule(static)
194 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
195 }
196 TMPMEMFREE(index_list);
197
198 /*Paso_Pattern_mis(out,mis_marker);*/
199 Paso_Pattern_greedy(out,mis_marker);
200 }
201
202 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
203 {
204 dim_t i,j,n;
205 index_t iptr;
206 double diag,eps_Aii,val;
207 double* diags;
208
209
210 Paso_Pattern *out=NULL;
211 Paso_IndexList* index_list=NULL;
212
213 n=A->numRows;
214 diags=MEMALLOC(n,double);
215
216 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
217 if (! Paso_checkPtr(index_list)) {
218 #pragma omp parallel for private(i) schedule(static)
219 for(i=0;i<A->pattern->numOutput;++i) {
220 index_list[i].extension=NULL;
221 index_list[i].n=0;
222 }
223 }
224
225 if (A->pattern->type & PATTERN_FORMAT_SYM) {
226 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
227 return;
228 }
229
230
231 #pragma omp parallel for private(i,iptr,diag) schedule(static)
232 for (i=0;i<n;++i) {
233 diag = 0;
234 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
235 if(A->pattern->index[iptr] == i){
236 diag+=A->val[iptr];
237 }
238 }
239 diags[i]=ABS(diag);
240 }
241
242
243 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
244 for (i=0;i<n;++i) {
245 if (mis_marker[i]==IS_AVAILABLE) {
246 eps_Aii = theta*theta*diags[i];
247 val=0.;
248 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
249 j=A->pattern->index[iptr];
250 val=A->val[iptr];
251 if(j!= i) {
252 if(val*val>=eps_Aii * diags[j]) {
253 Paso_IndexList_insertIndex(&(index_list[i]),j);
254 }
255 }
256 }
257 }
258 }
259
260 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
261
262 /* clean up */
263 if (index_list!=NULL) {
264 #pragma omp parallel for private(i) schedule(static)
265 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
266 }
267
268 TMPMEMFREE(index_list);
269 MEMFREE(diags);
270
271
272 /*Paso_Pattern_mis(out,mis_marker);*/
273 Paso_Pattern_greedy(out,mis_marker);
274
275 }
276
277 /* Greedy algorithm */
278 void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
279
280 dim_t i,j;
281 /*double sum;*/
282 index_t iptr;
283 bool_t passed=FALSE;
284 dim_t n=pattern->numOutput;
285
286 if (pattern->type & PATTERN_FORMAT_SYM) {
287 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
288 return;
289 }
290
291 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
292 #pragma omp parallel for private(i) schedule(static)
293 for (i=0;i<n;++i)
294 if(mis_marker[i]==IS_AVAILABLE)
295 mis_marker[i]=IS_IN_SET;
296
297
298 for (i=0;i<n;++i) {
299 if (mis_marker[i]==IS_IN_SET) {
300 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
301 j=pattern->index[iptr];
302 mis_marker[j]=IS_REMOVED;
303 }
304 }
305 }
306
307
308
309 for (i=0;i<n;i++) {
310 if (mis_marker[i]==IS_REMOVED) {
311 passed=TRUE;
312 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
313 j=pattern->index[iptr];
314 if (mis_marker[j]==IS_REMOVED) {
315 passed=TRUE;
316 }
317 else {
318 passed=FALSE;
319 break;
320 }
321 }
322 }
323 if (passed) mis_marker[i]=IS_IN_SET;
324 }
325
326 /* swap to TRUE/FALSE in mis_marker */
327 #pragma omp parallel for private(i) schedule(static)
328 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
329
330 }
331
332
333 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
334
335 dim_t i,j;
336 /*double sum;*/
337 index_t iptr;
338 index_t num_colors;
339 index_t* colorOf;
340 register index_t color;
341 bool_t passed=FALSE;
342 dim_t n=pattern->numOutput;
343
344
345 colorOf=MEMALLOC(n,index_t);
346
347 if (pattern->type & PATTERN_FORMAT_SYM) {
348 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
349 return;
350 }
351
352 Paso_Pattern_color(pattern,&num_colors,colorOf);
353
354 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
355 #pragma omp parallel for private(i) schedule(static)
356 for (i=0;i<n;++i)
357 if(mis_marker[i]==IS_AVAILABLE)
358 mis_marker[i]=IS_IN_SET;
359
360 #pragma omp barrier
361 for (color=0;color<num_colors;++color) {
362 #pragma omp parallel for schedule(static) private(i,iptr,j)
363 for (i=0;i<n;++i) {
364 if (colorOf[i]==color) {
365 if (mis_marker[i]==IS_IN_SET) {
366 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
367 j=pattern->index[iptr];
368 if (colorOf[j]<color)
369 mis_marker[j]=IS_REMOVED;
370 }
371 }
372 }
373 }
374 }
375
376
377 #pragma omp barrier
378 for (color=0;color<num_colors;++color) {
379 #pragma omp parallel for schedule(static) private(i,iptr,j)
380 for (i=0;i<n;i++) {
381 if (colorOf[i]==color) {
382 if (mis_marker[i]==IS_REMOVED) {
383 passed=TRUE;
384 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
385 j=pattern->index[iptr];
386 if (colorOf[j]<color && passed) {
387 if (mis_marker[j]==IS_REMOVED) {
388 passed=TRUE;
389 }
390 else {
391 passed=FALSE;
392 /*break;*/
393 }
394 }
395 }
396 }
397 if (passed) mis_marker[i]=IS_IN_SET;
398 }
399 }
400 }
401
402 /* swap to TRUE/FALSE in mis_marker */
403 #pragma omp parallel for private(i) schedule(static)
404 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
405
406 MEMFREE(colorOf);
407 }
408
409
410
411 #undef IS_AVAILABLE
412 #undef IS_IN_SET
413 #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26