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

  ViewVC Help
Powered by ViewVC 1.1.26