/[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 3010 - (show annotations)
Tue Apr 27 05:10:46 2010 UTC (9 years, 5 months ago) by artak
File MIME type: text/plain
File size: 41191 byte(s)
Preparation for AMG on Systems without cut using frobenius norm
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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
53 void Paso_Pattern_Read(char *fileName,dim_t n,index_t* mis_marker) {
54
55 dim_t i;
56 int scan_ret;
57
58 FILE *fileHandle_p = NULL;
59
60 fileHandle_p = fopen( fileName, "r" );
61 if( fileHandle_p == NULL )
62 {
63 Paso_setError(IO_ERROR, "Paso_Pattern_Read: Cannot open file for reading.");
64 return;
65 }
66
67 for (i=0;i<n;++i) {
68 scan_ret=fscanf( fileHandle_p, "%d\n", &mis_marker[i]);
69 if (scan_ret!=1) {
70 Paso_setError(IO_ERROR, "Paso_Pattern_Read: Cannot read line from file.");
71 return;
72 }
73 }
74
75 fclose(fileHandle_p);
76
77 /* swap to TRUE/FALSE in mis_marker */
78 #pragma omp parallel for private(i) schedule(static)
79 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=1);
80 }
81
82 void Paso_Pattern_Write(char *fileName,dim_t n,index_t* mis_marker) {
83
84 dim_t i;
85
86 FILE *fileHandle_p = NULL;
87
88 fileHandle_p = fopen( fileName, "w+" );
89 if( fileHandle_p == NULL )
90 {
91 Paso_setError(IO_ERROR, "Paso_Pattern_Write: Cannot open file for writing.");
92 return;
93 }
94
95 for (i=0;i<n;++i) {
96 fprintf(fileHandle_p, "%d\n", mis_marker[i]);
97 }
98
99 fclose(fileHandle_p);
100
101 }
102
103
104 void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
105
106 dim_t i,j,bi;
107 /*double sum;*/
108 index_t iptr;
109 /*index_t *index,*where_p;*/
110 bool_t passed=FALSE;
111 dim_t n=A->numRows;
112 dim_t block_size=A->row_block_size;
113 double *diags;
114 double fnorm;
115 diags=MEMALLOC(n,double);
116
117
118 if (A->pattern->type & PATTERN_FORMAT_SYM) {
119 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
120 return;
121 }
122
123 #pragma omp parallel for private(i) schedule(static)
124 for (i=0;i<n;++i)
125 if(mis_marker[i]==IS_AVAILABLE)
126 mis_marker[i]=IS_IN_C;
127
128 #pragma omp parallel for private(i,j,fnorm,bi) schedule(static)
129 for (i=0;i<n;++i) {
130 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
131 j=A->pattern->index[iptr];
132 if(i==j) {
133 if(block_size==1) {
134 diags[i]=A->val[iptr];
135 } else {
136 fnorm=0;
137 for(bi=0;bi<block_size*block_size;++bi)
138 {
139 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
140 }
141 fnorm=sqrt(fnorm);
142 diags[i]=fnorm;
143 }
144
145 }
146 }
147 }
148
149 /*This loop cannot be parallelized, as order matters here.*/
150 for (i=0;i<n;++i) {
151 if (mis_marker[i]==IS_IN_C) {
152 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
153 j=A->pattern->index[iptr];
154 if(block_size==1) {
155 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(diags[i])) {
156 mis_marker[j]=IS_IN_F;
157 }
158 } else {
159 if (j!=i) {
160 fnorm=0;
161 for(bi=0;bi<block_size*block_size;++bi)
162 {
163 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
164 }
165 fnorm=sqrt(fnorm);
166 if (fnorm>=threshold*diags[i]) {
167 mis_marker[j]=IS_IN_F;
168 }
169 }
170 }
171 }
172 }
173 }
174
175 /*This loop cannot be parallelized, as order matters here.*/
176 for (i=0;i<n;i++) {
177 if (mis_marker[i]==IS_IN_F) {
178 passed=TRUE;
179 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
180 j=A->pattern->index[iptr];
181 if (mis_marker[j]==IS_IN_C) {
182 if(block_size==1) {
183 if ((A->val[iptr]/diags[i])>=-threshold) {
184 passed=TRUE;
185 }
186 else {
187 passed=FALSE;
188 break;
189 }
190 } else {
191 fnorm=0;
192 for(bi=0;bi<block_size*block_size;++bi)
193 {
194 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
195 }
196 fnorm=sqrt(fnorm);
197 if ((fnorm/diags[i])<=threshold) {
198 passed=TRUE;
199 }
200 else {
201 passed=FALSE;
202 break;
203 }
204 }
205
206 }
207
208 }
209 if (passed) mis_marker[i]=IS_IN_C;
210 }
211 }
212
213 /* swap to TRUE/FALSE in mis_marker */
214 #pragma omp parallel for private(i) schedule(static)
215 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
216
217 MEMFREE(diags);
218 }
219
220
221 /*
222 * Ruge-Stueben strength of connection mask.
223 *
224 */
225 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
226 {
227 dim_t i,n,j,bi;
228 index_t iptr;
229 double threshold,max_offdiagonal;
230 double fnorm;
231 dim_t block_size=A->row_block_size;
232
233 Paso_Pattern *out=NULL;
234
235 Paso_IndexList* index_list=NULL;
236
237 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
238 if (! Paso_checkPtr(index_list)) {
239 #pragma omp parallel for private(i) schedule(static)
240 for(i=0;i<A->pattern->numOutput;++i) {
241 index_list[i].extension=NULL;
242 index_list[i].n=0;
243 }
244 }
245
246
247 n=A->numRows;
248 if (A->pattern->type & PATTERN_FORMAT_SYM) {
249 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
250 return;
251 }
252 /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j,fnorm,bi) schedule(static)*/
253 for (i=0;i<n;++i) {
254 if(mis_marker[i]==IS_AVAILABLE) {
255 max_offdiagonal = DBL_MIN;
256 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
257 if(A->pattern->index[iptr] != i){
258 if(block_size==1) {
259 max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
260 } else {
261 fnorm=0;
262 for(bi=0;bi<block_size*block_size;++bi)
263 {
264 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
265 }
266 fnorm=sqrt(fnorm);
267 max_offdiagonal = MAX(max_offdiagonal,fnorm);
268 }
269 }
270 }
271
272 threshold = theta*max_offdiagonal;
273 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
274 j=A->pattern->index[iptr];
275 if(block_size==1) {
276 if(ABS(A->val[iptr])>=threshold && i!=j) {
277 Paso_IndexList_insertIndex(&(index_list[i]),j);
278 /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
279 }
280 } else {
281 if (i!=j) {
282 fnorm=0;
283 for(bi=0;bi<block_size*block_size;++bi)
284 {
285 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
286 }
287 fnorm=sqrt(fnorm);
288 if(fnorm>=threshold) {
289 Paso_IndexList_insertIndex(&(index_list[i]),j);
290 /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
291 }
292 }
293 }
294 }
295 }
296 }
297
298
299 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
300
301 /* clean up */
302 if (index_list!=NULL) {
303 #pragma omp parallel for private(i) schedule(static)
304 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
305 }
306 TMPMEMFREE(index_list);
307
308 /*Paso_Pattern_mis(out,mis_marker);*/
309 Paso_Pattern_greedy(out,mis_marker);
310 Paso_Pattern_free(out);
311 }
312
313 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
314 {
315 dim_t i,j,n,bi;
316 index_t iptr;
317 double diag,eps_Aii,val;
318 double* diags;
319 double fnorm;
320 dim_t block_size=A->row_block_size;
321
322
323 Paso_Pattern *out=NULL;
324 Paso_IndexList* index_list=NULL;
325
326 n=A->numRows;
327 diags=MEMALLOC(n,double);
328
329 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
330 if (! Paso_checkPtr(index_list)) {
331 #pragma omp parallel for private(i) schedule(static)
332 for(i=0;i<A->pattern->numOutput;++i) {
333 index_list[i].extension=NULL;
334 index_list[i].n=0;
335 }
336 }
337
338 if (A->pattern->type & PATTERN_FORMAT_SYM) {
339 Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
340 return;
341 }
342
343
344 #pragma omp parallel for private(i,iptr,diag,fnorm,bi) schedule(static)
345 for (i=0;i<n;++i) {
346 diag = 0;
347 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
348 if(A->pattern->index[iptr] == i){
349 if(block_size==1) {
350 diag+=A->val[iptr];
351 } else {
352 fnorm=0;
353 for(bi=0;bi<block_size*block_size;++bi)
354 {
355 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
356 }
357 fnorm=sqrt(fnorm);
358 diag+=fnorm;
359 }
360
361 }
362 }
363 diags[i]=ABS(diag);
364 }
365
366
367 #pragma omp parallel for private(i,iptr,j,val,eps_Aii,fnorm,bi) schedule(static)
368 for (i=0;i<n;++i) {
369 if (mis_marker[i]==IS_AVAILABLE) {
370 eps_Aii = theta*theta*diags[i];
371 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
372 j=A->pattern->index[iptr];
373 if(block_size==1) {
374 val=A->val[iptr];
375 } else {
376 fnorm=0;
377 for(bi=0;bi<block_size*block_size;++bi)
378 {
379 fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
380 }
381 fnorm=sqrt(fnorm);
382 val=fnorm;
383 }
384
385 if((val*val)>=(eps_Aii*diags[j])) {
386 Paso_IndexList_insertIndex(&(index_list[i]),j);
387 }
388 }
389 }
390 }
391
392 out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
393
394 /* clean up */
395 if (index_list!=NULL) {
396 #pragma omp parallel for private(i) schedule(static)
397 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
398 }
399
400 TMPMEMFREE(index_list);
401 MEMFREE(diags);
402
403
404 /*Paso_Pattern_mis(out,mis_marker);*/
405 Paso_Pattern_greedy(out,mis_marker);
406 Paso_Pattern_free(out);
407
408 }
409
410 /* Greedy algorithm */
411 void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
412
413 dim_t i,j;
414 /*double sum;*/
415 index_t iptr;
416 bool_t passed=FALSE;
417 dim_t n=pattern->numOutput;
418
419 if (pattern->type & PATTERN_FORMAT_SYM) {
420 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
421 return;
422 }
423
424 #pragma omp parallel for private(i) schedule(static)
425 for (i=0;i<n;++i)
426 if(mis_marker[i]==IS_AVAILABLE)
427 mis_marker[i]=IS_IN_C;
428
429
430 for (i=0;i<n;++i) {
431 if (mis_marker[i]==IS_IN_C) {
432 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
433 j=pattern->index[iptr];
434 mis_marker[j]=IS_IN_F;
435 }
436 }
437 }
438
439
440 for (i=0;i<n;i++) {
441 if (mis_marker[i]==IS_IN_F) {
442 passed=TRUE;
443 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
444 j=pattern->index[iptr];
445 if (mis_marker[j]==IS_IN_F) {
446 passed=TRUE;
447 }
448 else {
449 passed=FALSE;
450 break;
451 }
452 }
453 if (passed) mis_marker[i]=IS_IN_C;
454 }
455 }
456
457 /* swap to TRUE/FALSE in mis_marker */
458 #pragma omp parallel for private(i) schedule(static)
459 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
460
461 }
462
463 void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
464
465 dim_t i,j;
466 /*double sum;*/
467 index_t iptr;
468 index_t num_colors;
469 index_t* colorOf;
470 register index_t color;
471 bool_t passed=FALSE;
472 dim_t n=pattern->numOutput;
473
474
475 colorOf=MEMALLOC(n,index_t);
476
477 if (pattern->type & PATTERN_FORMAT_SYM) {
478 Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
479 return;
480 }
481
482 Paso_Pattern_color(pattern,&num_colors,colorOf);
483
484 /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
485 #pragma omp parallel for private(i) schedule(static)
486 for (i=0;i<n;++i)
487 if(mis_marker[i]==IS_AVAILABLE)
488 mis_marker[i]=IS_IN_F;
489
490 #pragma omp barrier
491 for (color=0;color<num_colors;++color) {
492 #pragma omp parallel for schedule(static) private(i,iptr,j)
493 for (i=0;i<n;++i) {
494 if (colorOf[i]==color) {
495 if (mis_marker[i]==IS_IN_F) {
496 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
497 j=pattern->index[iptr];
498 if (colorOf[j]<color)
499 mis_marker[j]=IS_IN_C;
500 }
501 }
502 }
503 }
504 }
505
506
507 #pragma omp barrier
508 for (color=0;color<num_colors;++color) {
509 #pragma omp parallel for schedule(static) private(i,iptr,j)
510 for (i=0;i<n;i++) {
511 if (colorOf[i]==color) {
512 if (mis_marker[i]==IS_IN_C) {
513 passed=TRUE;
514 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
515 j=pattern->index[iptr];
516 if (colorOf[j]<color && passed) {
517 if (mis_marker[j]==IS_IN_C) {
518 passed=TRUE;
519 }
520 else {
521 passed=FALSE;
522 /*break;*/
523 }
524 }
525 }
526 if (passed) mis_marker[i]=IS_IN_F;
527 }
528
529 }
530 }
531 }
532
533 /* swap to TRUE/FALSE in mis_marker */
534 #pragma omp parallel for private(i) schedule(static)
535 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
536
537 MEMFREE(colorOf);
538 }
539
540 /*For testing */
541 void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
542
543 dim_t i,j=0,k;
544 double *theta;
545 index_t iptr;
546 dim_t n=A->numRows;
547 double rsum,diag=0;
548 index_t *AvADJ;
549 theta=MEMALLOC(n,double);
550 AvADJ=MEMALLOC(n,index_t);
551
552
553
554
555 if (A->pattern->type & PATTERN_FORMAT_SYM) {
556 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
557 return;
558 }
559
560
561 #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
562 for (i=0;i<n;++i) {
563 rsum=0;
564 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
565 j=A->pattern->index[iptr];
566 if(j!=i) {
567 rsum+=ABS(A->val[iptr]);
568 }
569 else {
570 diag=ABS(A->val[iptr]);
571 }
572 }
573 theta[i]=diag/rsum;
574 if(theta[i]>threshold) {
575 mis_marker[i]=IS_IN_F;
576 }
577 }
578
579 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
580 k=0;
581
582 for (i=0;i<n;++i) {
583 if(mis_marker[i]==IS_AVAILABLE) {
584 if(k==0) {
585 j=i;
586 k++;
587 }
588 if(theta[j]>theta[i]) {
589 j=i;
590 }
591 }
592 }
593 mis_marker[j]=IS_IN_C;
594
595 for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
596 k=A->pattern->index[iptr];
597 if(mis_marker[k]==IS_AVAILABLE) {
598 AvADJ[k]=1;
599 }
600 else {
601 AvADJ[k]=-1;
602 }
603
604 }
605
606 for (i=0;i<n;++i) {
607 if(AvADJ[i]) {
608 rsum=0;
609 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
610 k=A->pattern->index[iptr];
611 if(k!=i && mis_marker[k]!=IS_IN_C ) {
612 rsum+=ABS(A->val[iptr]);
613 }
614 if(j==i) {
615 diag=ABS(A->val[iptr]);
616 }
617 }
618 theta[i]=diag/rsum;
619 if(theta[i]>threshold) {
620 mis_marker[i]=IS_IN_F;
621 }
622 }
623 }
624
625
626 }
627
628 /* swap to TRUE/FALSE in mis_marker */
629 #pragma omp parallel for private(i) schedule(static)
630 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
631
632 MEMFREE(AvADJ);
633 MEMFREE(theta);
634 }
635
636
637 void Paso_Pattern_YS_plus(Paso_SparseMatrix* A, index_t* mis_marker, double alpha, double taw, double delta) {
638
639 dim_t i,j;
640 /*double sum;*/
641 index_t iptr,*index,*where_p,*diagptr;
642 double *rsum;
643 double sum;
644 dim_t n=A->numRows;
645 Paso_SparseMatrix* A_alpha;
646 diagptr=MEMALLOC(n,index_t);
647 rsum=MEMALLOC(n,double);
648
649 if (A->pattern->type & PATTERN_FORMAT_SYM) {
650 Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
651 return;
652 }
653
654 A_alpha=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern,1,1, FALSE);
655 #pragma omp parallel for private(i) schedule(static)
656 for (i=0;i<A->len;++i) {
657 A_alpha->val[i]=A->val[i];
658 }
659
660
661 #pragma omp parallel for private(i) schedule(static)
662 for (i=0;i<n;++i)
663 if(mis_marker[i]==IS_AVAILABLE)
664 mis_marker[i]=IS_IN_C;
665
666 /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
667 for (i=0;i<n;++i) {
668 diagptr[i]=A->pattern->ptr[i];
669 index=&(A->pattern->index[A->pattern->ptr[i]]);
670 where_p=(index_t*)bsearch(&i,
671 index,
672 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
673 sizeof(index_t),
674 Paso_comparIndex);
675 if (where_p==NULL) {
676 Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
677 } else {
678 diagptr[i]+=(index_t)(where_p-index);
679 }
680 }
681
682
683 j=0;
684 for (i=0;i<n;++i) {
685 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
686 j=A_alpha->pattern->index[iptr];
687 if(i==j) {
688 A_alpha->val[iptr]=0;
689 }
690 else {
691 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) ) {
692 A_alpha->val[iptr]=0;
693 }
694 }
695
696 }
697 }
698
699
700 for (i=0;i<n;++i) {
701 rsum[i]=0;
702 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
703 rsum[i]+=A_alpha->val[iptr];
704 }
705 }
706
707 #pragma omp parallel for private(i,j) schedule(static)
708 for (i=0;i<n;++i) {
709 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
710 j=A_alpha->pattern->index[iptr];
711 if (i==j) {
712 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr]+rsum[i];
713 }
714 else {
715 A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr];
716 }
717
718 }
719 }
720
721 /*This loop cannot be parallelized, as order matters here.*/
722 for (i=0;i<n;++i) {
723 if (mis_marker[i]==IS_IN_C) {
724 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
725 j=A->pattern->index[iptr];
726 if (j!=i && ABS(A->val[iptr])>=alpha*MIN(ABS(A->val[diagptr[i]]),ABS(A->val[diagptr[j]]))) {
727 mis_marker[j]=IS_IN_F;
728 }
729 }
730 }
731 }
732
733
734 /*This loop cannot be parallelized, as order matters here.*/
735 for (i=0;i<n;i++) {
736 if (mis_marker[i]==IS_IN_F) {
737 sum=0;
738 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
739 j=A_alpha->pattern->index[iptr];
740 if (mis_marker[j]==IS_IN_C) {
741 sum+=A_alpha->val[iptr];
742 }
743 }
744 if (ABS(sum)>=taw*ABS(A->val[diagptr[i]])) {
745 mis_marker[j]=IS_IN_FA;
746 }
747 }
748 }
749
750 /*This loop cannot be parallelized, as order matters here.*/
751 for (i=0;i<n;i++) {
752 if (mis_marker[i]!=IS_IN_C || mis_marker[i]!=IS_IN_FA) {
753 sum=0;
754 for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
755 j=A_alpha->pattern->index[iptr];
756 if (mis_marker[j]==IS_IN_C || mis_marker[j]==IS_IN_FA) {
757 sum+=A_alpha->val[iptr];
758 }
759 }
760 if (ABS(sum)>=delta*ABS(A->val[diagptr[i]])) {
761 mis_marker[j]=IS_IN_FB;
762 }
763 else {
764 mis_marker[j]=IS_IN_C;
765 }
766
767 }
768 }
769
770
771 /* swap to TRUE/FALSE in mis_marker */
772 #pragma omp parallel for private(i) schedule(static)
773 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_FA || mis_marker[i]==IS_IN_FB);
774
775 MEMFREE(diagptr);
776 MEMFREE(rsum);
777 Paso_SparseMatrix_free(A_alpha);
778 }
779
780
781 void Paso_Pattern_Standard(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
782 {
783 dim_t i,n,j,k;
784 index_t iptr,jptr;
785 /*index_t *index,*where_p;*/
786 double threshold,max_offdiagonal;
787 dim_t *lambda; /*mesure of importance */
788 dim_t maxlambda=0;
789 index_t index_maxlambda=0;
790 double time0=0;
791 bool_t verbose=0;
792
793 Paso_Pattern *S=NULL;
794 Paso_Pattern *ST=NULL;
795 Paso_IndexList* index_list=NULL;
796
797 /*dim_t lk;*/
798
799 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
800 if (! Paso_checkPtr(index_list)) {
801 #pragma omp parallel for private(i) schedule(static)
802 for(i=0;i<A->pattern->numOutput;++i) {
803 index_list[i].extension=NULL;
804 index_list[i].n=0;
805 }
806 }
807
808
809 n=A->numRows;
810 if (A->pattern->type & PATTERN_FORMAT_SYM) {
811 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
812 return;
813 }
814
815 time0=Paso_timer();
816 k=0;
817 /*S_i={j \in N_i; i strongly coupled to j}*/
818 #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)
819 for (i=0;i<n;++i) {
820 if(mis_marker[i]==IS_AVAILABLE) {
821 max_offdiagonal = DBL_MIN;
822 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
823 j=A->pattern->index[iptr];
824 if( j != i){
825 max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
826 }
827 }
828 threshold = theta*max_offdiagonal;
829 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
830 j=A->pattern->index[iptr];
831 if(ABS(A->val[iptr])>threshold && i!=j) {
832 Paso_IndexList_insertIndex(&(index_list[i]),j);
833 }
834 }
835 }
836 }
837
838 S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
839 ST=Paso_Pattern_getTranspose(S);
840
841 time0=Paso_timer()-time0;
842 if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
843
844 lambda=TMPMEMALLOC(n,dim_t);
845
846 #pragma omp parallel for private(i) schedule(static)
847 for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }
848
849 /*S_i={j \in N_i; i strongly coupled to j}*/
850
851 /*
852 #pragma omp parallel for private(i,iptr,lk) schedule(static)
853 for (i=0;i<n;++i) {
854 lk=0;
855 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
856 if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )
857 lk++;
858 }
859 #pragma omp critical
860 k+=lk;
861 if(k==0) {
862 mis_marker[i]=IS_IN_F;
863 }
864 }
865 */
866
867
868 k=0;
869 maxlambda=0;
870
871 time0=Paso_timer();
872
873 for (i=0;i<n;++i) {
874 if(mis_marker[i]==IS_AVAILABLE) {
875 lambda[i]=how_many(k,ST,FALSE); /* if every row is available then k and i are the same.*/
876 /*lambda[i]=how_many(i,S,TRUE);*/
877 /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/
878 k++;
879 if(maxlambda<lambda[i]) {
880 maxlambda=lambda[i];
881 index_maxlambda=i;
882 }
883 }
884 }
885
886 k=0;
887 time0=Paso_timer()-time0;
888 if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);
889
890 time0=Paso_timer();
891
892 /*Paso_Pattern_getReport(n,mis_marker);*/
893
894 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
895
896 if(index_maxlambda<0) {
897 break;
898 }
899
900 i=index_maxlambda;
901 if(mis_marker[i]==IS_AVAILABLE) {
902 mis_marker[index_maxlambda]=IS_IN_C;
903 lambda[index_maxlambda]=IS_NOT_AVAILABLE;
904 for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
905 j=ST->index[iptr];
906 if(mis_marker[j]==IS_AVAILABLE) {
907 mis_marker[j]=IS_IN_F;
908 lambda[j]=IS_NOT_AVAILABLE;
909 for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
910 k=S->index[jptr];
911 if(mis_marker[k]==IS_AVAILABLE) {
912 lambda[k]++;
913 }
914 }
915 }
916 }
917 for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
918 j=S->index[iptr];
919 if(mis_marker[j]==IS_AVAILABLE) {
920 lambda[j]--;
921 }
922 }
923 }
924
925 /* Used when transpose of S is not available */
926 /*
927 for (i=0;i<n;++i) {
928 if(mis_marker[i]==IS_AVAILABLE) {
929 if (i==index_maxlambda) {
930 mis_marker[index_maxlambda]=IS_IN_C;
931 lambda[index_maxlambda]=-1;
932 for (j=0;j<n;++j) {
933 if(mis_marker[j]==IS_AVAILABLE) {
934 index=&(S->index[S->ptr[j]]);
935 where_p=(index_t*)bsearch(&i,
936 index,
937 S->ptr[j + 1]-S->ptr[j],
938 sizeof(index_t),
939 Paso_comparIndex);
940 if (where_p!=NULL) {
941 mis_marker[j]=IS_IN_F;
942 lambda[j]=-1;
943 for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
944 k=S->index[iptr];
945 if(mis_marker[k]==IS_AVAILABLE) {
946 lambda[k]++;
947 }
948 }
949 }
950
951 }
952 }
953
954 }
955 }
956 }
957 */
958 index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);
959 }
960
961 time0=Paso_timer()-time0;
962 if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
963
964 /*Paso_Pattern_getReport(n,mis_marker);*/
965
966 #pragma omp parallel for private(i) schedule(static)
967 for (i=0;i<n;++i)
968 if(mis_marker[i]==IS_AVAILABLE) {
969 mis_marker[i]=IS_IN_F;
970 }
971
972 /*Paso_Pattern_getReport(n,mis_marker);*/
973
974 TMPMEMFREE(lambda);
975
976 /* clean up */
977 if (index_list!=NULL) {
978 #pragma omp parallel for private(i) schedule(static)
979 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
980 }
981
982 TMPMEMFREE(index_list);
983 Paso_Pattern_free(S);
984
985 /* swap to TRUE/FALSE in mis_marker */
986 #pragma omp parallel for private(i) schedule(static)
987 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
988
989 }
990
991 void Paso_Pattern_getReport(dim_t n,index_t* mis_marker) {
992
993 dim_t i,av,nav,f,c;
994
995 av=0;
996 nav=0;
997 f=0;
998 c=0;
999
1000 for (i=0;i<n;i++) {
1001 if(mis_marker[i]==IS_NOT_AVAILABLE) {
1002 nav++;
1003 }
1004 else if (mis_marker[i]==IS_AVAILABLE) {
1005 av++;
1006 }
1007 else if (mis_marker[i]==IS_IN_F) {
1008 f++;
1009 }
1010 else if (mis_marker[i]==IS_IN_C) {
1011 c++;
1012 }
1013 }
1014
1015 printf("Available=%d, Not Available = %d, In F = %d, In C = %d \n",av,nav,f,c);
1016
1017 }
1018
1019
1020 /*Used in Paso_Pattern_RS_MI*/
1021
1022 /*Computes how many connections unknown i have in S.
1023 bool_t transpose - TRUE if we want to compute how many strong connection of i in S^T, FALSE otherwise.
1024 Note that if we first transpose S and then call method on S^T then "transpose" should be set to FALSE.
1025 */
1026 dim_t how_many(dim_t i,Paso_Pattern * S, bool_t transpose) {
1027 dim_t j,n;
1028 dim_t total,ltotal;
1029 index_t *index,*where_p;
1030
1031 /*index_t iptr;*/
1032 total=0;
1033 ltotal=0;
1034
1035 n=S->numOutput;
1036
1037 if(transpose) {
1038 #pragma omp parallel
1039 {
1040 ltotal=0;
1041 #pragma omp for private(j,index,where_p,ltotal) schedule(static)
1042 for (j=0;j<n;++j) {
1043 index=&(S->index[S->ptr[j]]);
1044 where_p=(index_t*)bsearch(&i,
1045 index,
1046 S->ptr[j + 1]-S->ptr[j],
1047 sizeof(index_t),
1048 Paso_comparIndex);
1049 if (where_p!=NULL) {
1050 ltotal++;
1051 }
1052 }
1053 }
1054 #pragma omp critical
1055 {
1056 total+=ltotal;
1057 }
1058
1059 }
1060 else {
1061 total=S->ptr[i+1]-S->ptr[i];
1062 /*#pragma omp parallel for private(iptr) schedule(static)*/
1063 /*for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
1064 if(S->index[iptr]!=i && mis_marker[S->index[iptr]]==IS_AVAILABLE)
1065 total++;
1066 }*/
1067
1068 }
1069
1070 if (total==0) total=IS_NOT_AVAILABLE;
1071
1072 return total;
1073 }
1074
1075 dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {
1076 dim_t i;
1077 dim_t max=-1;
1078 dim_t argmax=-1;
1079
1080 #ifdef _OPENMP
1081 dim_t lmax=-1;
1082 dim_t li=-1;
1083 #pragma omp parallel private(i,lmax,li)
1084 {
1085 lmax=-1;
1086 li=-1;
1087 #pragma omp for schedule(static)
1088 for (i=0;i<n;++i) {
1089 if(lmax<lambda[i] && lambda[i]!=mask){
1090 lmax=lambda[i];
1091 li=i;
1092 }
1093 }
1094 #pragma omp critical
1095 {
1096 if (max<=lmax) {
1097 if(max==lmax && argmax>li)
1098 {
1099 argmax=li;
1100 }
1101 if (max < lmax)
1102 {
1103 max=lmax;
1104 argmax=li;
1105 }
1106 }
1107 }
1108 }
1109 #else
1110 for (i=0;i<n;++i) {
1111 if(max<lambda[i] && lambda[i]!=mask){
1112 max=lambda[i];
1113 argmax=i;
1114 }
1115 }
1116 #endif
1117
1118 return argmax;
1119 }
1120
1121
1122 Paso_Pattern* Paso_Pattern_getTranspose(Paso_Pattern* P){
1123
1124 Paso_Pattern *outpattern=NULL;
1125
1126 Paso_IndexList* index_list=NULL;
1127 dim_t C=P->numInput;
1128 dim_t F=P->numOutput-C;
1129 dim_t n=C+F;
1130 dim_t i,j;
1131 index_t iptr;
1132
1133 index_list=TMPMEMALLOC(C,Paso_IndexList);
1134 if (! Paso_checkPtr(index_list)) {
1135 #pragma omp parallel for private(i) schedule(static)
1136 for(i=0;i<C;++i) {
1137 index_list[i].extension=NULL;
1138 index_list[i].n=0;
1139 }
1140 }
1141
1142
1143 /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
1144 for (i=0;i<n;++i) {
1145 for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {
1146 j=P->index[iptr];
1147 Paso_IndexList_insertIndex(&(index_list[j]),i);
1148 }
1149 }
1150
1151 outpattern=Paso_IndexList_createPattern(0, C,index_list,0,n,0);
1152
1153 /* clean up */
1154 if (index_list!=NULL) {
1155 #pragma omp parallel for private(i) schedule(static)
1156 for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
1157 }
1158 TMPMEMFREE(index_list);
1159
1160 return outpattern;
1161 }
1162
1163
1164 /************** BLOCK COARSENENING *********************/
1165
1166 void Paso_Pattern_Standard_Block(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
1167 {
1168 dim_t i,n,j,k;
1169 index_t iptr,jptr;
1170 /*index_t *index,*where_p;*/
1171 double threshold,max_offdiagonal;
1172 dim_t *lambda; /*mesure of importance */
1173 dim_t maxlambda=0;
1174 index_t index_maxlambda=0;
1175 double time0=0;
1176 bool_t verbose=0;
1177 dim_t n_block=A->row_block_size;
1178
1179 double fnorm=0;
1180 dim_t bi;
1181
1182 Paso_Pattern *S=NULL;
1183 Paso_Pattern *ST=NULL;
1184 Paso_IndexList* index_list=NULL;
1185
1186 /*dim_t lk;*/
1187
1188 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
1189 if (! Paso_checkPtr(index_list)) {
1190 #pragma omp parallel for private(i) schedule(static)
1191 for(i=0;i<A->pattern->numOutput;++i) {
1192 index_list[i].extension=NULL;
1193 index_list[i].n=0;
1194 }
1195 }
1196
1197
1198 n=A->numRows;
1199 if (A->pattern->type & PATTERN_FORMAT_SYM) {
1200 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
1201 return;
1202 }
1203
1204 time0=Paso_timer();
1205 k=0;
1206 /*Paso_Pattern_getReport(n,mis_marker);*/
1207 /*printf("Blocks %d %d\n",n_block,A->len);*/
1208
1209 /*S_i={j \in N_i; i strongly coupled to j}*/
1210 #pragma omp parallel for private(i,bi,fnorm,iptr,max_offdiagonal,threshold,j) schedule(static)
1211 for (i=0;i<n;++i) {
1212 if(mis_marker[i]==IS_AVAILABLE) {
1213 max_offdiagonal = DBL_MIN;
1214 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1215 j=A->pattern->index[iptr];
1216 if( j != i){
1217 fnorm=0;
1218 for(bi=0;bi<n_block*n_block;++bi)
1219 {
1220 fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
1221 }
1222 fnorm=sqrt(fnorm);
1223 max_offdiagonal = MAX(max_offdiagonal,fnorm);
1224 }
1225 }
1226 threshold = theta*max_offdiagonal;
1227 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1228 j=A->pattern->index[iptr];
1229 if( j != i){
1230 fnorm=0;
1231 for(bi=0;bi<n_block*n_block;++bi)
1232 {
1233 fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
1234 }
1235 fnorm=sqrt(fnorm);
1236 if(fnorm>=threshold) {
1237 Paso_IndexList_insertIndex(&(index_list[i]),j);
1238 }
1239 }
1240 }
1241 }
1242 }
1243
1244 S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
1245 ST=Paso_Pattern_getTranspose(S);
1246
1247 /*printf("Patterns len %d %d\n",S->len,ST->len);*/
1248
1249 time0=Paso_timer()-time0;
1250 if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
1251
1252 lambda=TMPMEMALLOC(n,dim_t);
1253
1254 #pragma omp parallel for private(i) schedule(static)
1255 for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }
1256
1257 /*S_i={j \in N_i; i strongly coupled to j}*/
1258
1259 /*
1260 #pragma omp parallel for private(i,iptr,lk) schedule(static)
1261 for (i=0;i<n;++i) {
1262 lk=0;
1263 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1264 if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )
1265 lk++;
1266 }
1267 #pragma omp critical
1268 k+=lk;
1269 if(k==0) {
1270 mis_marker[i]=IS_IN_F;
1271 }
1272 }
1273 */
1274
1275
1276 k=0;
1277 maxlambda=0;
1278
1279 time0=Paso_timer();
1280
1281 for (i=0;i<n;++i) {
1282 if(mis_marker[i]==IS_AVAILABLE) {
1283 lambda[i]=how_many(k,ST,FALSE); /* if every row is available then k and i are the same.*/
1284 /*lambda[i]=how_many(i,S,TRUE);*/
1285 /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/
1286 k++;
1287 if(maxlambda<lambda[i]) {
1288 maxlambda=lambda[i];
1289 index_maxlambda=i;
1290 }
1291 }
1292 }
1293
1294 k=0;
1295 time0=Paso_timer()-time0;
1296 if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);
1297
1298 time0=Paso_timer();
1299
1300 /*Paso_Pattern_getReport(n,mis_marker);*/
1301
1302 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
1303
1304 if(index_maxlambda<0) {
1305 break;
1306 }
1307
1308 i=index_maxlambda;
1309 if(mis_marker[i]==IS_AVAILABLE) {
1310 mis_marker[index_maxlambda]=IS_IN_C;
1311 lambda[index_maxlambda]=IS_NOT_AVAILABLE;
1312 for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
1313 j=ST->index[iptr];
1314 if(mis_marker[j]==IS_AVAILABLE) {
1315 mis_marker[j]=IS_IN_F;
1316 lambda[j]=IS_NOT_AVAILABLE;
1317 for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
1318 k=S->index[jptr];
1319 if(mis_marker[k]==IS_AVAILABLE) {
1320 lambda[k]++;
1321 }
1322 }
1323 }
1324 }
1325 for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
1326 j=S->index[iptr];
1327 if(mis_marker[j]==IS_AVAILABLE) {
1328 lambda[j]--;
1329 }
1330 }
1331 }
1332
1333 /* Used when transpose of S is not available */
1334 /*
1335 for (i=0;i<n;++i) {
1336 if(mis_marker[i]==IS_AVAILABLE) {
1337 if (i==index_maxlambda) {
1338 mis_marker[index_maxlambda]=IS_IN_C;
1339 lambda[index_maxlambda]=-1;
1340 for (j=0;j<n;++j) {
1341 if(mis_marker[j]==IS_AVAILABLE) {
1342 index=&(S->index[S->ptr[j]]);
1343 where_p=(index_t*)bsearch(&i,
1344 index,
1345 S->ptr[j + 1]-S->ptr[j],
1346 sizeof(index_t),
1347 Paso_comparIndex);
1348 if (where_p!=NULL) {
1349 mis_marker[j]=IS_IN_F;
1350 lambda[j]=-1;
1351 for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
1352 k=S->index[iptr];
1353 if(mis_marker[k]==IS_AVAILABLE) {
1354 lambda[k]++;
1355 }
1356 }
1357 }
1358
1359 }
1360 }
1361
1362 }
1363 }
1364 }
1365 */
1366 index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);
1367 }
1368
1369 time0=Paso_timer()-time0;
1370 if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
1371
1372 /*Paso_Pattern_getReport(n,mis_marker);*/
1373
1374 #pragma omp parallel for private(i) schedule(static)
1375 for (i=0;i<n;++i)
1376 if(mis_marker[i]==IS_AVAILABLE) {
1377 mis_marker[i]=IS_IN_F;
1378 }
1379
1380 /*Paso_Pattern_getReport(n,mis_marker);*/
1381
1382 TMPMEMFREE(lambda);
1383
1384 /* clean up */
1385 if (index_list!=NULL) {
1386 #pragma omp parallel for private(i) schedule(static)
1387 for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
1388 }
1389
1390 TMPMEMFREE(index_list);
1391 Paso_Pattern_free(S);
1392
1393 /* swap to TRUE/FALSE in mis_marker */
1394 #pragma omp parallel for private(i) schedule(static)
1395 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
1396
1397 }
1398
1399
1400
1401 #undef IS_AVAILABLE
1402 #undef IS_NOT_AVAILABLE
1403 #undef IS_IN_F
1404 #undef IS_IN_C
1405
1406 #undef IS_UNDECIDED
1407 #undef IS_STRONG
1408 #undef IS_WEAK
1409
1410 #undef IS_IN_FB
1411 #undef IS_IN_FB
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421

  ViewVC Help
Powered by ViewVC 1.1.26