/[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 2807 - (show annotations)
Mon Dec 7 00:02:55 2009 UTC (9 years, 10 months ago) by artak
File MIME type: text/plain
File size: 29963 byte(s)
Post and pre smoothing parameters are added for AMLI. Little bit reengineering according to Axelson's original AMLI paper.
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,*index,*where_p,*diagptr;
57 bool_t passed=FALSE;
58 dim_t n=A->numRows;
59 diagptr=MEMALLOC(n,index_t);
60
61 if (A->pattern->type & PATTERN_FORMAT_SYM) {
62 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
63 return;
64 }
65
66 #pragma omp parallel for private(i) schedule(static)
67 for (i=0;i<n;++i)
68 if(mis_marker[i]==IS_AVAILABLE)
69 mis_marker[i]=IS_IN_C;
70
71 #pragma omp parallel for private(i,index,where_p) schedule(static)
72 for (i=0;i<n;++i) {
73 diagptr[i]=A->pattern->ptr[i];
74 index=&(A->pattern->index[A->pattern->ptr[i]]);
75 where_p=(index_t*)bsearch(&i,
76 index,
77 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
78 sizeof(index_t),
79 Paso_comparIndex);
80 if (where_p==NULL) {
81 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
82 } else {
83 diagptr[i]+=(index_t)(where_p-index);
84 }
85 }
86
87 /*This loop cannot be parallelized, as order matters here.*/
88 for (i=0;i<n;++i) {
89 if (mis_marker[i]==IS_IN_C) {
90 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
91 j=A->pattern->index[iptr];
92 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
93 mis_marker[j]=IS_IN_F;
94 }
95 }
96 }
97 }
98
99
100
101 /*This loop cannot be parallelized, as order matters here.*/
102 for (i=0;i<n;i++) {
103 if (mis_marker[i]==IS_IN_F) {
104 passed=TRUE;
105 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
106 j=A->pattern->index[iptr];
107 if (mis_marker[j]==IS_IN_C) {
108 if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
109 passed=TRUE;
110 }
111 else {
112 passed=FALSE;
113 break;
114 }
115 }
116 }
117 if (passed) mis_marker[i]=IS_IN_C;
118 }
119 }
120
121 /* swap to TRUE/FALSE in mis_marker */
122 #pragma omp parallel for private(i) schedule(static)
123 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
124
125 MEMFREE(diagptr);
126 }
127
128
129 /*
130 * Ruge-Stueben strength of connection mask.
131 *
132 */
133 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
134 {
135 dim_t i,n,j;
136 index_t iptr;
137 double threshold,max_offdiagonal;
138
139 Paso_Pattern *out=NULL;
140
141 Paso_IndexList* index_list=NULL;
142
143 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
144 if (! Paso_checkPtr(index_list)) {
145 #pragma omp parallel for private(i) schedule(static)
146 for(i=0;i<A->pattern->numOutput;++i) {
147 index_list[i].extension=NULL;
148 index_list[i].n=0;
149 }
150 }
151
152
153 n=A->numRows;
154 if (A->pattern->type & PATTERN_FORMAT_SYM) {
155 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
156 return;
157 }
158 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
159 for (i=0;i<n;++i) {
160 if(mis_marker[i]==IS_AVAILABLE) {
161 max_offdiagonal = DBL_MIN;
162 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
163 if(A->pattern->index[iptr] != i){
164 max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
165 }
166 }
167
168 threshold = theta*max_offdiagonal;
169 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170 j=A->pattern->index[iptr];
171 if((-A->val[iptr])>=threshold) {
172 Paso_IndexList_insertIndex(&(index_list[i]),j);
173 Paso_IndexList_insertIndex(&(index_list[j]),i);
174 }
175 }
176 }
177 }
178
179
180 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
181
182 /* clean up */
183 if (index_list!=NULL) {
184 #pragma omp parallel for private(i) schedule(static)
185 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
186 }
187 TMPMEMFREE(index_list);
188
189 /*Paso_Pattern_mis(out,mis_marker);*/
190 Paso_Pattern_greedy(out,mis_marker);
191 Paso_Pattern_free(out);
192 }
193
194 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
195 {
196 dim_t i,j,n;
197 index_t iptr;
198 double diag,eps_Aii,val;
199 double* diags;
200
201
202 Paso_Pattern *out=NULL;
203 Paso_IndexList* index_list=NULL;
204
205 n=A->numRows;
206 diags=MEMALLOC(n,double);
207
208 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
209 if (! Paso_checkPtr(index_list)) {
210 #pragma omp parallel for private(i) schedule(static)
211 for(i=0;i<A->pattern->numOutput;++i) {
212 index_list[i].extension=NULL;
213 index_list[i].n=0;
214 }
215 }
216
217 if (A->pattern->type & PATTERN_FORMAT_SYM) {
218 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
219 return;
220 }
221
222
223 #pragma omp parallel for private(i,iptr,diag) schedule(static)
224 for (i=0;i<n;++i) {
225 diag = 0;
226 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
227 if(A->pattern->index[iptr] == i){
228 diag+=A->val[iptr];
229 }
230 }
231 diags[i]=ABS(diag);
232 }
233
234
235 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
236 for (i=0;i<n;++i) {
237 if (mis_marker[i]==IS_AVAILABLE) {
238 eps_Aii = theta*theta*diags[i];
239 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
240 j=A->pattern->index[iptr];
241 val=A->val[iptr];
242 if((val*val)>=(eps_Aii*diags[j])) {
243 Paso_IndexList_insertIndex(&(index_list[i]),j);
244 }
245 }
246 }
247 }
248
249 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
250
251 /* clean up */
252 if (index_list!=NULL) {
253 #pragma omp parallel for private(i) schedule(static)
254 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
255 }
256
257 TMPMEMFREE(index_list);
258 MEMFREE(diags);
259
260
261 /*Paso_Pattern_mis(out,mis_marker);*/
262 Paso_Pattern_greedy(out,mis_marker);
263 Paso_Pattern_free(out);
264
265 }
266
267 /* Greedy algorithm */
268 void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
269
270 dim_t i,j;
271 /*double sum;*/
272 index_t iptr;
273 bool_t passed=FALSE;
274 dim_t n=pattern->numOutput;
275
276 if (pattern->type & PATTERN_FORMAT_SYM) {
277 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
278 return;
279 }
280
281 #pragma omp parallel for private(i) schedule(static)
282 for (i=0;i<n;++i)
283 if(mis_marker[i]==IS_AVAILABLE)
284 mis_marker[i]=IS_IN_C;
285
286
287 for (i=0;i<n;++i) {
288 if (mis_marker[i]==IS_IN_C) {
289 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
290 j=pattern->index[iptr];
291 mis_marker[j]=IS_IN_F;
292 }
293 }
294 }
295
296
297
298 for (i=0;i<n;i++) {
299 if (mis_marker[i]==IS_IN_F) {
300 passed=TRUE;
301 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
302 j=pattern->index[iptr];
303 if (mis_marker[j]==IS_IN_F) {
304 passed=TRUE;
305 }
306 else {
307 passed=FALSE;
308 break;
309 }
310 }
311 if (passed) mis_marker[i]=IS_IN_C;
312 }
313 }
314
315 /* swap to TRUE/FALSE in mis_marker */
316 #pragma omp parallel for private(i) schedule(static)
317 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
318
319 }
320
321
322 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
323
324 dim_t i,j;
325 /*double sum;*/
326 index_t iptr;
327 index_t num_colors;
328 index_t* colorOf;
329 register index_t color;
330 bool_t passed=FALSE;
331 dim_t n=pattern->numOutput;
332
333
334 colorOf=MEMALLOC(n,index_t);
335
336 if (pattern->type & PATTERN_FORMAT_SYM) {
337 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
338 return;
339 }
340
341 Paso_Pattern_color(pattern,&num_colors,colorOf);
342
343 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
344 #pragma omp parallel for private(i) schedule(static)
345 for (i=0;i<n;++i)
346 if(mis_marker[i]==IS_AVAILABLE)
347 mis_marker[i]=IS_IN_F;
348
349 #pragma omp barrier
350 for (color=0;color<num_colors;++color) {
351 #pragma omp parallel for schedule(static) private(i,iptr,j)
352 for (i=0;i<n;++i) {
353 if (colorOf[i]==color) {
354 if (mis_marker[i]==IS_IN_F) {
355 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
356 j=pattern->index[iptr];
357 if (colorOf[j]<color)
358 mis_marker[j]=IS_IN_C;
359 }
360 }
361 }
362 }
363 }
364
365
366 #pragma omp barrier
367 for (color=0;color<num_colors;++color) {
368 #pragma omp parallel for schedule(static) private(i,iptr,j)
369 for (i=0;i<n;i++) {
370 if (colorOf[i]==color) {
371 if (mis_marker[i]==IS_IN_C) {
372 passed=TRUE;
373 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
374 j=pattern->index[iptr];
375 if (colorOf[j]<color && passed) {
376 if (mis_marker[j]==IS_IN_C) {
377 passed=TRUE;
378 }
379 else {
380 passed=FALSE;
381 /*break;*/
382 }
383 }
384 }
385 if (passed) mis_marker[i]=IS_IN_F;
386 }
387
388 }
389 }
390 }
391
392 /* swap to TRUE/FALSE in mis_marker */
393 #pragma omp parallel for private(i) schedule(static)
394 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
395
396 MEMFREE(colorOf);
397 }
398
399 /*For testing */
400 void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
401
402 dim_t i,j=0,k;
403 double *theta;
404 index_t iptr;
405 dim_t n=A->numRows;
406 double rsum,diag=0;
407 index_t *AvADJ;
408 theta=MEMALLOC(n,double);
409 AvADJ=MEMALLOC(n,index_t);
410
411
412
413
414 if (A->pattern->type & PATTERN_FORMAT_SYM) {
415 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
416 return;
417 }
418
419
420 #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
421 for (i=0;i<n;++i) {
422 rsum=0;
423 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
424 j=A->pattern->index[iptr];
425 if(j!=i) {
426 rsum+=ABS(A->val[iptr]);
427 }
428 else {
429 diag=ABS(A->val[iptr]);
430 }
431 }
432 theta[i]=diag/rsum;
433 if(theta[i]>threshold) {
434 mis_marker[i]=IS_IN_F;
435 }
436 }
437
438 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
439 k=0;
440
441 for (i=0;i<n;++i) {
442 if(mis_marker[i]==IS_AVAILABLE) {
443 if(k==0) {
444 j=i;
445 k++;
446 }
447 if(theta[j]>theta[i]) {
448 j=i;
449 }
450 }
451 }
452 mis_marker[j]=IS_IN_C;
453
454 for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
455 k=A->pattern->index[iptr];
456 if(mis_marker[k]==IS_AVAILABLE) {
457 AvADJ[k]=1;
458 }
459 else {
460 AvADJ[k]=-1;
461 }
462
463 }
464
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,jptr;
644 /*index_t *index,*where_p;*/
645 double threshold,max_offdiagonal;
646 dim_t *lambda; /*mesure of importance */
647 /*bool_t breakloop=FALSE;*/
648 dim_t maxlambda=0;
649 index_t index_maxlambda=0;
650 double time0=0;
651 bool_t verbose=1;
652
653 Paso_Pattern *S=NULL;
654 Paso_Pattern *ST=NULL;
655 Paso_IndexList* index_list=NULL;
656
657 index_t* counter;
658 /*dim_t lk;*/
659
660 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
661 if (! Paso_checkPtr(index_list)) {
662 #pragma omp parallel for private(i) schedule(static)
663 for(i=0;i<A->pattern->numOutput;++i) {
664 index_list[i].extension=NULL;
665 index_list[i].n=0;
666 }
667 }
668
669
670 n=A->numRows;
671 if (A->pattern->type & PATTERN_FORMAT_SYM) {
672 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
673 return;
674 }
675
676 time0=Paso_timer();
677 k=0;
678 /*S_i={j \in N_i; i strongly coupled to j}*/
679 #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)
680 for (i=0;i<n;++i) {
681 if(mis_marker[i]==IS_AVAILABLE) {
682 max_offdiagonal = DBL_MIN;
683 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
684 if(A->pattern->index[iptr] != i){
685 if(A->val[iptr]<0) {
686 max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
687 }
688 }
689 }
690 threshold = theta*max_offdiagonal;
691 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
692 j=A->pattern->index[iptr];
693 if((-A->val[iptr])>=threshold) {
694 Paso_IndexList_insertIndex(&(index_list[i]),j);
695 }
696 }
697 }
698 }
699
700
701 S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
702 ST=Paso_Pattern_getTranspose(S);
703
704 time0=Paso_timer()-time0;
705 if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
706
707 lambda=TMPMEMALLOC(n,dim_t);
708
709 #pragma omp parallel for private(i) schedule(static)
710 for (i=0;i<n;++i) { lambda[i]=-1; }
711
712 /*S_i={j \in N_i; i strongly coupled to j}*/
713
714 /*
715 #pragma omp parallel for private(i,iptr,lk) schedule(static)
716 for (i=0;i<n;++i) {
717 lk=0;
718 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
719 if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )
720 lk++;
721 }
722 #pragma omp critical
723 k+=lk;
724 if(k==0) {
725 mis_marker[i]=IS_IN_F;
726 }
727 }
728 */
729
730
731 k=0;
732 maxlambda=0;
733
734 time0=Paso_timer();
735
736 for (i=0;i<n;++i) {
737 if(mis_marker[i]==IS_AVAILABLE) {
738 lambda[i]=how_many(k,ST,FALSE);
739 /*lambda[i]=how_many(i,S,TRUE);*/
740 /*printf("lambda[%d]=%d, k=%d ",i,lambda[i],k);*/
741 k++;
742 if(maxlambda<lambda[i]) {
743 maxlambda=lambda[i];
744 index_maxlambda=i;
745 }
746 }
747 }
748
749 k=0;
750 time0=Paso_timer()-time0;
751 if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);
752
753 time0=Paso_timer();
754
755 /*Paso_Pattern_getReport(n,mis_marker);*/
756 counter=TMPMEMALLOC(n,index_t);
757
758 #pragma omp parallel for private(i) schedule(static)
759 for (i = 0; i < n; ++i) counter[i]=(mis_marker[i]==IS_AVAILABLE);
760
761 Paso_Util_cumsum(n,counter);
762
763
764 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
765
766 if(index_maxlambda<0) {
767 break;
768 }
769
770 i=index_maxlambda;
771 if(mis_marker[i]==IS_AVAILABLE) {
772 mis_marker[index_maxlambda]=IS_IN_C;
773 lambda[index_maxlambda]=-1;
774 for (iptr=ST->ptr[counter[i]];iptr<ST->ptr[counter[i]+1]; ++iptr) {
775 j=ST->index[iptr];
776 if(mis_marker[j]==IS_AVAILABLE) {
777 mis_marker[j]=IS_IN_F;
778 lambda[j]=-1;
779 for (jptr=S->ptr[counter[j]];jptr<S->ptr[counter[j]+1]; ++jptr) {
780 k=S->index[jptr];
781 if(mis_marker[k]==IS_AVAILABLE) {
782 lambda[k]++;
783 }
784 }
785 }
786 }
787 }
788
789 /* Used when transpose of S is not available */
790 /*
791 for (i=0;i<n;++i) {
792 if(mis_marker[i]==IS_AVAILABLE) {
793 if (i==index_maxlambda) {
794 mis_marker[index_maxlambda]=IS_IN_C;
795 lambda[index_maxlambda]=-1;
796 for (j=0;j<n;++j) {
797 if(mis_marker[j]==IS_AVAILABLE) {
798 index=&(S->index[S->ptr[j]]);
799 where_p=(index_t*)bsearch(&i,
800 index,
801 S->ptr[j + 1]-S->ptr[j],
802 sizeof(index_t),
803 Paso_comparIndex);
804 if (where_p!=NULL) {
805 mis_marker[j]=IS_IN_F;
806 lambda[j]=-1;
807 for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
808 k=S->index[iptr];
809 if(mis_marker[k]==IS_AVAILABLE) {
810 lambda[k]++;
811 }
812 }
813 }
814
815 }
816 }
817
818 }
819 }
820 }
821 */
822 index_maxlambda=arg_max(n,lambda, -1);
823 }
824
825 time0=Paso_timer()-time0;
826 if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
827
828 /*Paso_Pattern_getReport(n,mis_marker);*/
829
830 #pragma omp parallel for private(i) schedule(static)
831 for (i=0;i<n;++i)
832 if(mis_marker[i]==IS_AVAILABLE || mis_marker[i]==IS_NOT_AVAILABLE) {
833 mis_marker[i]=IS_IN_F;
834 }
835
836 /*Paso_Pattern_getReport(n,mis_marker);*/
837
838 TMPMEMFREE(lambda);
839
840 /* clean up */
841 if (index_list!=NULL) {
842 #pragma omp parallel for private(i) schedule(static)
843 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
844 }
845
846 TMPMEMFREE(index_list);
847 Paso_Pattern_free(S);
848
849 /* swap to TRUE/FALSE in mis_marker */
850 #pragma omp parallel for private(i) schedule(static)
851 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
852
853 }
854
855 void Paso_Pattern_RS_MI_Aggressive(Paso_SparseMatrix* A, index_t* mis_marker, double theta) {
856 dim_t i;
857
858 Paso_Pattern_RS_MI(A, mis_marker, theta);
859
860 #pragma omp parallel for private(i) schedule(static)
861 for (i=0;i<A->numRows;i++) {
862 if(!mis_marker[i]) {
863 mis_marker[i]=IS_AVAILABLE;
864 }
865 else {
866 mis_marker[i]=IS_NOT_AVAILABLE;
867 }
868 }
869
870 Paso_Pattern_RS_MI(A, mis_marker, theta);
871 }
872
873 void Paso_Pattern_getReport(dim_t n,index_t* mis_marker) {
874
875 dim_t i,av,nav,f,c;
876
877 av=0;
878 nav=0;
879 f=0;
880 c=0;
881
882 for (i=0;i<n;i++) {
883 if(mis_marker[i]==IS_NOT_AVAILABLE) {
884 nav++;
885 }
886 else if (mis_marker[i]==IS_AVAILABLE) {
887 av++;
888 }
889 else if (mis_marker[i]==IS_IN_F) {
890 f++;
891 }
892 else if (mis_marker[i]==IS_IN_C) {
893 c++;
894 }
895 }
896
897 printf("Available=%d, Not Available = %d, In F = %d, In C = %d \n",av,nav,f,c);
898
899 }
900
901
902 /*Used in Paso_Pattern_RS_MI*/
903
904 /*Computes how many connections unknown i has in S.
905 bool_t transpose - TRUE if we want to compute how many strong connection of i in S^T, FALSE otherwise.
906 Note that if we first transpose S and then call method on S^T then "transpose" should be set to FALSE.
907 */
908 dim_t how_many(dim_t i,Paso_Pattern * S, bool_t transpose) {
909 dim_t j,n;
910 dim_t total,ltotal;
911 index_t *index,*where_p;
912
913 /*index_t iptr;*/
914 total=0;
915 ltotal=0;
916
917 n=S->numOutput;
918
919 if(transpose) {
920 #pragma omp parallel
921 {
922 ltotal=0;
923 #pragma omp for private(j,index,where_p,ltotal) schedule(static)
924 for (j=0;j<n;++j) {
925 index=&(S->index[S->ptr[j]]);
926 where_p=(index_t*)bsearch(&i,
927 index,
928 S->ptr[j + 1]-S->ptr[j],
929 sizeof(index_t),
930 Paso_comparIndex);
931 if (where_p!=NULL) {
932 ltotal++;
933 }
934 }
935 }
936 #pragma omp critical
937 {
938 total+=ltotal;
939 }
940
941 }
942 else {
943 total=S->ptr[i+1]-S->ptr[i];
944 /*#pragma omp parallel for private(iptr) schedule(static)*/
945 /*for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
946 if(S->index[iptr]!=i && mis_marker[S->index[iptr]]==IS_AVAILABLE)
947 total++;
948 }*/
949
950 }
951 return total;
952 }
953
954 dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {
955 dim_t i;
956 dim_t max=0;
957 dim_t argmax=-1;
958
959 #ifdef _OPENMP
960 dim_t lmax=0;
961 dim_t li=-1;
962 #pragma omp parallel private(i,lmax,li)
963 {
964 lmax=0;
965 li=-1;
966 #pragma omp for schedule(static)
967 for (i=0;i<n;++i) {
968 if(lmax<lambda[i] && lambda[i]!=mask){
969 lmax=lambda[i];
970 li=i;
971 }
972 }
973 #pragma omp critical
974 {
975 if (max<=lmax) {
976 if(max==lmax && argmax>li)
977 {
978 argmax=li;
979 }
980 if (max < lmax)
981 {
982 max=lmax;
983 argmax=li;
984 }
985 }
986 }
987 }
988 #else
989 for (i=0;i<n;++i) {
990 if(max<lambda[i] && lambda[i]!=mask){
991 max=lambda[i];
992 argmax=i;
993 }
994 }
995 #endif
996
997 return argmax;
998 }
999
1000
1001 Paso_Pattern* Paso_Pattern_getTranspose(Paso_Pattern* P){
1002
1003 Paso_Pattern *outpattern=NULL;
1004
1005 Paso_IndexList* index_list=NULL;
1006 dim_t C=P->numInput;
1007 dim_t F=P->numOutput-C;
1008 dim_t n=C+F;
1009 dim_t i,j;
1010 index_t iptr;
1011
1012 index_list=TMPMEMALLOC(C,Paso_IndexList);
1013 if (! Paso_checkPtr(index_list)) {
1014 #pragma omp parallel for private(i) schedule(static)
1015 for(i=0;i<C;++i) {
1016 index_list[i].extension=NULL;
1017 index_list[i].n=0;
1018 }
1019 }
1020
1021
1022 /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
1023 for (i=0;i<n;++i) {
1024 for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {
1025 j=P->index[iptr];
1026 Paso_IndexList_insertIndex(&(index_list[j]),i);
1027 }
1028 }
1029
1030 outpattern=Paso_IndexList_createPattern(0, C,index_list,0,n,0);
1031
1032 /* clean up */
1033 if (index_list!=NULL) {
1034 #pragma omp parallel for private(i) schedule(static)
1035 for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
1036 }
1037 TMPMEMFREE(index_list);
1038
1039 return outpattern;
1040 }
1041
1042
1043 #undef IS_AVAILABLE
1044 #undef IS_NOT_AVAILABLE
1045 #undef IS_IN_F
1046 #undef IS_IN_C
1047
1048 #undef IS_UNDECIDED
1049 #undef IS_STRONG
1050 #undef IS_WEAK
1051
1052 #undef IS_IN_FB
1053 #undef IS_IN_FB
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063

  ViewVC Help
Powered by ViewVC 1.1.26