/[escript]/trunk/paso/src/Pattern_coupling.c
ViewVC logotype

Annotation of /trunk/paso/src/Pattern_coupling.c

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26