/[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 2728 - (show annotations)
Thu Oct 22 00:13:10 2009 UTC (10 years, 4 months ago) by artak
File MIME type: text/plain
File size: 13787 byte(s)
minor
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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_F -3 /* in F (strong) */
40 #define IS_IN_C -4 /* in C (weak) */
41
42 void Paso_Pattern_YS(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_C;
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[A->pattern->ptr[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_C) {
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_IN_F;
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_IN_F) {
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_C) {
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_C;
110 }
111 }
112
113 /* swap to TRUE/FALSE in mis_marker */
114 #pragma omp parallel for private(i) schedule(static)
115 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
116
117 MEMFREE(diagptr);
118 }
119
120
121 /*
122 * Ruge-Stueben strength of connection mask.
123 *
124 */
125 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
126 {
127 dim_t i,n,j;
128 index_t iptr;
129 double threshold,max_offdiagonal;
130
131 Paso_Pattern *out=NULL;
132
133 Paso_IndexList* index_list=NULL;
134
135 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
136 if (! Paso_checkPtr(index_list)) {
137 #pragma omp parallel for private(i) schedule(static)
138 for(i=0;i<A->pattern->numOutput;++i) {
139 index_list[i].extension=NULL;
140 index_list[i].n=0;
141 }
142 }
143
144
145 n=A->numRows;
146 if (A->pattern->type & PATTERN_FORMAT_SYM) {
147 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
148 return;
149 }
150 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
151 for (i=0;i<n;++i) {
152 if(mis_marker[i]==IS_AVAILABLE) {
153 max_offdiagonal = DBL_MIN;
154 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
155 if(A->pattern->index[iptr] != i){
156 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
157 }
158 }
159
160 threshold = theta*max_offdiagonal;
161 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
162 j=A->pattern->index[iptr];
163 if((-A->val[iptr])>=threshold) {
164 Paso_IndexList_insertIndex(&(index_list[i]),j);
165 Paso_IndexList_insertIndex(&(index_list[j]),i);
166 }
167 }
168 }
169 }
170
171
172 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
173
174 /* clean up */
175 if (index_list!=NULL) {
176 #pragma omp parallel for private(i) schedule(static)
177 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
178 }
179 TMPMEMFREE(index_list);
180
181 /*Paso_Pattern_mis(out,mis_marker);*/
182 Paso_Pattern_greedy(out,mis_marker);
183 Paso_Pattern_free(out);
184 }
185
186 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
187 {
188 dim_t i,j,n;
189 index_t iptr;
190 double diag,eps_Aii,val;
191 double* diags;
192
193
194 Paso_Pattern *out=NULL;
195 Paso_IndexList* index_list=NULL;
196
197 n=A->numRows;
198 diags=MEMALLOC(n,double);
199
200 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
201 if (! Paso_checkPtr(index_list)) {
202 #pragma omp parallel for private(i) schedule(static)
203 for(i=0;i<A->pattern->numOutput;++i) {
204 index_list[i].extension=NULL;
205 index_list[i].n=0;
206 }
207 }
208
209 if (A->pattern->type & PATTERN_FORMAT_SYM) {
210 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
211 return;
212 }
213
214
215 #pragma omp parallel for private(i,iptr,diag) schedule(static)
216 for (i=0;i<n;++i) {
217 diag = 0;
218 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
219 if(A->pattern->index[iptr] == i){
220 diag+=A->val[iptr];
221 }
222 }
223 diags[i]=ABS(diag);
224 }
225
226
227 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
228 for (i=0;i<n;++i) {
229 if (mis_marker[i]==IS_AVAILABLE) {
230 eps_Aii = theta*theta*diags[i];
231 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
232 j=A->pattern->index[iptr];
233 val=A->val[iptr];
234 if((val*val)>=(eps_Aii*diags[j])) {
235 Paso_IndexList_insertIndex(&(index_list[i]),j);
236 }
237 }
238 }
239 }
240
241 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
242
243 /* clean up */
244 if (index_list!=NULL) {
245 #pragma omp parallel for private(i) schedule(static)
246 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
247 }
248
249 TMPMEMFREE(index_list);
250 MEMFREE(diags);
251
252
253 /*Paso_Pattern_mis(out,mis_marker);*/
254 Paso_Pattern_greedy(out,mis_marker);
255 Paso_Pattern_free(out);
256
257 }
258
259 /* Greedy algorithm */
260 void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
261
262 dim_t i,j;
263 /*double sum;*/
264 index_t iptr;
265 bool_t passed=FALSE;
266 dim_t n=pattern->numOutput;
267
268 if (pattern->type & PATTERN_FORMAT_SYM) {
269 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
270 return;
271 }
272
273 #pragma omp parallel for private(i) schedule(static)
274 for (i=0;i<n;++i)
275 if(mis_marker[i]==IS_AVAILABLE)
276 mis_marker[i]=IS_IN_C;
277
278
279 for (i=0;i<n;++i) {
280 if (mis_marker[i]==IS_IN_C) {
281 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
282 j=pattern->index[iptr];
283 mis_marker[j]=IS_IN_F;
284 }
285 }
286 }
287
288
289
290 for (i=0;i<n;i++) {
291 if (mis_marker[i]==IS_IN_F) {
292 passed=TRUE;
293 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
294 j=pattern->index[iptr];
295 if (mis_marker[j]==IS_IN_F) {
296 passed=TRUE;
297 }
298 else {
299 passed=FALSE;
300 break;
301 }
302 }
303 if (passed) mis_marker[i]=IS_IN_C;
304 }
305 }
306
307 /* swap to TRUE/FALSE in mis_marker */
308 #pragma omp parallel for private(i) schedule(static)
309 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
310
311 }
312
313
314 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
315
316 dim_t i,j;
317 /*double sum;*/
318 index_t iptr;
319 index_t num_colors;
320 index_t* colorOf;
321 register index_t color;
322 bool_t passed=FALSE;
323 dim_t n=pattern->numOutput;
324
325
326 colorOf=MEMALLOC(n,index_t);
327
328 if (pattern->type & PATTERN_FORMAT_SYM) {
329 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
330 return;
331 }
332
333 Paso_Pattern_color(pattern,&num_colors,colorOf);
334
335 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
336 #pragma omp parallel for private(i) schedule(static)
337 for (i=0;i<n;++i)
338 if(mis_marker[i]==IS_AVAILABLE)
339 mis_marker[i]=IS_IN_F;
340
341 #pragma omp barrier
342 for (color=0;color<num_colors;++color) {
343 #pragma omp parallel for schedule(static) private(i,iptr,j)
344 for (i=0;i<n;++i) {
345 if (colorOf[i]==color) {
346 if (mis_marker[i]==IS_IN_F) {
347 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
348 j=pattern->index[iptr];
349 if (colorOf[j]<color)
350 mis_marker[j]=IS_IN_C;
351 }
352 }
353 }
354 }
355 }
356
357
358 #pragma omp barrier
359 for (color=0;color<num_colors;++color) {
360 #pragma omp parallel for schedule(static) private(i,iptr,j)
361 for (i=0;i<n;i++) {
362 if (colorOf[i]==color) {
363 if (mis_marker[i]==IS_IN_C) {
364 passed=TRUE;
365 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
366 j=pattern->index[iptr];
367 if (colorOf[j]<color && passed) {
368 if (mis_marker[j]==IS_IN_C) {
369 passed=TRUE;
370 }
371 else {
372 passed=FALSE;
373 /*break;*/
374 }
375 }
376 }
377 if (passed) mis_marker[i]=IS_IN_F;
378 }
379
380 }
381 }
382 }
383
384 /* swap to TRUE/FALSE in mis_marker */
385 #pragma omp parallel for private(i) schedule(static)
386 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
387
388 MEMFREE(colorOf);
389 }
390
391 /*For testing */
392 void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
393
394 dim_t i,j,k;
395 double *theta;
396 index_t iptr;
397 dim_t n=A->numRows;
398 double rsum,diag=0;
399 index_t *AvADJ;
400 theta=MEMALLOC(n,double);
401 AvADJ=MEMALLOC(n,index_t);
402
403
404
405
406 if (A->pattern->type & PATTERN_FORMAT_SYM) {
407 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
408 return;
409 }
410
411 j=0;
412
413 for (i=0;i<n;++i) {
414 rsum=0;
415 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
416 j=A->pattern->index[iptr];
417 if(j!=i) {
418 rsum+=ABS(A->val[iptr]);
419 }
420 else {
421 diag=ABS(A->val[iptr]);
422 }
423 }
424 theta[i]=diag/rsum;
425 if(theta[i]>threshold) {
426 mis_marker[i]=IS_IN_F;
427 }
428 }
429
430 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
431 k=0;
432 for (i=0;i<n;++i) {
433 if(mis_marker[i]==IS_AVAILABLE) {
434 if(k==0) {
435 j=i;
436 k++;
437 }
438 if(theta[j]>theta[i]) {
439 j=i;
440 }
441 }
442 }
443 mis_marker[j]=IS_IN_C;
444
445 for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
446 k=A->pattern->index[iptr];
447 if(mis_marker[k]==IS_AVAILABLE) {
448 AvADJ[k]=1;
449 }
450 else {
451 AvADJ[k]=-1;
452 }
453
454 }
455
456 for (i=0;i<n;++i) {
457 if(AvADJ[i]) {
458 rsum=0;
459 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
460 k=A->pattern->index[iptr];
461 if(k!=i && mis_marker[k]!=IS_IN_C ) {
462 rsum+=ABS(A->val[iptr]);
463 }
464 if(j==i) {
465 diag=ABS(A->val[iptr]);
466 }
467 }
468 theta[i]=diag/rsum;
469 if(theta[i]>threshold) {
470 mis_marker[i]=IS_IN_F;
471 }
472 }
473 }
474
475
476 }
477
478 /* swap to TRUE/FALSE in mis_marker */
479 #pragma omp parallel for private(i) schedule(static)
480 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
481
482 MEMFREE(AvADJ);
483 MEMFREE(theta);
484 }
485
486 #undef IS_AVAILABLE
487 #undef IS_IN_F
488 #undef IS_IN_C
489
490
491
492
493
494
495
496
497
498

  ViewVC Help
Powered by ViewVC 1.1.26