/[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 2828 - (show annotations)
Tue Dec 22 01:24:40 2009 UTC (9 years, 7 months ago) by artak
File MIME type: text/plain
File size: 28987 byte(s)
Smoother for AMG now can be selected from python. Now only Jacobi and Gauss-Seidel are available as smoothers.
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_NOT_AVAILABLE -2
40 #define IS_IN_F -3 /* in F (strong) */
41 #define IS_IN_C -4 /* in C (weak) */
42
43
44 #define IS_UNDECIDED -1
45 #define IS_STRONG -2
46 #define IS_WEAK -3
47
48
49 #define IS_IN_FA -5 /* test */
50 #define IS_IN_FB -6 /* test */
51
52 void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
53
54 dim_t i,j;
55 /*double sum;*/
56 index_t iptr;
57 /*index_t *index,*where_p;*/
58 bool_t passed=FALSE;
59 dim_t n=A->numRows;
60 double *diags;
61 diags=MEMALLOC(n,double);
62
63 if (A->pattern->type & PATTERN_FORMAT_SYM) {
64 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
65 return;
66 }
67
68 #pragma omp parallel for private(i) schedule(static)
69 for (i=0;i<n;++i)
70 if(mis_marker[i]==IS_AVAILABLE)
71 mis_marker[i]=IS_IN_C;
72
73 #pragma omp parallel for private(i,j) schedule(static)
74 for (i=0;i<n;++i) {
75 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
76 j=A->pattern->index[iptr];
77 if(i==j) {
78 diags[i]=A->val[iptr];
79 }
80 }
81 }
82
83 /*This loop cannot be parallelized, as order matters here.*/
84 for (i=0;i<n;++i) {
85 if (mis_marker[i]==IS_IN_C) {
86 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
87 j=A->pattern->index[iptr];
88 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(diags[i])) {
89 mis_marker[j]=IS_IN_F;
90 }
91 }
92 }
93 }
94
95
96
97 /*This loop cannot be parallelized, as order matters here.*/
98 for (i=0;i<n;i++) {
99 if (mis_marker[i]==IS_IN_F) {
100 passed=TRUE;
101 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
102 j=A->pattern->index[iptr];
103 if (mis_marker[j]==IS_IN_C) {
104 if ((A->val[iptr]/diags[i])>=-threshold) {
105 passed=TRUE;
106 }
107 else {
108 passed=FALSE;
109 break;
110 }
111 }
112 }
113 if (passed) mis_marker[i]=IS_IN_C;
114 }
115 }
116
117 /* swap to TRUE/FALSE in mis_marker */
118 #pragma omp parallel for private(i) schedule(static)
119 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
120
121 MEMFREE(diags);
122 }
123
124
125 /*
126 * Ruge-Stueben strength of connection mask.
127 *
128 */
129 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
130 {
131 dim_t i,n,j;
132 index_t iptr;
133 double threshold,max_offdiagonal;
134
135 Paso_Pattern *out=NULL;
136
137 Paso_IndexList* index_list=NULL;
138
139 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
140 if (! Paso_checkPtr(index_list)) {
141 #pragma omp parallel for private(i) schedule(static)
142 for(i=0;i<A->pattern->numOutput;++i) {
143 index_list[i].extension=NULL;
144 index_list[i].n=0;
145 }
146 }
147
148
149 n=A->numRows;
150 if (A->pattern->type & PATTERN_FORMAT_SYM) {
151 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
152 return;
153 }
154 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
155 for (i=0;i<n;++i) {
156 if(mis_marker[i]==IS_AVAILABLE) {
157 max_offdiagonal = DBL_MIN;
158 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
159 if(A->pattern->index[iptr] != i){
160 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
161 }
162 }
163
164 threshold = theta*max_offdiagonal;
165 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
166 j=A->pattern->index[iptr];
167 if((-A->val[iptr])>=threshold) {
168 Paso_IndexList_insertIndex(&(index_list[i]),j);
169 /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
170 }
171 }
172 }
173 }
174
175
176 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
177
178 /* clean up */
179 if (index_list!=NULL) {
180 #pragma omp parallel for private(i) schedule(static)
181 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
182 }
183 TMPMEMFREE(index_list);
184
185 /*Paso_Pattern_mis(out,mis_marker);*/
186 Paso_Pattern_greedy(out,mis_marker);
187 Paso_Pattern_free(out);
188 }
189
190 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
191 {
192 dim_t i,j,n;
193 index_t iptr;
194 double diag,eps_Aii,val;
195 double* diags;
196
197
198 Paso_Pattern *out=NULL;
199 Paso_IndexList* index_list=NULL;
200
201 n=A->numRows;
202 diags=MEMALLOC(n,double);
203
204 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
205 if (! Paso_checkPtr(index_list)) {
206 #pragma omp parallel for private(i) schedule(static)
207 for(i=0;i<A->pattern->numOutput;++i) {
208 index_list[i].extension=NULL;
209 index_list[i].n=0;
210 }
211 }
212
213 if (A->pattern->type & PATTERN_FORMAT_SYM) {
214 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
215 return;
216 }
217
218
219 #pragma omp parallel for private(i,iptr,diag) schedule(static)
220 for (i=0;i<n;++i) {
221 diag = 0;
222 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
223 if(A->pattern->index[iptr] == i){
224 diag+=A->val[iptr];
225 }
226 }
227 diags[i]=ABS(diag);
228 }
229
230
231 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
232 for (i=0;i<n;++i) {
233 if (mis_marker[i]==IS_AVAILABLE) {
234 eps_Aii = theta*theta*diags[i];
235 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
236 j=A->pattern->index[iptr];
237 val=A->val[iptr];
238 if((val*val)>=(eps_Aii*diags[j])) {
239 Paso_IndexList_insertIndex(&(index_list[i]),j);
240 }
241 }
242 }
243 }
244
245 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
246
247 /* clean up */
248 if (index_list!=NULL) {
249 #pragma omp parallel for private(i) schedule(static)
250 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
251 }
252
253 TMPMEMFREE(index_list);
254 MEMFREE(diags);
255
256
257 /*Paso_Pattern_mis(out,mis_marker);*/
258 Paso_Pattern_greedy(out,mis_marker);
259 Paso_Pattern_free(out);
260
261 }
262
263 /* Greedy algorithm */
264 void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
265
266 dim_t i,j;
267 /*double sum;*/
268 index_t iptr;
269 bool_t passed=FALSE;
270 dim_t n=pattern->numOutput;
271
272 if (pattern->type & PATTERN_FORMAT_SYM) {
273 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
274 return;
275 }
276
277 #pragma omp parallel for private(i) schedule(static)
278 for (i=0;i<n;++i)
279 if(mis_marker[i]==IS_AVAILABLE)
280 mis_marker[i]=IS_IN_C;
281
282
283 for (i=0;i<n;++i) {
284 if (mis_marker[i]==IS_IN_C) {
285 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
286 j=pattern->index[iptr];
287 mis_marker[j]=IS_IN_F;
288 }
289 }
290 }
291
292
293 for (i=0;i<n;i++) {
294 if (mis_marker[i]==IS_IN_F) {
295 passed=TRUE;
296 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
297 j=pattern->index[iptr];
298 if (mis_marker[j]==IS_IN_F) {
299 passed=TRUE;
300 }
301 else {
302 passed=FALSE;
303 break;
304 }
305 }
306 if (passed) mis_marker[i]=IS_IN_C;
307 }
308 }
309
310 /* swap to TRUE/FALSE in mis_marker */
311 #pragma omp parallel for private(i) schedule(static)
312 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
313
314 }
315
316 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
317
318 dim_t i,j;
319 /*double sum;*/
320 index_t iptr;
321 index_t num_colors;
322 index_t* colorOf;
323 register index_t color;
324 bool_t passed=FALSE;
325 dim_t n=pattern->numOutput;
326
327
328 colorOf=MEMALLOC(n,index_t);
329
330 if (pattern->type & PATTERN_FORMAT_SYM) {
331 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
332 return;
333 }
334
335 Paso_Pattern_color(pattern,&num_colors,colorOf);
336
337 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
338 #pragma omp parallel for private(i) schedule(static)
339 for (i=0;i<n;++i)
340 if(mis_marker[i]==IS_AVAILABLE)
341 mis_marker[i]=IS_IN_F;
342
343 #pragma omp barrier
344 for (color=0;color<num_colors;++color) {
345 #pragma omp parallel for schedule(static) private(i,iptr,j)
346 for (i=0;i<n;++i) {
347 if (colorOf[i]==color) {
348 if (mis_marker[i]==IS_IN_F) {
349 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
350 j=pattern->index[iptr];
351 if (colorOf[j]<color)
352 mis_marker[j]=IS_IN_C;
353 }
354 }
355 }
356 }
357 }
358
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_C) {
366 passed=TRUE;
367 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
368 j=pattern->index[iptr];
369 if (colorOf[j]<color && passed) {
370 if (mis_marker[j]==IS_IN_C) {
371 passed=TRUE;
372 }
373 else {
374 passed=FALSE;
375 /*break;*/
376 }
377 }
378 }
379 if (passed) mis_marker[i]=IS_IN_F;
380 }
381
382 }
383 }
384 }
385
386 /* swap to TRUE/FALSE in mis_marker */
387 #pragma omp parallel for private(i) schedule(static)
388 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
389
390 MEMFREE(colorOf);
391 }
392
393 /*For testing */
394 void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
395
396 dim_t i,j=0,k;
397 double *theta;
398 index_t iptr;
399 dim_t n=A->numRows;
400 double rsum,diag=0;
401 index_t *AvADJ;
402 theta=MEMALLOC(n,double);
403 AvADJ=MEMALLOC(n,index_t);
404
405
406
407
408 if (A->pattern->type & PATTERN_FORMAT_SYM) {
409 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
410 return;
411 }
412
413
414 #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
415 for (i=0;i<n;++i) {
416 rsum=0;
417 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
418 j=A->pattern->index[iptr];
419 if(j!=i) {
420 rsum+=ABS(A->val[iptr]);
421 }
422 else {
423 diag=ABS(A->val[iptr]);
424 }
425 }
426 theta[i]=diag/rsum;
427 if(theta[i]>threshold) {
428 mis_marker[i]=IS_IN_F;
429 }
430 }
431
432 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
433 k=0;
434
435 for (i=0;i<n;++i) {
436 if(mis_marker[i]==IS_AVAILABLE) {
437 if(k==0) {
438 j=i;
439 k++;
440 }
441 if(theta[j]>theta[i]) {
442 j=i;
443 }
444 }
445 }
446 mis_marker[j]=IS_IN_C;
447
448 for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
449 k=A->pattern->index[iptr];
450 if(mis_marker[k]==IS_AVAILABLE) {
451 AvADJ[k]=1;
452 }
453 else {
454 AvADJ[k]=-1;
455 }
456
457 }
458
459 for (i=0;i<n;++i) {
460 if(AvADJ[i]) {
461 rsum=0;
462 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
463 k=A->pattern->index[iptr];
464 if(k!=i && mis_marker[k]!=IS_IN_C ) {
465 rsum+=ABS(A->val[iptr]);
466 }
467 if(j==i) {
468 diag=ABS(A->val[iptr]);
469 }
470 }
471 theta[i]=diag/rsum;
472 if(theta[i]>threshold) {
473 mis_marker[i]=IS_IN_F;
474 }
475 }
476 }
477
478
479 }
480
481 /* swap to TRUE/FALSE in mis_marker */
482 #pragma omp parallel for private(i) schedule(static)
483 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
484
485 MEMFREE(AvADJ);
486 MEMFREE(theta);
487 }
488
489
490 void Paso_Pattern_YS_plus(Paso_SparseMatrix* A, index_t* mis_marker, double alpha, double taw, double delta) {
491
492 dim_t i,j;
493 /*double sum;*/
494 index_t iptr,*index,*where_p,*diagptr;
495 double *rsum;
496 double sum;
497 dim_t n=A->numRows;
498 Paso_SparseMatrix* A_alpha;
499 diagptr=MEMALLOC(n,index_t);
500 rsum=MEMALLOC(n,double);
501
502 if (A->pattern->type & PATTERN_FORMAT_SYM) {
503 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
504 return;
505 }
506
507 A_alpha=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern,1,1, FALSE);
508 #pragma omp parallel for private(i) schedule(static)
509 for (i=0;i<A->len;++i) {
510 A_alpha->val[i]=A->val[i];
511 }
512
513
514 #pragma omp parallel for private(i) schedule(static)
515 for (i=0;i<n;++i)
516 if(mis_marker[i]==IS_AVAILABLE)
517 mis_marker[i]=IS_IN_C;
518
519 /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
520 for (i=0;i<n;++i) {
521 diagptr[i]=A->pattern->ptr[i];
522 index=&(A->pattern->index[A->pattern->ptr[i]]);
523 where_p=(index_t*)bsearch(&i,
524 index,
525 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
526 sizeof(index_t),
527 Paso_comparIndex);
528 if (where_p==NULL) {
529 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
530 } else {
531 diagptr[i]+=(index_t)(where_p-index);
532 }
533 }
534
535
536 j=0;
537 for (i=0;i<n;++i) {
538 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
539 j=A_alpha->pattern->index[iptr];
540 if(i==j) {
541 A_alpha->val[iptr]=0;
542 }
543 else {
544 if( !(ABS(A_alpha->val[iptr])<alpha*MIN(ABS(A_alpha->val[diagptr[i]]),ABS(A_alpha->val[diagptr[j]])) || A_alpha->val[iptr]*A_alpha->val[diagptr[i]]>0 || A_alpha->val[iptr]*A_alpha->val[diagptr[i]]<0) ) {
545 A_alpha->val[iptr]=0;
546 }
547 }
548
549 }
550 }
551
552
553 for (i=0;i<n;++i) {
554 rsum[i]=0;
555 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
556 rsum[i]+=A_alpha->val[iptr];
557 }
558 }
559
560 #pragma omp parallel for private(i,j) schedule(static)
561 for (i=0;i<n;++i) {
562 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
563 j=A_alpha->pattern->index[iptr];
564 if (i==j) {
565 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr]+rsum[i];
566 }
567 else {
568 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr];
569 }
570
571 }
572 }
573
574 /*This loop cannot be parallelized, as order matters here.*/
575 for (i=0;i<n;++i) {
576 if (mis_marker[i]==IS_IN_C) {
577 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
578 j=A->pattern->index[iptr];
579 if (j!=i && ABS(A->val[iptr])>=alpha*MIN(ABS(A->val[diagptr[i]]),ABS(A->val[diagptr[j]]))) {
580 mis_marker[j]=IS_IN_F;
581 }
582 }
583 }
584 }
585
586
587 /*This loop cannot be parallelized, as order matters here.*/
588 for (i=0;i<n;i++) {
589 if (mis_marker[i]==IS_IN_F) {
590 sum=0;
591 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
592 j=A_alpha->pattern->index[iptr];
593 if (mis_marker[j]==IS_IN_C) {
594 sum+=A_alpha->val[iptr];
595 }
596 }
597 if (ABS(sum)>=taw*ABS(A->val[diagptr[i]])) {
598 mis_marker[j]=IS_IN_FA;
599 }
600 }
601 }
602
603 /*This loop cannot be parallelized, as order matters here.*/
604 for (i=0;i<n;i++) {
605 if (mis_marker[i]!=IS_IN_C || mis_marker[i]!=IS_IN_FA) {
606 sum=0;
607 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
608 j=A_alpha->pattern->index[iptr];
609 if (mis_marker[j]==IS_IN_C || mis_marker[j]==IS_IN_FA) {
610 sum+=A_alpha->val[iptr];
611 }
612 }
613 if (ABS(sum)>=delta*ABS(A->val[diagptr[i]])) {
614 mis_marker[j]=IS_IN_FB;
615 }
616 else {
617 mis_marker[j]=IS_IN_C;
618 }
619
620 }
621 }
622
623
624 /* swap to TRUE/FALSE in mis_marker */
625 #pragma omp parallel for private(i) schedule(static)
626 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_FA || mis_marker[i]==IS_IN_FB);
627
628 MEMFREE(diagptr);
629 MEMFREE(rsum);
630 Paso_SparseMatrix_free(A_alpha);
631 }
632
633
634 void Paso_Pattern_Standard(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
635 {
636 dim_t i,n,j,k;
637 index_t iptr,jptr;
638 /*index_t *index,*where_p;*/
639 double threshold,max_offdiagonal;
640 dim_t *lambda; /*mesure of importance */
641 dim_t maxlambda=0;
642 index_t index_maxlambda=0;
643 double time0=0;
644 bool_t verbose=0;
645
646 Paso_Pattern *S=NULL;
647 Paso_Pattern *ST=NULL;
648 Paso_IndexList* index_list=NULL;
649
650 /*dim_t lk;*/
651
652 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
653 if (! Paso_checkPtr(index_list)) {
654 #pragma omp parallel for private(i) schedule(static)
655 for(i=0;i<A->pattern->numOutput;++i) {
656 index_list[i].extension=NULL;
657 index_list[i].n=0;
658 }
659 }
660
661
662 n=A->numRows;
663 if (A->pattern->type & PATTERN_FORMAT_SYM) {
664 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
665 return;
666 }
667
668 time0=Paso_timer();
669 k=0;
670 /*S_i={j \in N_i; i strongly coupled to j}*/
671 #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)
672 for (i=0;i<n;++i) {
673 if(mis_marker[i]==IS_AVAILABLE) {
674 max_offdiagonal = DBL_MIN;
675 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
676 j=A->pattern->index[iptr];
677 if( j != i){
678 if(A->val[iptr]<0) {
679 max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
680 }
681 }
682 }
683 threshold = theta*max_offdiagonal;
684 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
685 j=A->pattern->index[iptr];
686 if(ABS(A->val[iptr])>=threshold) {
687 Paso_IndexList_insertIndex(&(index_list[i]),j);
688 }
689 }
690 }
691 }
692
693 S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
694 ST=Paso_Pattern_getTranspose(S);
695
696
697 time0=Paso_timer()-time0;
698 if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
699
700 lambda=TMPMEMALLOC(n,dim_t);
701
702 #pragma omp parallel for private(i) schedule(static)
703 for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }
704
705 /*S_i={j \in N_i; i strongly coupled to j}*/
706
707 /*
708 #pragma omp parallel for private(i,iptr,lk) schedule(static)
709 for (i=0;i<n;++i) {
710 lk=0;
711 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
712 if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )
713 lk++;
714 }
715 #pragma omp critical
716 k+=lk;
717 if(k==0) {
718 mis_marker[i]=IS_IN_F;
719 }
720 }
721 */
722
723
724 k=0;
725 maxlambda=0;
726
727 time0=Paso_timer();
728
729 for (i=0;i<n;++i) {
730 if(mis_marker[i]==IS_AVAILABLE) {
731 lambda[i]=how_many(k,ST,FALSE); /* if every row is available then k and i are the same.*/
732 /*lambda[i]=how_many(i,S,TRUE);*/
733 /*printf("lambda[%d]=%d, k=%d ",i,lambda[i],k);*/
734 k++;
735 if(maxlambda<lambda[i]) {
736 maxlambda=lambda[i];
737 index_maxlambda=i;
738 }
739 }
740 }
741
742 k=0;
743 time0=Paso_timer()-time0;
744 if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);
745
746 time0=Paso_timer();
747
748 /*Paso_Pattern_getReport(n,mis_marker);*/
749
750 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
751
752 if(index_maxlambda<0) {
753 break;
754 }
755
756 i=index_maxlambda;
757 if(mis_marker[i]==IS_AVAILABLE) {
758 mis_marker[index_maxlambda]=IS_IN_C;
759 lambda[index_maxlambda]=IS_NOT_AVAILABLE;
760 for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
761 j=ST->index[iptr];
762 if(mis_marker[j]==IS_AVAILABLE) {
763 mis_marker[j]=IS_IN_F;
764 lambda[j]=IS_NOT_AVAILABLE;
765 for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
766 k=S->index[jptr];
767 if(mis_marker[k]==IS_AVAILABLE) {
768 lambda[k]++;
769 }
770 }
771 }
772 }
773 }
774
775 /* Used when transpose of S is not available */
776 /*
777 for (i=0;i<n;++i) {
778 if(mis_marker[i]==IS_AVAILABLE) {
779 if (i==index_maxlambda) {
780 mis_marker[index_maxlambda]=IS_IN_C;
781 lambda[index_maxlambda]=-1;
782 for (j=0;j<n;++j) {
783 if(mis_marker[j]==IS_AVAILABLE) {
784 index=&(S->index[S->ptr[j]]);
785 where_p=(index_t*)bsearch(&i,
786 index,
787 S->ptr[j + 1]-S->ptr[j],
788 sizeof(index_t),
789 Paso_comparIndex);
790 if (where_p!=NULL) {
791 mis_marker[j]=IS_IN_F;
792 lambda[j]=-1;
793 for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
794 k=S->index[iptr];
795 if(mis_marker[k]==IS_AVAILABLE) {
796 lambda[k]++;
797 }
798 }
799 }
800
801 }
802 }
803
804 }
805 }
806 }
807 */
808 index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);
809 }
810
811 time0=Paso_timer()-time0;
812 if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
813
814 /*Paso_Pattern_getReport(n,mis_marker);*/
815
816 #pragma omp parallel for private(i) schedule(static)
817 for (i=0;i<n;++i)
818 if(mis_marker[i]==IS_AVAILABLE) {
819 mis_marker[i]=IS_IN_F;
820 }
821
822 /*Paso_Pattern_getReport(n,mis_marker);*/
823
824 TMPMEMFREE(lambda);
825
826 /* clean up */
827 if (index_list!=NULL) {
828 #pragma omp parallel for private(i) schedule(static)
829 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
830 }
831
832 TMPMEMFREE(index_list);
833 Paso_Pattern_free(S);
834
835 /* swap to TRUE/FALSE in mis_marker */
836 #pragma omp parallel for private(i) schedule(static)
837 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
838
839 }
840
841 void Paso_Pattern_getReport(dim_t n,index_t* mis_marker) {
842
843 dim_t i,av,nav,f,c;
844
845 av=0;
846 nav=0;
847 f=0;
848 c=0;
849
850 for (i=0;i<n;i++) {
851 if(mis_marker[i]==IS_NOT_AVAILABLE) {
852 nav++;
853 }
854 else if (mis_marker[i]==IS_AVAILABLE) {
855 av++;
856 }
857 else if (mis_marker[i]==IS_IN_F) {
858 f++;
859 }
860 else if (mis_marker[i]==IS_IN_C) {
861 c++;
862 }
863 }
864
865 printf("Available=%d, Not Available = %d, In F = %d, In C = %d \n",av,nav,f,c);
866
867 }
868
869
870 /*Used in Paso_Pattern_RS_MI*/
871
872 /*Computes how many connections unknown i have in S.
873 bool_t transpose - TRUE if we want to compute how many strong connection of i in S^T, FALSE otherwise.
874 Note that if we first transpose S and then call method on S^T then "transpose" should be set to FALSE.
875 */
876 dim_t how_many(dim_t i,Paso_Pattern * S, bool_t transpose) {
877 dim_t j,n;
878 dim_t total,ltotal;
879 index_t *index,*where_p;
880
881 /*index_t iptr;*/
882 total=0;
883 ltotal=0;
884
885 n=S->numOutput;
886
887 if(transpose) {
888 #pragma omp parallel
889 {
890 ltotal=0;
891 #pragma omp for private(j,index,where_p,ltotal) schedule(static)
892 for (j=0;j<n;++j) {
893 index=&(S->index[S->ptr[j]]);
894 where_p=(index_t*)bsearch(&i,
895 index,
896 S->ptr[j + 1]-S->ptr[j],
897 sizeof(index_t),
898 Paso_comparIndex);
899 if (where_p!=NULL) {
900 ltotal++;
901 }
902 }
903 }
904 #pragma omp critical
905 {
906 total+=ltotal;
907 }
908
909 }
910 else {
911 total=S->ptr[i+1]-S->ptr[i];
912 /*#pragma omp parallel for private(iptr) schedule(static)*/
913 /*for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
914 if(S->index[iptr]!=i && mis_marker[S->index[iptr]]==IS_AVAILABLE)
915 total++;
916 }*/
917
918 }
919
920 if (total==0) total=IS_NOT_AVAILABLE;
921
922 return total;
923 }
924
925 dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {
926 dim_t i;
927 dim_t max=-1;
928 dim_t argmax=-1;
929
930 #ifdef _OPENMP
931 dim_t lmax=-1;
932 dim_t li=-1;
933 #pragma omp parallel private(i,lmax,li)
934 {
935 lmax=-1;
936 li=-1;
937 #pragma omp for schedule(static)
938 for (i=0;i<n;++i) {
939 if(lmax<lambda[i] && lambda[i]!=mask){
940 lmax=lambda[i];
941 li=i;
942 }
943 }
944 #pragma omp critical
945 {
946 if (max<=lmax) {
947 if(max==lmax && argmax>li)
948 {
949 argmax=li;
950 }
951 if (max < lmax)
952 {
953 max=lmax;
954 argmax=li;
955 }
956 }
957 }
958 }
959 #else
960 for (i=0;i<n;++i) {
961 if(max<lambda[i] && lambda[i]!=mask){
962 max=lambda[i];
963 argmax=i;
964 }
965 }
966 #endif
967
968 return argmax;
969 }
970
971
972 Paso_Pattern* Paso_Pattern_getTranspose(Paso_Pattern* P){
973
974 Paso_Pattern *outpattern=NULL;
975
976 Paso_IndexList* index_list=NULL;
977 dim_t C=P->numInput;
978 dim_t F=P->numOutput-C;
979 dim_t n=C+F;
980 dim_t i,j;
981 index_t iptr;
982
983 index_list=TMPMEMALLOC(C,Paso_IndexList);
984 if (! Paso_checkPtr(index_list)) {
985 #pragma omp parallel for private(i) schedule(static)
986 for(i=0;i<C;++i) {
987 index_list[i].extension=NULL;
988 index_list[i].n=0;
989 }
990 }
991
992
993 /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
994 for (i=0;i<n;++i) {
995 for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {
996 j=P->index[iptr];
997 Paso_IndexList_insertIndex(&(index_list[j]),i);
998 }
999 }
1000
1001 outpattern=Paso_IndexList_createPattern(0, C,index_list,0,n,0);
1002
1003 /* clean up */
1004 if (index_list!=NULL) {
1005 #pragma omp parallel for private(i) schedule(static)
1006 for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
1007 }
1008 TMPMEMFREE(index_list);
1009
1010 return outpattern;
1011 }
1012
1013
1014 #undef IS_AVAILABLE
1015 #undef IS_NOT_AVAILABLE
1016 #undef IS_IN_F
1017 #undef IS_IN_C
1018
1019 #undef IS_UNDECIDED
1020 #undef IS_STRONG
1021 #undef IS_WEAK
1022
1023 #undef IS_IN_FB
1024 #undef IS_IN_FB
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034

  ViewVC Help
Powered by ViewVC 1.1.26