/[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 2760 - (show annotations)
Thu Nov 19 05:22:45 2009 UTC (10 years, 3 months ago) by artak
File MIME type: text/plain
File size: 25150 byte(s)
The first version of the new AMG preconditioner. It need a lot of polishing for efficiency. Old AMG now called AMLI preconditioner.
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 #pragma omp parallel for private(i,j,k) schedule(static)
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 #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
465 for (i=0;i<n;++i) {
466 if(AvADJ[i]) {
467 rsum=0;
468 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
469 k=A->pattern->index[iptr];
470 if(k!=i && mis_marker[k]!=IS_IN_C ) {
471 rsum+=ABS(A->val[iptr]);
472 }
473 if(j==i) {
474 diag=ABS(A->val[iptr]);
475 }
476 }
477 theta[i]=diag/rsum;
478 if(theta[i]>threshold) {
479 mis_marker[i]=IS_IN_F;
480 }
481 }
482 }
483
484
485 }
486
487 /* swap to TRUE/FALSE in mis_marker */
488 #pragma omp parallel for private(i) schedule(static)
489 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
490
491 MEMFREE(AvADJ);
492 MEMFREE(theta);
493 }
494
495
496 void Paso_Pattern_YS_plus(Paso_SparseMatrix* A, index_t* mis_marker, double alpha, double taw, double delta) {
497
498 dim_t i,j;
499 /*double sum;*/
500 index_t iptr,*index,*where_p,*diagptr;
501 double *rsum;
502 double sum;
503 dim_t n=A->numRows;
504 Paso_SparseMatrix* A_alpha;
505 diagptr=MEMALLOC(n,index_t);
506 rsum=MEMALLOC(n,double);
507
508 if (A->pattern->type & PATTERN_FORMAT_SYM) {
509 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
510 return;
511 }
512
513 A_alpha=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern,1,1, FALSE);
514 #pragma omp parallel for private(i) schedule(static)
515 for (i=0;i<A->len;++i) {
516 A_alpha->val[i]=A->val[i];
517 }
518
519
520 #pragma omp parallel for private(i) schedule(static)
521 for (i=0;i<n;++i)
522 if(mis_marker[i]==IS_AVAILABLE)
523 mis_marker[i]=IS_IN_C;
524
525 /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
526 for (i=0;i<n;++i) {
527 diagptr[i]=A->pattern->ptr[i];
528 index=&(A->pattern->index[A->pattern->ptr[i]]);
529 where_p=(index_t*)bsearch(&i,
530 index,
531 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
532 sizeof(index_t),
533 Paso_comparIndex);
534 if (where_p==NULL) {
535 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
536 } else {
537 diagptr[i]+=(index_t)(where_p-index);
538 }
539 }
540
541
542 j=0;
543 for (i=0;i<n;++i) {
544 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
545 j=A_alpha->pattern->index[iptr];
546 if(i==j) {
547 A_alpha->val[iptr]=0;
548 }
549 else {
550 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) ) {
551 A_alpha->val[iptr]=0;
552 }
553 }
554
555 }
556 }
557
558
559 for (i=0;i<n;++i) {
560 rsum[i]=0;
561 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
562 rsum[i]+=A_alpha->val[iptr];
563 }
564 }
565
566 #pragma omp parallel for private(i,j) schedule(static)
567 for (i=0;i<n;++i) {
568 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
569 j=A_alpha->pattern->index[iptr];
570 if (i==j) {
571 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr]+rsum[i];
572 }
573 else {
574 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr];
575 }
576
577 }
578 }
579
580 /*This loop cannot be parallelized, as order matters here.*/
581 for (i=0;i<n;++i) {
582 if (mis_marker[i]==IS_IN_C) {
583 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
584 j=A->pattern->index[iptr];
585 if (j!=i && ABS(A->val[iptr])>=alpha*MIN(ABS(A->val[diagptr[i]]),ABS(A->val[diagptr[j]]))) {
586 mis_marker[j]=IS_IN_F;
587 }
588 }
589 }
590 }
591
592
593 /*This loop cannot be parallelized, as order matters here.*/
594 for (i=0;i<n;i++) {
595 if (mis_marker[i]==IS_IN_F) {
596 sum=0;
597 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
598 j=A_alpha->pattern->index[iptr];
599 if (mis_marker[j]==IS_IN_C) {
600 sum+=A_alpha->val[iptr];
601 }
602 }
603 if (ABS(sum)>=taw*ABS(A->val[diagptr[i]])) {
604 mis_marker[j]=IS_IN_FA;
605 }
606 }
607 }
608
609 /*This loop cannot be parallelized, as order matters here.*/
610 for (i=0;i<n;i++) {
611 if (mis_marker[i]!=IS_IN_C || mis_marker[i]!=IS_IN_FA) {
612 sum=0;
613 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
614 j=A_alpha->pattern->index[iptr];
615 if (mis_marker[j]==IS_IN_C || mis_marker[j]==IS_IN_FA) {
616 sum+=A_alpha->val[iptr];
617 }
618 }
619 if (ABS(sum)>=delta*ABS(A->val[diagptr[i]])) {
620 mis_marker[j]=IS_IN_FB;
621 }
622 else {
623 mis_marker[j]=IS_IN_C;
624 }
625
626 }
627 }
628
629
630 /* swap to TRUE/FALSE in mis_marker */
631 #pragma omp parallel for private(i) schedule(static)
632 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_FA || mis_marker[i]==IS_IN_FB);
633
634 MEMFREE(diagptr);
635 MEMFREE(rsum);
636 Paso_SparseMatrix_free(A_alpha);
637 }
638
639
640 void Paso_Pattern_RS_MI(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
641 {
642 dim_t i,n,j,k;
643 index_t iptr,*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
650 Paso_Pattern *S=NULL;
651 Paso_IndexList* index_list=NULL;
652
653 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
654 if (! Paso_checkPtr(index_list)) {
655 #pragma omp parallel for private(i) schedule(static)
656 for(i=0;i<A->pattern->numOutput;++i) {
657 index_list[i].extension=NULL;
658 index_list[i].n=0;
659 }
660 }
661
662
663 n=A->numRows;
664 if (A->pattern->type & PATTERN_FORMAT_SYM) {
665 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
666 return;
667 }
668
669 /*S_i={j \in N_i; i strongly coupled to j}*/
670 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
671 for (i=0;i<n;++i) {
672 if(mis_marker[i]==IS_AVAILABLE) {
673 max_offdiagonal = DBL_MIN;
674 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
675 if(A->pattern->index[iptr] != i){
676 if(A->val[iptr]<0) {
677 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
678 }
679 }
680 }
681
682 threshold = theta*max_offdiagonal;
683 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
684 j=A->pattern->index[iptr];
685 if((-A->val[iptr])>=threshold) {
686 Paso_IndexList_insertIndex(&(index_list[i]),j);
687 }
688 }
689 }
690 }
691
692
693 S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
694
695 lambda=TMPMEMALLOC(n,dim_t);
696
697 for (i=0;i<n;++i) {
698 lambda[i]=-1;
699 }
700
701 /*S_i={j \in N_i; i strongly coupled to j}*/
702
703 k=0;
704 maxlambda=0;
705
706 for (i=0;i<n;++i) {
707 if(mis_marker[i]==IS_AVAILABLE) {
708 lambda[i]=how_many(i,S,TRUE);
709 /*printf("lambda[%d]=%d, ",i,lambda[i]);*/
710 if(maxlambda<lambda[i]) {
711 maxlambda=lambda[i];
712 index_maxlambda=i;
713 }
714 }
715 }
716
717 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
718
719 index_maxlambda=arg_max(n,lambda, -1);
720 if(index_maxlambda<0) {
721 break;
722 }
723
724 for (i=0;i<n;++i) {
725 if(mis_marker[i]==IS_AVAILABLE) {
726 if (i==index_maxlambda) {
727 mis_marker[index_maxlambda]=IS_IN_C;
728 lambda[index_maxlambda]=-1;
729 for (j=0;j<n;++j) {
730 if(mis_marker[j]==IS_AVAILABLE) {
731 index=&(S->index[S->ptr[j]]);
732 where_p=(index_t*)bsearch(&i,
733 index,
734 S->ptr[j + 1]-S->ptr[j],
735 sizeof(index_t),
736 Paso_comparIndex);
737 if (where_p!=NULL) {
738 mis_marker[j]=IS_IN_F;
739 lambda[j]=-1;
740 for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
741 k=S->index[iptr];
742 if(mis_marker[k]==IS_AVAILABLE) {
743 lambda[k]++;
744 }
745 }
746 }
747
748 }
749 }
750
751 }
752 }
753 }
754
755 }
756
757 for (i=0;i<n;++i)
758 if(mis_marker[i]==IS_AVAILABLE)
759 mis_marker[i]=IS_IN_F;
760
761 /*update lambdas*/
762 /*for (i=0;i<n;++i) {
763 if(mis_marker[i]==IS_AVAILABLE) {
764 lambda[i]=how_many(n,S_T[i], IS_STRONG, mis_marker, IS_AVAILABLE)+2*how_many(n,S_T[i], IS_STRONG, mis_marker, IS_IN_F);
765 if(maxlambda<lambda[i]) {
766 maxlambda=lambda[i];
767 index_maxlambda=i;
768 }
769 }
770 if(lambda[i]==0) {
771 breakloop=TRUE;
772 break;
773 }
774 }
775 if(breakloop) {
776 break;
777 }
778
779 for (i=0;i<n;++i) {
780 if(mis_marker[i]==IS_AVAILABLE) {
781 mis_marker[index_maxlambda]=IS_IN_C;
782 }
783
784 for (j=0;j<n;++j) {
785 if(S_T_[i][j]=IS_STRONG && mis_marker[i]==IS_AVAILABLE) {
786 mis_marker[j]==IS__IN_F;
787 }
788 }
789 }
790
791 }
792 */
793
794 TMPMEMFREE(lambda);
795
796 /* clean up */
797 if (index_list!=NULL) {
798 #pragma omp parallel for private(i) schedule(static)
799 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
800 }
801
802 TMPMEMFREE(index_list);
803 Paso_Pattern_free(S);
804
805
806 /* swap to TRUE/FALSE in mis_marker */
807 #pragma omp parallel for private(i) schedule(static)
808 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
809
810 }
811
812 /*Used in Paso_Pattern_RS_MI*/
813
814 dim_t how_many(dim_t i,Paso_Pattern * S, bool_t transpose) {
815 dim_t j,n;
816 dim_t total;
817 index_t iptr,*index,*where_p;
818 total=0;
819
820 n=S->numOutput;
821
822 if(transpose) {
823 for (j=0;j<n;++j) {
824 index=&(S->index[S->ptr[j]]);
825 where_p=(index_t*)bsearch(&i,
826 index,
827 S->ptr[j + 1]-S->ptr[j],
828 sizeof(index_t),
829 Paso_comparIndex);
830 if (where_p!=NULL) {
831 total++;
832 }
833 }
834 }
835 else {
836 for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
837 total++;
838 }
839
840 }
841 return total;
842 }
843
844 dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {
845 dim_t i;
846 dim_t max=0;
847 dim_t argmax=-1;
848 for (i=0;i<n;++i) {
849 if(max<lambda[i] && lambda[i]!=mask){
850 argmax=i;
851 max=lambda[i];
852 }
853 }
854 return argmax;
855 }
856
857
858 /*dim_t how_many(dim_t n,dim_t* S_i, int value1, dim_t* addedSet, int value2) {
859 dim_t j;
860 dim_t total;
861 total=0;
862 for (j=0;j<n;++j) {
863 if(S_i[j]==value1 && addedSet[j]==value2)
864 total++;
865 }
866 return total;
867 }
868 */
869
870 #undef IS_AVAILABLE
871 #undef IS_IN_F
872 #undef IS_IN_C
873
874 #undef IS_UNDECIDED
875 #undef IS_STRONG
876 #undef IS_WEAK
877
878 #undef IS_IN_FB
879 #undef IS_IN_FB
880
881
882
883
884
885
886
887
888
889

  ViewVC Help
Powered by ViewVC 1.1.26