/[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 2699 - (hide annotations)
Wed Sep 30 05:43:20 2009 UTC (10 years, 4 months ago) by artak
File MIME type: text/plain
File size: 12168 byte(s)
Direct solver switched to UMFPACK, till problem with MKL solved. When all unknows are eliminated then we switch to Jacobi preconditiner.
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 2442 #define IS_IN_SET -3 /* Week connection */
40     #define IS_REMOVED -4 /* strong */
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 2686 mis_marker[i]=IS_REMOVED;
60 jfenwick 1851
61 artak 2107 #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 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 2686 if (mis_marker[i]==IS_REMOVED) {
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 2686 mis_marker[j]=IS_IN_SET;
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 2686 if (mis_marker[i]==IS_IN_SET) {
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 2686 if (mis_marker[j]==IS_REMOVED) {
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 2686 if (passed) mis_marker[i]=IS_REMOVED;
110 artak 1902 }
111 artak 2107 }
112 artak 2381 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113 artak 2122 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114 artak 2442 /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115 artak 2122 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 artak 2381 if(ABS(sum)<1.e-25)
124 artak 2122 mis_marker[i]=IS_IN_SET;
125     }
126     }
127 artak 2442 */
128 artak 2122
129 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
130     #pragma omp parallel for private(i) schedule(static)
131 artak 2107 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
132    
133     MEMFREE(diagptr);
134 jfenwick 1851 }
135 artak 1890
136 artak 2381 /*
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 artak 2686 double threshold,max_offdiagonal;
145 artak 2381
146     Paso_Pattern *out=NULL;
147    
148     Paso_IndexList* index_list=NULL;
149 artak 1890
150 artak 2381 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 artak 2699 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
166 artak 2381 for (i=0;i<n;++i) {
167     if(mis_marker[i]==IS_AVAILABLE) {
168 artak 2686 max_offdiagonal = DBL_MIN;
169 artak 2381 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170     if(A->pattern->index[iptr] != i){
171 artak 2686 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
172 artak 2381 }
173     }
174 artak 2475
175 artak 2686 threshold = theta*max_offdiagonal;
176 artak 2381 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
177     j=A->pattern->index[iptr];
178 artak 2686 if((-A->val[iptr])>=threshold) {
179 artak 2381 Paso_IndexList_insertIndex(&(index_list[i]),j);
180 artak 2686 Paso_IndexList_insertIndex(&(index_list[j]),i);
181 artak 2381 }
182     }
183     }
184     }
185    
186 artak 2475
187 artak 2381 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
188    
189     /* clean up */
190     if (index_list!=NULL) {
191     #pragma omp parallel for private(i) schedule(static)
192     for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
193     }
194     TMPMEMFREE(index_list);
195    
196 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
197     Paso_Pattern_greedy(out,mis_marker);
198 gross 2551 Paso_Pattern_free(out);
199 artak 2381 }
200    
201     void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
202     {
203     dim_t i,j,n;
204     index_t iptr;
205     double diag,eps_Aii,val;
206     double* diags;
207    
208    
209     Paso_Pattern *out=NULL;
210     Paso_IndexList* index_list=NULL;
211    
212     n=A->numRows;
213     diags=MEMALLOC(n,double);
214    
215     index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
216     if (! Paso_checkPtr(index_list)) {
217     #pragma omp parallel for private(i) schedule(static)
218     for(i=0;i<A->pattern->numOutput;++i) {
219     index_list[i].extension=NULL;
220     index_list[i].n=0;
221     }
222     }
223    
224     if (A->pattern->type & PATTERN_FORMAT_SYM) {
225 artak 2447 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
226 artak 2381 return;
227     }
228    
229    
230 artak 2552 #pragma omp parallel for private(i,iptr,diag) schedule(static)
231 artak 2381 for (i=0;i<n;++i) {
232     diag = 0;
233     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
234 artak 2447 if(A->pattern->index[iptr] == i){
235 artak 2381 diag+=A->val[iptr];
236     }
237     }
238     diags[i]=ABS(diag);
239     }
240    
241    
242     #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
243     for (i=0;i<n;++i) {
244     if (mis_marker[i]==IS_AVAILABLE) {
245     eps_Aii = theta*theta*diags[i];
246     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
247     j=A->pattern->index[iptr];
248     val=A->val[iptr];
249     if(j!= i) {
250 artak 2662 if((val*val)>=(eps_Aii*diags[j])) {
251 artak 2381 Paso_IndexList_insertIndex(&(index_list[i]),j);
252     }
253     }
254     }
255     }
256     }
257    
258     out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
259    
260     /* clean up */
261     if (index_list!=NULL) {
262     #pragma omp parallel for private(i) schedule(static)
263     for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
264     }
265    
266     TMPMEMFREE(index_list);
267     MEMFREE(diags);
268    
269 artak 2475
270 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
271     Paso_Pattern_greedy(out,mis_marker);
272 artak 2555 Paso_Pattern_free(out);
273 artak 2381
274 artak 2450 }
275 artak 2381
276 artak 2450 /* Greedy algorithm */
277     void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
278    
279     dim_t i,j;
280     /*double sum;*/
281     index_t iptr;
282     bool_t passed=FALSE;
283     dim_t n=pattern->numOutput;
284    
285     if (pattern->type & PATTERN_FORMAT_SYM) {
286     Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
287     return;
288     }
289    
290 artak 2699 /* We do not need this loop if we set IS_REMOVED=IS_AVAILABLE. */
291 artak 2450 #pragma omp parallel for private(i) schedule(static)
292     for (i=0;i<n;++i)
293     if(mis_marker[i]==IS_AVAILABLE)
294 artak 2686 mis_marker[i]=IS_REMOVED;
295 artak 2450
296    
297     for (i=0;i<n;++i) {
298 artak 2686 if (mis_marker[i]==IS_REMOVED) {
299 artak 2450 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
300     j=pattern->index[iptr];
301 artak 2686 mis_marker[j]=IS_IN_SET;
302 artak 2450 }
303     }
304     }
305    
306    
307    
308     for (i=0;i<n;i++) {
309 artak 2686 if (mis_marker[i]==IS_IN_SET) {
310 artak 2450 passed=TRUE;
311     for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
312     j=pattern->index[iptr];
313 artak 2686 if (mis_marker[j]==IS_IN_SET) {
314 artak 2450 passed=TRUE;
315     }
316     else {
317     passed=FALSE;
318     break;
319     }
320     }
321 artak 2699 if (passed) mis_marker[i]=IS_REMOVED;
322 artak 2450 }
323     }
324    
325     /* swap to TRUE/FALSE in mis_marker */
326     #pragma omp parallel for private(i) schedule(static)
327     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
328    
329 artak 2381 }
330    
331    
332 artak 2475 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
333 artak 2381
334 artak 2475 dim_t i,j;
335     /*double sum;*/
336     index_t iptr;
337     index_t num_colors;
338     index_t* colorOf;
339     register index_t color;
340     bool_t passed=FALSE;
341     dim_t n=pattern->numOutput;
342 artak 2381
343 artak 2475
344     colorOf=MEMALLOC(n,index_t);
345    
346     if (pattern->type & PATTERN_FORMAT_SYM) {
347     Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
348     return;
349     }
350    
351     Paso_Pattern_color(pattern,&num_colors,colorOf);
352    
353     /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
354     #pragma omp parallel for private(i) schedule(static)
355     for (i=0;i<n;++i)
356     if(mis_marker[i]==IS_AVAILABLE)
357     mis_marker[i]=IS_IN_SET;
358    
359     #pragma omp barrier
360     for (color=0;color<num_colors;++color) {
361     #pragma omp parallel for schedule(static) private(i,iptr,j)
362     for (i=0;i<n;++i) {
363     if (colorOf[i]==color) {
364     if (mis_marker[i]==IS_IN_SET) {
365     for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
366     j=pattern->index[iptr];
367     if (colorOf[j]<color)
368     mis_marker[j]=IS_REMOVED;
369     }
370     }
371     }
372     }
373     }
374    
375    
376     #pragma omp barrier
377     for (color=0;color<num_colors;++color) {
378     #pragma omp parallel for schedule(static) private(i,iptr,j)
379     for (i=0;i<n;i++) {
380     if (colorOf[i]==color) {
381     if (mis_marker[i]==IS_REMOVED) {
382     passed=TRUE;
383     for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
384     j=pattern->index[iptr];
385     if (colorOf[j]<color && passed) {
386     if (mis_marker[j]==IS_REMOVED) {
387     passed=TRUE;
388     }
389     else {
390     passed=FALSE;
391     /*break;*/
392     }
393     }
394     }
395     }
396     if (passed) mis_marker[i]=IS_IN_SET;
397     }
398     }
399     }
400    
401     /* swap to TRUE/FALSE in mis_marker */
402     #pragma omp parallel for private(i) schedule(static)
403     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
404    
405     MEMFREE(colorOf);
406     }
407    
408 jfenwick 1851 #undef IS_AVAILABLE
409 artak 1975 #undef IS_IN_SET
410     #undef IS_REMOVED
411 artak 2686
412    
413    
414    
415    
416    
417    
418    
419    
420    

  ViewVC Help
Powered by ViewVC 1.1.26