/[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 2699 - (show 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
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 /* Week connection */
40 #define IS_REMOVED -4 /* strong */
41
42 void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43
44 dim_t i,j;
45 /*double sum;*/
46 index_t iptr,*index,*where_p,*diagptr;
47 bool_t passed=FALSE;
48 dim_t n=A->numRows;
49 diagptr=MEMALLOC(n,index_t);
50
51 if (A->pattern->type & PATTERN_FORMAT_SYM) {
52 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
53 return;
54 }
55
56 #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_REMOVED;
60
61 #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[A->pattern->ptr[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 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72 } 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_REMOVED) {
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_IN_SET;
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_IN_SET) {
96 passed=TRUE;
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_REMOVED) {
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) mis_marker[i]=IS_REMOVED;
110 }
111 }
112 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114 /*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115 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 if(ABS(sum)<1.e-25)
124 mis_marker[i]=IS_IN_SET;
125 }
126 }
127 */
128
129 /* swap to TRUE/FALSE in mis_marker */
130 #pragma omp parallel for private(i) schedule(static)
131 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
132
133 MEMFREE(diagptr);
134 }
135
136 /*
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,max_offdiagonal;
145
146 Paso_Pattern *out=NULL;
147
148 Paso_IndexList* index_list=NULL;
149
150 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 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
166 for (i=0;i<n;++i) {
167 if(mis_marker[i]==IS_AVAILABLE) {
168 max_offdiagonal = DBL_MIN;
169 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170 if(A->pattern->index[iptr] != i){
171 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
172 }
173 }
174
175 threshold = theta*max_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 Paso_IndexList_insertIndex(&(index_list[i]),j);
180 Paso_IndexList_insertIndex(&(index_list[j]),i);
181 }
182 }
183 }
184 }
185
186
187 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 /*Paso_Pattern_mis(out,mis_marker);*/
197 Paso_Pattern_greedy(out,mis_marker);
198 Paso_Pattern_free(out);
199 }
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 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
226 return;
227 }
228
229
230 #pragma omp parallel for private(i,iptr,diag) schedule(static)
231 for (i=0;i<n;++i) {
232 diag = 0;
233 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
234 if(A->pattern->index[iptr] == i){
235 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 if((val*val)>=(eps_Aii*diags[j])) {
251 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
270 /*Paso_Pattern_mis(out,mis_marker);*/
271 Paso_Pattern_greedy(out,mis_marker);
272 Paso_Pattern_free(out);
273
274 }
275
276 /* 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 /* We do not need this loop if we set IS_REMOVED=IS_AVAILABLE. */
291 #pragma omp parallel for private(i) schedule(static)
292 for (i=0;i<n;++i)
293 if(mis_marker[i]==IS_AVAILABLE)
294 mis_marker[i]=IS_REMOVED;
295
296
297 for (i=0;i<n;++i) {
298 if (mis_marker[i]==IS_REMOVED) {
299 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
300 j=pattern->index[iptr];
301 mis_marker[j]=IS_IN_SET;
302 }
303 }
304 }
305
306
307
308 for (i=0;i<n;i++) {
309 if (mis_marker[i]==IS_IN_SET) {
310 passed=TRUE;
311 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
312 j=pattern->index[iptr];
313 if (mis_marker[j]==IS_IN_SET) {
314 passed=TRUE;
315 }
316 else {
317 passed=FALSE;
318 break;
319 }
320 }
321 if (passed) mis_marker[i]=IS_REMOVED;
322 }
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 }
330
331
332 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
333
334 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
343
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 #undef IS_AVAILABLE
409 #undef IS_IN_SET
410 #undef IS_REMOVED
411
412
413
414
415
416
417
418
419
420

  ViewVC Help
Powered by ViewVC 1.1.26