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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2381 by artak, Tue Apr 14 03:46:59 2009 UTC revision 2764 by artak, Thu Nov 19 23:48:05 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 36  Line 36 
36  /***************************************************************/  /***************************************************************/
37    
38  #define IS_AVAILABLE -1  #define IS_AVAILABLE -1
39  #define IS_IN_SET -3  #define IS_IN_F -3   /* in F (strong) */
40  #define IS_REMOVED -4  #define IS_IN_C -4  /* in C (weak) */
41                        
42    
43  void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  #define IS_UNDECIDED -1                      
44    #define IS_STRONG -2
45    #define IS_WEAK -3
46    
47    
48    #define IS_IN_FA -5  /* test */
49    #define IS_IN_FB -6  /* test */
50    
51    void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
52    
53    dim_t i,j;    dim_t i,j;
54    /*double threshold=0.05;*/    /*double sum;*/
   double sum;  
55    index_t iptr,*index,*where_p,*diagptr;    index_t iptr,*index,*where_p,*diagptr;
56    bool_t passed=FALSE;    bool_t passed=FALSE;
57    dim_t n=A->numRows;    dim_t n=A->numRows;
# Line 57  void Paso_Pattern_coup(Paso_SparseMatrix Line 65  void Paso_Pattern_coup(Paso_SparseMatrix
65     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
66     for (i=0;i<n;++i)     for (i=0;i<n;++i)
67          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
68                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_IN_C;
69    
70      #pragma omp parallel for private(i,index,where_p) schedule(static)      /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
71      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
72           diagptr[i]=A->pattern->ptr[i];           diagptr[i]=A->pattern->ptr[i];
73           index=&(A->pattern->index[diagptr[i]]);           index=&(A->pattern->index[A->pattern->ptr[i]]);
74           where_p=(index_t*)bsearch(&i,           where_p=(index_t*)bsearch(&i,
75                                  index,                                  index,
76                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],                                  A->pattern->ptr[i + 1]-A->pattern->ptr[i],
# Line 75  void Paso_Pattern_coup(Paso_SparseMatrix Line 83  void Paso_Pattern_coup(Paso_SparseMatrix
83          }          }
84      }      }
85            
   
   
86      /*This loop cannot be parallelized, as order matters here.*/      /*This loop cannot be parallelized, as order matters here.*/
87      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
88        if (mis_marker[i]==IS_IN_SET) {        if (mis_marker[i]==IS_IN_C) {
89          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
90               j=A->pattern->index[iptr];               j=A->pattern->index[iptr];
91               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
92                  mis_marker[j]=IS_REMOVED;                  mis_marker[j]=IS_IN_F;
93               }               }
94          }          }
95        }        }
# Line 93  void Paso_Pattern_coup(Paso_SparseMatrix Line 99  void Paso_Pattern_coup(Paso_SparseMatrix
99            
100        /*This loop cannot be parallelized, as order matters here.*/        /*This loop cannot be parallelized, as order matters here.*/
101      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
102          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_IN_F) {
103               passed=TRUE;
104             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
105                j=A->pattern->index[iptr];                j=A->pattern->index[iptr];
106                if (mis_marker[j]==IS_IN_SET) {                if (mis_marker[j]==IS_IN_C) {
107                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {                  if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
108                      passed=TRUE;                      passed=TRUE;
109                  }                  }
# Line 105  void Paso_Pattern_coup(Paso_SparseMatrix Line 112  void Paso_Pattern_coup(Paso_SparseMatrix
112                      break;                      break;
113                  }                  }
114                }                }
             }  
             if (passed) {  
                mis_marker[i]=IS_IN_SET;  
                passed=FALSE;  
             }  
         }  
     }  
   
     /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/  
     /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */  
     #pragma omp parallel for private(i,iptr,j,sum) schedule(static)  
     for (i=0;i<n;i++) {  
         if (mis_marker[i]==IS_REMOVED) {  
            sum=0;  
            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
              j=A->pattern->index[iptr];  
              if (mis_marker[j]==IS_REMOVED)  
                 sum+=A->val[iptr];  
115             }             }
116             if(ABS(sum)<1.e-25)             if (passed) mis_marker[i]=IS_IN_C;
              mis_marker[i]=IS_IN_SET;  
117          }          }
118      }      }
119    
   
