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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2728 - (hide 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 jfenwick 1851
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 jfenwick 1851 * 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 phornby 1914 #include "PasoUtil.h"
32 phornby 1913 #include "Pattern_coupling.h"
33 artak 2381 #include <limits.h>
34 jfenwick 1851
35    
36     /***************************************************************/
37    
38     #define IS_AVAILABLE -1
39 artak 2726 #define IS_IN_F -3 /* in F (strong) */
40     #define IS_IN_C -4 /* in C (weak) */
41 jfenwick 1851
42 artak 2652 void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43 jfenwick 1851
44 artak 1881 dim_t i,j;
45 artak 2442 /*double sum;*/
46 artak 2107 index_t iptr,*index,*where_p,*diagptr;
47     bool_t passed=FALSE;
48 artak 1954 dim_t n=A->numRows;
49 artak 2107 diagptr=MEMALLOC(n,index_t);
50 artak 2122
51 jfenwick 1851 if (A->pattern->type & PATTERN_FORMAT_SYM) {
52 artak 2159 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
53 jfenwick 1851 return;
54     }
55    
56 artak 2107 #pragma omp parallel for private(i) schedule(static)
57     for (i=0;i<n;++i)
58     if(mis_marker[i]==IS_AVAILABLE)
59 artak 2726 mis_marker[i]=IS_IN_C;
60 jfenwick 1851
61 artak 2726 /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
62 artak 2107 for (i=0;i<n;++i) {
63     diagptr[i]=A->pattern->ptr[i];
64 artak 2686 index=&(A->pattern->index[A->pattern->ptr[i]]);
65 artak 2107 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 artak 2159 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72 artak 2107 } 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 artak 2726 if (mis_marker[i]==IS_IN_C) {
82 artak 2107 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 artak 2726 mis_marker[j]=IS_IN_F;
86 artak 2107 }
87     }
88     }
89     }
90    
91    
92    
93     /*This loop cannot be parallelized, as order matters here.*/
94     for (i=0;i<n;i++) {
95 artak 2726 if (mis_marker[i]==IS_IN_F) {
96 artak 2442 passed=TRUE;
97 artak 2107 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98     j=A->pattern->index[iptr];
99 artak 2726 if (mis_marker[j]==IS_IN_C) {
100 artak 2107 if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101     passed=TRUE;
102 artak 1902 }
103 artak 2107 else {
104     passed=FALSE;
105     break;
106 artak 1902 }
107 artak 2307 }
108 artak 2442 }
109 artak 2726 if (passed) mis_marker[i]=IS_IN_C;
110 artak 1902 }
111 artak 2107 }
112 artak 2122
113 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
114     #pragma omp parallel for private(i) schedule(static)
115 artak 2726 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
116 artak 2107
117     MEMFREE(diagptr);
118 jfenwick 1851 }
119 artak 1890
120 artak 2726
121 artak 2381 /*
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 artak 2686 double threshold,max_offdiagonal;
130 artak 2381
131     Paso_Pattern *out=NULL;
132    
133     Paso_IndexList* index_list=NULL;
134 artak 1890
135 artak 2381 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 artak 2699 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
151 artak 2381 for (i=0;i<n;++i) {
152     if(mis_marker[i]==IS_AVAILABLE) {
153 artak 2686 max_offdiagonal = DBL_MIN;
154 artak 2381 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
155     if(A->pattern->index[iptr] != i){
156 artak 2686 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
157 artak 2381 }
158     }
159 artak 2475
160 artak 2686 threshold = theta*max_offdiagonal;
161 artak 2381 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
162     j=A->pattern->index[iptr];
163 artak 2686 if((-A->val[iptr])>=threshold) {
164 artak 2381 Paso_IndexList_insertIndex(&(index_list[i]),j);
165 artak 2686 Paso_IndexList_insertIndex(&(index_list[j]),i);
166 artak 2381 }
167     }
168     }
169     }
170    
171 artak 2475
172 artak 2381 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 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
182     Paso_Pattern_greedy(out,mis_marker);
183 gross 2551 Paso_Pattern_free(out);
184 artak 2381 }
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 artak 2447 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
211 artak 2381 return;
212     }
213    
214    
215 artak 2552 #pragma omp parallel for private(i,iptr,diag) schedule(static)
216 artak 2381 for (i=0;i<n;++i) {
217     diag = 0;
218     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
219 artak 2447 if(A->pattern->index[iptr] == i){
220 artak 2381 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 artak 2662 if((val*val)>=(eps_Aii*diags[j])) {
235 artak 2381 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 artak 2475
253 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
254     Paso_Pattern_greedy(out,mis_marker);
255 artak 2555 Paso_Pattern_free(out);
256 artak 2381
257 artak 2450 }
258 artak 2381
259 artak 2450 /* 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 artak 2726 mis_marker[i]=IS_IN_C;
277 artak 2450
278    
279     for (i=0;i<n;++i) {
280 artak 2726 if (mis_marker[i]==IS_IN_C) {
281 artak 2450 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
282     j=pattern->index[iptr];
283 artak 2726 mis_marker[j]=IS_IN_F;
284 artak 2450 }
285     }
286     }
287    
288    
289    
290     for (i=0;i<n;i++) {
291 artak 2726 if (mis_marker[i]==IS_IN_F) {
292 artak 2450 passed=TRUE;
293     for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
294     j=pattern->index[iptr];
295 artak 2726 if (mis_marker[j]==IS_IN_F) {
296 artak 2450 passed=TRUE;
297     }
298     else {
299     passed=FALSE;
300     break;
301     }
302     }
303 artak 2726 if (passed) mis_marker[i]=IS_IN_C;
304 artak 2450 }
305     }
306    
307     /* swap to TRUE/FALSE in mis_marker */
308     #pragma omp parallel for private(i) schedule(static)
309 artak 2726 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
310 artak 2450
311 artak 2381 }
312    
313    
314 artak 2475 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
315 artak 2381
316 artak 2475 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 artak 2381
325 artak 2475
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 artak 2726 mis_marker[i]=IS_IN_F;
340 artak 2475
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 artak 2726 if (mis_marker[i]==IS_IN_F) {
347 artak 2475 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
348     j=pattern->index[iptr];
349     if (colorOf[j]<color)
350 artak 2726 mis_marker[j]=IS_IN_C;
351 artak 2475 }
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 artak 2726 if (mis_marker[i]==IS_IN_C) {
364 artak 2475 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 artak 2726 if (mis_marker[j]==IS_IN_C) {
369 artak 2475 passed=TRUE;
370     }
371     else {
372     passed=FALSE;
373     /*break;*/
374     }
375     }
376     }
377 artak 2726 if (passed) mis_marker[i]=IS_IN_F;
378 artak 2475 }
379 artak 2704
380 artak 2475 }
381     }
382     }
383    
384     /* swap to TRUE/FALSE in mis_marker */
385     #pragma omp parallel for private(i) schedule(static)
386 artak 2726 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
387 artak 2475
388     MEMFREE(colorOf);
389     }
390    
391 artak 2726 /*For testing */
392     void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
393    
394 artak 2728 dim_t i,j,k;
395 artak 2726 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 artak 2728 j=0;
412 artak 2726
413     for (i=0;i<n;++i) {
414 artak 2728 rsum=0;
415 artak 2726 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 jfenwick 1851 #undef IS_AVAILABLE
487 artak 2726 #undef IS_IN_F
488     #undef IS_IN_C
489 artak 2686
490    
491    
492    
493    
494    
495    
496    
497    
498    

  ViewVC Help
Powered by ViewVC 1.1.26