/[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 2381 - (show annotations)
Tue Apr 14 03:46:59 2009 UTC (10 years, 5 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
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 #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_SET -3
40 #define IS_REMOVED -4
41
42 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43
44 dim_t i,j;
45 /*double threshold=0.05;*/
46 double sum;
47 index_t iptr,*index,*where_p,*diagptr;
48 bool_t passed=FALSE;
49 dim_t n=A->numRows;
50 diagptr=MEMALLOC(n,index_t);
51
52 if (A->pattern->type & PATTERN_FORMAT_SYM) {
53 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
54 return;
55 }
56
57 #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
62 #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 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73 } 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 }
103 else {
104 passed=FALSE;
105 break;
106 }
107 }
108 }
109 if (passed) {
110 mis_marker[i]=IS_IN_SET;
111 passed=FALSE;
112 }
113 }
114 }
115
116 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
118 #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119 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 if(ABS(sum)<1.e-25)
128 mis_marker[i]=IS_IN_SET;
129 }
130 }
131
132
133 /* swap to TRUE/FALSE in mis_marker */
134 #pragma omp parallel for private(i) schedule(static)
135 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
136
137 MEMFREE(diagptr);
138 }
139
140 /*
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
154 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 #undef IS_AVAILABLE
283 #undef IS_IN_SET
284 #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26