120       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
121       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
122       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
123            
124       MEMFREE(diagptr);       MEMFREE(diagptr);
125  }  }
126    
127    
128  /*  /*
129   * Ruge-Stueben strength of connection mask.   * Ruge-Stueben strength of connection mask.
130   *   *
# Line 145  void Paso_Pattern_RS(Paso_SparseMatrix* Line 133  void Paso_Pattern_RS(Paso_SparseMatrix*
133  {  {
134    dim_t i,n,j;    dim_t i,n,j;
135    index_t iptr;    index_t iptr;
136    double threshold,min_offdiagonal;    double threshold,max_offdiagonal;
137        
138    Paso_Pattern *out=NULL;    Paso_Pattern *out=NULL;
139        
# Line 166  void Paso_Pattern_RS(Paso_SparseMatrix* Line 154  void Paso_Pattern_RS(Paso_SparseMatrix*
154      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
155      return;      return;
156    }    }
157        /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
     #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)  
158      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
159        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
160          min_offdiagonal = DBL_MAX;          max_offdiagonal = DBL_MIN;
161          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
162              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
163                  min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);                    max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
164              }              }
165          }          }
166            
167          threshold = theta*min_offdiagonal;          threshold = theta*max_offdiagonal;
168          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
169              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
170              if(A->val[iptr]<=threshold) {              if((-A->val[iptr])>=threshold) {
                if(j!=i) {  
171                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
172                  }                  Paso_IndexList_insertIndex(&(index_list[j]),i);
173              }              }
174          }          }
175         }         }
176        }        }
177            
178      
179      out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);      out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
180            
181       /* clean up */       /* clean up */
# Line 198  void Paso_Pattern_RS(Paso_SparseMatrix* Line 185  void Paso_Pattern_RS(Paso_SparseMatrix*
185       }       }
186      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
187    
188      Paso_Pattern_mis(out,mis_marker);      /*Paso_Pattern_mis(out,mis_marker);*/
189        Paso_Pattern_greedy(out,mis_marker);
190        Paso_Pattern_free(out);
191  }  }
192    
193  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
# Line 225  void Paso_Pattern_Aggregiation(Paso_Spar Line 214  void Paso_Pattern_Aggregiation(Paso_Spar
214      }      }
215            
216    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
217      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
218      return;      return;
219    }    }
220    
# Line 234  void Paso_Pattern_Aggregiation(Paso_Spar Line 223  void Paso_Pattern_Aggregiation(Paso_Spar
223        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
224          diag = 0;          diag = 0;
225          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
226              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] == i){
227                  diag+=A->val[iptr];                  diag+=A->val[iptr];
228              }              }
229          }          }
# Line 246  void Paso_Pattern_Aggregiation(Paso_Spar Line 235  void Paso_Pattern_Aggregiation(Paso_Spar
235       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
236         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
237          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
         val=0.;  
238          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
239              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
240              val=A->val[iptr];              val=A->val[iptr];
241              if(j!= i) {                if((val*val)>=(eps_Aii*diags[j])) {
               if(val*val>=eps_Aii * diags[j]) {  
242                 Paso_IndexList_insertIndex(&(index_list[i]),j);                 Paso_IndexList_insertIndex(&(index_list[i]),j);
243                }                }
             }  
