/[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 2475 - (hide annotations)
Wed Jun 17 01:48:46 2009 UTC (10 years, 8 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 jfenwick 1851
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 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 1931 void Paso_Pattern_coup(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     mis_marker[i]=IS_IN_SET;
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     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 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     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 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     if (mis_marker[j]==IS_IN_SET) {
100     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     if (passed) mis_marker[i]=IS_IN_SET;
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     double threshold,min_offdiagonal;
145    
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 2475 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold,j) schedule(static)
166 artak 2381 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 artak 2475
175 artak 2381 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 artak 2442 Paso_IndexList_insertIndex(&(index_list[j]),i);
182 artak 2381 }
183     }
184     }
185     }
186     }
187    
188 artak 2475
189 artak 2381 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 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
199     Paso_Pattern_greedy(out,mis_marker);
200 artak 2381 }
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 artak 2447 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
227 artak 2381 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 artak 2447 if(A->pattern->index[iptr] == i){
236 artak 2381 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 artak 2475
272 artak 2450 /*Paso_Pattern_mis(out,mis_marker);*/
273     Paso_Pattern_greedy(out,mis_marker);
274 artak 2381
275 artak 2450 }
276 artak 2381
277 artak 2450 /* 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 artak 2475 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
292 artak 2450 #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 artak 2381 }
331    
332    
333 artak 2475 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
334 artak 2381
335 artak 2475 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 artak 2381
344 artak 2475
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 jfenwick 1851 #undef IS_AVAILABLE
412 artak 1975 #undef IS_IN_SET
413     #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26