/[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 2381 - (hide annotations)
Tue Apr 14 03:46:59 2009 UTC (10 years, 10 months ago) by artak
File MIME type: text/plain
File size: 8369 byte(s)
Ruge Stuben and Aggregiation is added for testing performance of AMG solver.
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 1975 #define IS_IN_SET -3
40     #define IS_REMOVED -4
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 1931 /*double threshold=0.05;*/
46 artak 2122 double sum;
47 artak 2107 index_t iptr,*index,*where_p,*diagptr;
48     bool_t passed=FALSE;
49 artak 1954 dim_t n=A->numRows;
50 artak 2107 diagptr=MEMALLOC(n,index_t);
51 artak 2122
52 jfenwick 1851 if (A->pattern->type & PATTERN_FORMAT_SYM) {
53 artak 2159 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
54 jfenwick 1851 return;
55     }
56    
57 artak 2107 #pragma omp parallel for private(i) schedule(static)
58     for (i=0;i<n;++i)
59     if(mis_marker[i]==IS_AVAILABLE)
60     mis_marker[i]=IS_IN_SET;
61 jfenwick 1851
62 artak 2107 #pragma omp parallel for private(i,index,where_p) schedule(static)
63     for (i=0;i<n;++i) {
64     diagptr[i]=A->pattern->ptr[i];
65     index=&(A->pattern->index[diagptr[i]]);
66     where_p=(index_t*)bsearch(&i,
67     index,
68     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69     sizeof(index_t),
70     Paso_comparIndex);
71     if (where_p==NULL) {
72 artak 2159 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73 artak 2107 } else {
74     diagptr[i]+=(index_t)(where_p-index);
75     }
76     }
77    
78    
79    
80     /*This loop cannot be parallelized, as order matters here.*/
81     for (i=0;i<n;++i) {
82     if (mis_marker[i]==IS_IN_SET) {
83     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
84     j=A->pattern->index[iptr];
85     if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
86     mis_marker[j]=IS_REMOVED;
87     }
88     }
89     }
90     }
91    
92    
93    
94     /*This loop cannot be parallelized, as order matters here.*/
95     for (i=0;i<n;i++) {
96     if (mis_marker[i]==IS_REMOVED) {
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 artak 1902 }
103 artak 2107 else {
104     passed=FALSE;
105     break;
106 artak 1902 }
107 artak 2307 }
108 artak 1902 }
109 artak 2107 if (passed) {
110     mis_marker[i]=IS_IN_SET;
111     passed=FALSE;
112     }
113 artak 1902 }
114 artak 2107 }
115    
116 artak 2381 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117 artak 2122 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
118 caltinay 2309 #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119 artak 2122 for (i=0;i<n;i++) {
120     if (mis_marker[i]==IS_REMOVED) {
121     sum=0;
122     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
123     j=A->pattern->index[iptr];
124     if (mis_marker[j]==IS_REMOVED)
125     sum+=A->val[iptr];
126     }
127 artak 2381 if(ABS(sum)<1.e-25)
128 artak 2122 mis_marker[i]=IS_IN_SET;
129     }
130     }
131    
132    
133 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
134     #pragma omp parallel for private(i) schedule(static)
135 artak 2107 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
136    
137     MEMFREE(diagptr);
138 jfenwick 1851 }
139 artak 1890
140 artak 2381 /*
141     * Ruge-Stueben strength of connection mask.
142     *
143     */
144     void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
145     {
146     dim_t i,n,j;
147     index_t iptr;
148     double threshold,min_offdiagonal;
149    
150     Paso_Pattern *out=NULL;
151    
152     Paso_IndexList* index_list=NULL;
153 artak 1890
154 artak 2381 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
155     if (! Paso_checkPtr(index_list)) {
156     #pragma omp parallel for private(i) schedule(static)
157     for(i=0;i<A->pattern->numOutput;++i) {
158     index_list[i].extension=NULL;
159     index_list[i].n=0;
160     }
161     }
162    
163    
164     n=A->numRows;
165     if (A->pattern->type & PATTERN_FORMAT_SYM) {
166     Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
167     return;
168     }
169    
170     #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
171     for (i=0;i<n;++i) {
172     if(mis_marker[i]==IS_AVAILABLE) {
173     min_offdiagonal = DBL_MAX;
174     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
175     if(A->pattern->index[iptr] != i){
176     min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
177     }
178     }
179    
180     threshold = theta*min_offdiagonal;
181     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
182     j=A->pattern->index[iptr];
183     if(A->val[iptr]<=threshold) {
184     if(j!=i) {
185     Paso_IndexList_insertIndex(&(index_list[i]),j);
186     }
187     }
188     }
189     }
190     }
191    
192     out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
193    
194     /* clean up */
195     if (index_list!=NULL) {
196     #pragma omp parallel for private(i) schedule(static)
197     for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
198     }
199     TMPMEMFREE(index_list);
200    
201     Paso_Pattern_mis(out,mis_marker);
202     }
203    
204     void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
205     {
206     dim_t i,j,n;
207     index_t iptr;
208     double diag,eps_Aii,val;
209     double* diags;
210    
211    
212     Paso_Pattern *out=NULL;
213     Paso_IndexList* index_list=NULL;
214    
215     n=A->numRows;
216     diags=MEMALLOC(n,double);
217    
218     index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
219     if (! Paso_checkPtr(index_list)) {
220     #pragma omp parallel for private(i) schedule(static)
221     for(i=0;i<A->pattern->numOutput;++i) {
222     index_list[i].extension=NULL;
223     index_list[i].n=0;
224     }
225     }
226    
227     if (A->pattern->type & PATTERN_FORMAT_SYM) {
228     Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
229     return;
230     }
231    
232    
233     #pragma omp parallel for private(i,iptr,diag) schedule(static)
234     for (i=0;i<n;++i) {
235     diag = 0;
236     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
237     if(A->pattern->index[iptr] != i){
238     diag+=A->val[iptr];
239     }
240     }
241     diags[i]=ABS(diag);
242     }
243    
244    
245     #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
246     for (i=0;i<n;++i) {
247     if (mis_marker[i]==IS_AVAILABLE) {
248     eps_Aii = theta*theta*diags[i];
249     val=0.;
250     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
251     j=A->pattern->index[iptr];
252     val=A->val[iptr];
253     if(j!= i) {
254     if(val*val>=eps_Aii * diags[j]) {
255     Paso_IndexList_insertIndex(&(index_list[i]),j);
256     }
257     }
258     }
259     }
260     }
261    
262     out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
263    
264     /* clean up */
265     if (index_list!=NULL) {
266     #pragma omp parallel for private(i) schedule(static)
267     for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
268     }
269    
270     TMPMEMFREE(index_list);
271     MEMFREE(diags);
272    
273     Paso_Pattern_mis(out,mis_marker);
274    
275    
276     }
277    
278    
279    
280    
281    
282 jfenwick 1851 #undef IS_AVAILABLE
283 artak 1975 #undef IS_IN_SET
284     #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26