244          }          }
245         }         }
246       }       }
# Line 270  void Paso_Pattern_Aggregiation(Paso_Spar Line 256  void Paso_Pattern_Aggregiation(Paso_Spar
256      TMPMEMFREE(index_list);      TMPMEMFREE(index_list);
257      MEMFREE(diags);      MEMFREE(diags);
258            
259      Paso_Pattern_mis(out,mis_marker);      
260        /*Paso_Pattern_mis(out,mis_marker);*/
261        Paso_Pattern_greedy(out,mis_marker);
262        Paso_Pattern_free(out);
263    
264    }
265    
266    /* Greedy algorithm */
267    void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
268    
269      dim_t i,j;
270      /*double sum;*/
271      index_t iptr;
272      bool_t passed=FALSE;
273      dim_t n=pattern->numOutput;
274    
275      if (pattern->type & PATTERN_FORMAT_SYM) {
276        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
277        return;
278      }
279      
280       #pragma omp parallel for private(i) schedule(static)
281       for (i=0;i<n;++i)
282            if(mis_marker[i]==IS_AVAILABLE)
283                        mis_marker[i]=IS_IN_C;
284    
285    
286        for (i=0;i<n;++i) {
287          if (mis_marker[i]==IS_IN_C) {
288            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
289                 j=pattern->index[iptr];
290                 mis_marker[j]=IS_IN_F;
291            }
292          }
293        }
294        
295        
296        
297        for (i=0;i<n;i++) {
298            if (mis_marker[i]==IS_IN_F) {
299               passed=TRUE;
300               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
301                  j=pattern->index[iptr];
302                    if (mis_marker[j]==IS_IN_F) {
303                        passed=TRUE;
304                    }
305                    else {
306                        passed=FALSE;
307                        break;
308                    }
309                  }
310               if (passed) mis_marker[i]=IS_IN_C;
311               }
312            }
313    
314         /* swap to TRUE/FALSE in mis_marker */
315         #pragma omp parallel for private(i) schedule(static)
316         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
317        
318    }
319    
320    
321    void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
322    
323      dim_t i,j;
324      /*double sum;*/
325      index_t iptr;
326      index_t num_colors;
327      index_t* colorOf;
328      register index_t color;
329      bool_t passed=FALSE;
330      dim_t n=pattern->numOutput;
331    
332      
333      colorOf=MEMALLOC(n,index_t);
334    
335      if (pattern->type & PATTERN_FORMAT_SYM) {
336        Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
337        return;
338      }
339      
340       Paso_Pattern_color(pattern,&num_colors,colorOf);
341      
342       /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
343       #pragma omp parallel for private(i) schedule(static)
344       for (i=0;i<n;++i)
345            if(mis_marker[i]==IS_AVAILABLE)
346                        mis_marker[i]=IS_IN_F;
347    
348       #pragma omp barrier
349       for (color=0;color<num_colors;++color) {
350        #pragma omp parallel for schedule(static) private(i,iptr,j)
351        for (i=0;i<n;++i) {
352         if (colorOf[i]==color) {  
353          if (mis_marker[i]==IS_IN_F) {
354            for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
355                 j=pattern->index[iptr];
356                 if (colorOf[j]<color)
357                  mis_marker[j]=IS_IN_C;
358            }
359          }
360         }
361        }
362       }
363        
364        
365       #pragma omp barrier
366       for (color=0;color<num_colors;++color) {
367       #pragma omp parallel for schedule(static) private(i,iptr,j)
368        for (i=0;i<n;i++) {
369          if (colorOf[i]==color) {  
370            if (mis_marker[i]==IS_IN_C) {
371               passed=TRUE;
372               for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
373                  j=pattern->index[iptr];
374                   if (colorOf[j]<color && passed) {
375                    if (mis_marker[j]==IS_IN_C) {
376                        passed=TRUE;
377                    }
378                    else {
379                        passed=FALSE;
380                        /*break;*/
381                    }
382                  }
383               }
384               if (passed) mis_marker[i]=IS_IN_F;
385               }
386              
387            }
388        }
389       }
390    
391         /* swap to TRUE/FALSE in mis_marker */
392         #pragma omp parallel for private(i) schedule(static)
393         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
394        
395        MEMFREE(colorOf);
396    }
397    
398    /*For testing */
399    void Paso_Pattern_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
400    
401      dim_t i,j=0,k;
402      double *theta;
403      index_t iptr;
404      dim_t n=A->numRows;
405      double rsum,diag=0;
406      index_t *AvADJ;
407      theta=MEMALLOC(n,double);
408      AvADJ=MEMALLOC(n,index_t);
409    
410    
411      
412    
413      if (A->pattern->type & PATTERN_FORMAT_SYM) {
414        Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
415        return;
416      }
417      
418    
419        #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
420        for (i=0;i<n;++i) {
421            rsum=0;
422            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
423                j=A->pattern->index[iptr];
424                if(j!=i) {
425                  rsum+=ABS(A->val[iptr]);    
426                }
427                else {
428                    diag=ABS(A->val[iptr]);
429                }
430            }
431            theta[i]=diag/rsum;
432            if(theta[i]>threshold) {
433                mis_marker[i]=IS_IN_F;
434            }
435        }
436        
437        while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
438             k=0;
439            
440             for (i=0;i<n;++i) {
441               if(mis_marker[i]==IS_AVAILABLE) {
442                    if(k==0) {
443                        j=i;
444                        k++;
445                    }
446                    if(theta[j]>theta[i]) {
447                        j=i;
448                    }
449                }
450             }
451             mis_marker[j]=IS_IN_C;
452            
453             for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {
454                k=A->pattern->index[iptr];
455                if(mis_marker[k]==IS_AVAILABLE) {
456                   AvADJ[k]=1;
457                }
458                else {
459                    AvADJ[k]=-1;
460                }
461                
462             }
463                
464            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    
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  #undef IS_AVAILABLE  #undef IS_AVAILABLE
870  #undef IS_IN_SET  #undef IS_IN_F
871  #undef IS_REMOVED  #undef IS_IN_C
872    
873    #undef IS_UNDECIDED
874    #undef IS_STRONG
875    #undef IS_WEAK
876    
877    #undef IS_IN_FB
878    #undef IS_IN_FB
879    
880    
881    
882    
883    
884    
885    
886    
887    
888    

Legend:
Removed from v.2381  
changed lines
  Added in v.2764

  ViewVC Help
Powered by ViewVC 1.1.26