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

  ViewVC Help
Powered by ViewVC 1.1.26