/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2760 by artak, Thu Nov 19 05:22:45 2009 UTC
# Line 36  Line 36 
36  /***************************************************************/  /***************************************************************/
37    
38  #define IS_AVAILABLE -1  #define IS_AVAILABLE -1
39  #define IS_IN_SET -3   /* Week connection */  #define IS_IN_F -3   /* in F (strong) */
40  #define IS_REMOVED -4  /* strong */  #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 sum;*/    /*double sum;*/
# Line 56  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 74  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 92  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;             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 106  void Paso_Pattern_coup(Paso_SparseMatrix Line 113  void Paso_Pattern_coup(Paso_SparseMatrix
113                  }                  }
114                }                }
115             }             }
116             if (passed) mis_marker[i]=IS_IN_SET;             if (passed) mis_marker[i]=IS_IN_C;
117          }          }
118      }      }
     /* 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];  
            }  
            if(ABS(sum)<1.e-25)  
              mis_marker[i]=IS_IN_SET;  
         }  
     }  
     */  
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 141  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 162  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,min_offdiagonal,threshold,j) schedule(static)      /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) 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);*/                  Paso_IndexList_insertIndex(&(index_list[j]),i);
                 }  
173              }              }
174          }          }
175         }         }
# Line 197  void Paso_Pattern_RS(Paso_SparseMatrix* Line 187  void Paso_Pattern_RS(Paso_SparseMatrix*
187    
188      /*Paso_Pattern_mis(out,mis_marker);*/      /*Paso_Pattern_mis(out,mis_marker);*/
189      Paso_Pattern_greedy(out,mis_marker);      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 228  void Paso_Pattern_Aggregiation(Paso_Spar Line 219  void Paso_Pattern_Aggregiation(Paso_Spar
219    }    }
220    
221    
222      #pragma omp parallel for private(i,iptr) reduction(+:diag) schedule(static)      #pragma omp parallel for private(i,iptr,diag) schedule(static)
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) {
# Line 244  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 271  void Paso_Pattern_Aggregiation(Paso_Spar Line 259  void Paso_Pattern_Aggregiation(Paso_Spar
259            
260      /*Paso_Pattern_mis(out,mis_marker);*/      /*Paso_Pattern_mis(out,mis_marker);*/
261      Paso_Pattern_greedy(out,mis_marker);      Paso_Pattern_greedy(out,mis_marker);
262        Paso_Pattern_free(out);
263    
264  }  }
265    
# Line 288  void Paso_Pattern_greedy(Paso_Pattern* p Line 277  void Paso_Pattern_greedy(Paso_Pattern* p
277      return;      return;
278    }    }
279        
    /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */  
280     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
281     for (i=0;i<n;++i)     for (i=0;i<n;++i)
282          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
283                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_IN_C;
284    
285    
286      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
287        if (mis_marker[i]==IS_IN_SET) {        if (mis_marker[i]==IS_IN_C) {
288          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
289               j=pattern->index[iptr];               j=pattern->index[iptr];
290               mis_marker[j]=IS_REMOVED;               mis_marker[j]=IS_IN_F;
291          }          }
292        }        }
293      }      }
# Line 307  void Paso_Pattern_greedy(Paso_Pattern* p Line 295  void Paso_Pattern_greedy(Paso_Pattern* p
295            
296            
297      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
298          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_IN_F) {
299             passed=TRUE;             passed=TRUE;
300             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
301                j=pattern->index[iptr];                j=pattern->index[iptr];
302                  if (mis_marker[j]==IS_REMOVED) {                  if (mis_marker[j]==IS_IN_F) {
303                      passed=TRUE;                      passed=TRUE;
304                  }                  }
305                  else {                  else {
# Line 319  void Paso_Pattern_greedy(Paso_Pattern* p Line 307  void Paso_Pattern_greedy(Paso_Pattern* p
307                      break;                      break;
308                  }                  }
309                }                }
310               if (passed) mis_marker[i]=IS_IN_C;
311             }             }
            if (passed) mis_marker[i]=IS_IN_SET;  
312          }          }
313    
314       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
315       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
316       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);
317            
318  }  }
319    
# Line 355  void Paso_Pattern_greedy_color(Paso_Patt Line 343  void Paso_Pattern_greedy_color(Paso_Patt
343     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
344     for (i=0;i<n;++i)     for (i=0;i<n;++i)
345          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
346                      mis_marker[i]=IS_IN_SET;                      mis_marker[i]=IS_IN_F;
347    
348     #pragma omp barrier     #pragma omp barrier
349     for (color=0;color<num_colors;++color) {     for (color=0;color<num_colors;++color) {
350      #pragma omp parallel for schedule(static) private(i,iptr,j)      #pragma omp parallel for schedule(static) private(i,iptr,j)
351      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
352       if (colorOf[i]==color) {         if (colorOf[i]==color) {  
353        if (mis_marker[i]==IS_IN_SET) {        if (mis_marker[i]==IS_IN_F) {
354          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {          for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
355               j=pattern->index[iptr];               j=pattern->index[iptr];
356               if (colorOf[j]<color)               if (colorOf[j]<color)
357                mis_marker[j]=IS_REMOVED;                mis_marker[j]=IS_IN_C;
358          }          }
359        }        }
360       }       }
# Line 379  void Paso_Pattern_greedy_color(Paso_Patt Line 367  void Paso_Pattern_greedy_color(Paso_Patt
367     #pragma omp parallel for schedule(static) private(i,iptr,j)     #pragma omp parallel for schedule(static) private(i,iptr,j)
368      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
369        if (colorOf[i]==color) {          if (colorOf[i]==color) {  
370          if (mis_marker[i]==IS_REMOVED) {          if (mis_marker[i]==IS_IN_C) {
371             passed=TRUE;             passed=TRUE;
372             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {             for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
373                j=pattern->index[iptr];                j=pattern->index[iptr];
374                 if (colorOf[j]<color && passed) {                 if (colorOf[j]<color && passed) {
375                  if (mis_marker[j]==IS_REMOVED) {                  if (mis_marker[j]==IS_IN_C) {
376                      passed=TRUE;                      passed=TRUE;
377                  }                  }
378                  else {                  else {
# Line 393  void Paso_Pattern_greedy_color(Paso_Patt Line 381  void Paso_Pattern_greedy_color(Paso_Patt
381                  }                  }
382                }                }
383             }             }
384               if (passed) mis_marker[i]=IS_IN_F;
385             }             }
386             if (passed) mis_marker[i]=IS_IN_SET;            
387          }          }
388      }      }
389     }     }
390    
391       /* swap to TRUE/FALSE in mis_marker */       /* swap to TRUE/FALSE in mis_marker */
392       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
393       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);
394            
395      MEMFREE(colorOf);      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             #pragma omp parallel for private(i,j,k) schedule(static)
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             #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)
465            for (i=0;i<n;++i) {
466                if(AvADJ[i]) {
467                    rsum=0;
468                    for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
469                        k=A->pattern->index[iptr];
470                        if(k!=i && mis_marker[k]!=IS_IN_C ) {
471                          rsum+=ABS(A->val[iptr]);    
472                        }
473                        if(j==i) {
474                            diag=ABS(A->val[iptr]);
475                        }
476                    }
477                    theta[i]=diag/rsum;
478                    if(theta[i]>threshold) {
479                       mis_marker[i]=IS_IN_F;
480                    }
481                }
482            }
483            
484            
485        }
486    
487         /* swap to TRUE/FALSE in mis_marker */
488         #pragma omp parallel for private(i) schedule(static)
489         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
490        
491         MEMFREE(AvADJ);
492         MEMFREE(theta);
493    }
494    
495    
496    void Paso_Pattern_YS_plus(Paso_SparseMatrix* A, index_t* mis_marker, double alpha, double taw, double delta) {
497    
498      dim_t i,j;
499      /*double sum;*/
500      index_t iptr,*index,*where_p,*diagptr;
501      double *rsum;
502      double sum;
503      dim_t n=A->numRows;
504      Paso_SparseMatrix* A_alpha;
505      diagptr=MEMALLOC(n,index_t);
506      rsum=MEMALLOC(n,double);
507    
508      if (A->pattern->type & PATTERN_FORMAT_SYM) {
509        Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
510        return;
511      }
512      
513        A_alpha=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern,1,1, FALSE);
514        #pragma omp parallel for private(i) schedule(static)
515        for (i=0;i<A->len;++i) {
516             A_alpha->val[i]=A->val[i];
517        }
518        
519    
520       #pragma omp parallel for private(i) schedule(static)
521       for (i=0;i<n;++i)
522            if(mis_marker[i]==IS_AVAILABLE)
523                        mis_marker[i]=IS_IN_C;
524    
525        /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/
526        for (i=0;i<n;++i) {
527             diagptr[i]=A->pattern->ptr[i];
528             index=&(A->pattern->index[A->pattern->ptr[i]]);
529             where_p=(index_t*)bsearch(&i,
530                                    index,
531                                    A->pattern->ptr[i + 1]-A->pattern->ptr[i],
532                                    sizeof(index_t),
533                                    Paso_comparIndex);
534            if (where_p==NULL) {
535                Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
536            } else {
537                    diagptr[i]+=(index_t)(where_p-index);
538            }
539        }
540        
541    
542        j=0;
543        for (i=0;i<n;++i) {
544            for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
545                j=A_alpha->pattern->index[iptr];
546                if(i==j) {
547                    A_alpha->val[iptr]=0;
548                }
549                else {
550                    if( !(ABS(A_alpha->val[iptr])<alpha*MIN(ABS(A_alpha->val[diagptr[i]]),ABS(A_alpha->val[diagptr[j]])) || A_alpha->val[iptr]*A_alpha->val[diagptr[i]]>0 || A_alpha->val[iptr]*A_alpha->val[diagptr[i]]<0) ) {
551                    A_alpha->val[iptr]=0;
552                    }
553                }
554                
555            }
556        }
557        
558        
559        for (i=0;i<n;++i) {
560            rsum[i]=0;
561            for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
562             rsum[i]+=A_alpha->val[iptr];  
563            }
564        }
565        
566        #pragma omp parallel for private(i,j) schedule(static)
567        for (i=0;i<n;++i) {
568            for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
569                j=A_alpha->pattern->index[iptr];
570                if (i==j) {
571                    A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr]+rsum[i];  
572                }
573                else {
574                    A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr];
575                }
576            
577            }
578        }
579    
580        /*This loop cannot be parallelized, as order matters here.*/
581        for (i=0;i<n;++i) {
582          if (mis_marker[i]==IS_IN_C) {
583            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
584                 j=A->pattern->index[iptr];
585                 if (j!=i && ABS(A->val[iptr])>=alpha*MIN(ABS(A->val[diagptr[i]]),ABS(A->val[diagptr[j]]))) {
586                    mis_marker[j]=IS_IN_F;
587                 }
588            }
589          }
590        }
591        
592          
593          /*This loop cannot be parallelized, as order matters here.*/
594        for (i=0;i<n;i++) {
595            if (mis_marker[i]==IS_IN_F) {
596                sum=0;
597               for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
598                  j=A_alpha->pattern->index[iptr];
599                  if (mis_marker[j]==IS_IN_C) {
600                    sum+=A_alpha->val[iptr];
601                  }
602               }
603               if (ABS(sum)>=taw*ABS(A->val[diagptr[i]])) {
604                   mis_marker[j]=IS_IN_FA;
605                }
606            }
607        }
608        
609          /*This loop cannot be parallelized, as order matters here.*/
610        for (i=0;i<n;i++) {
611            if (mis_marker[i]!=IS_IN_C || mis_marker[i]!=IS_IN_FA) {
612               sum=0;
613               for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {
614                  j=A_alpha->pattern->index[iptr];
615                  if (mis_marker[j]==IS_IN_C || mis_marker[j]==IS_IN_FA) {
616                    sum+=A_alpha->val[iptr];
617                  }
618               }
619               if (ABS(sum)>=delta*ABS(A->val[diagptr[i]])) {
620                   mis_marker[j]=IS_IN_FB;
621                }
622                else {
623                    mis_marker[j]=IS_IN_C;
624                }
625                
626            }
627        }
628        
629      
630         /* swap to TRUE/FALSE in mis_marker */
631         #pragma omp parallel for private(i) schedule(static)
632         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_FA || mis_marker[i]==IS_IN_FB);
633        
634         MEMFREE(diagptr);
635         MEMFREE(rsum);
636         Paso_SparseMatrix_free(A_alpha);
637    }
638    
639    
640    void Paso_Pattern_RS_MI(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
641    {
642      dim_t i,n,j,k;
643      index_t iptr,*index,*where_p;
644      double threshold,max_offdiagonal;
645      dim_t *lambda;   /*mesure of importance */
646      /*bool_t breakloop=FALSE;*/
647      dim_t maxlambda=0;
648      index_t index_maxlambda=0;
649      
650      Paso_Pattern *S=NULL;
651      Paso_IndexList* index_list=NULL;
652    
653     index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
654       if (! Paso_checkPtr(index_list)) {
655            #pragma omp parallel for private(i) schedule(static)
656            for(i=0;i<A->pattern->numOutput;++i) {
657                 index_list[i].extension=NULL;
658                 index_list[i].n=0;
659            }
660        }
661      
662      
663      n=A->numRows;
664      if (A->pattern->type & PATTERN_FORMAT_SYM) {
665        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
666        return;
667      }
668      
669        /*S_i={j \in N_i; i strongly coupled to j}*/
670        /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
671        for (i=0;i<n;++i) {
672          if(mis_marker[i]==IS_AVAILABLE) {
673            max_offdiagonal = DBL_MIN;
674            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
675                if(A->pattern->index[iptr] != i){
676                    if(A->val[iptr]<0) {
677                      max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
678                    }
679                }
680            }
681            
682            threshold = theta*max_offdiagonal;
683            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
684                j=A->pattern->index[iptr];
685                if((-A->val[iptr])>=threshold) {
686                    Paso_IndexList_insertIndex(&(index_list[i]),j);
687                }
688            }
689           }
690          }
691        
692      
693        S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
694      
695      lambda=TMPMEMALLOC(n,dim_t);
696      
697      for (i=0;i<n;++i) {
698         lambda[i]=-1;
699       }
700      
701      /*S_i={j \in N_i; i strongly coupled to j}*/
702    
703      k=0;
704      maxlambda=0;
705      
706        for (i=0;i<n;++i) {
707          if(mis_marker[i]==IS_AVAILABLE) {
708            lambda[i]=how_many(i,S,TRUE);
709            /*printf("lambda[%d]=%d, ",i,lambda[i]);*/
710            if(maxlambda<lambda[i]) {
711                maxlambda=lambda[i];
712                index_maxlambda=i;
713            }
714          }
715        }
716      
717      while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
718        
719        index_maxlambda=arg_max(n,lambda, -1);
720        if(index_maxlambda<0) {
721            break;
722        }
723    
724        for (i=0;i<n;++i) {
725            if(mis_marker[i]==IS_AVAILABLE) {
726                if (i==index_maxlambda) {
727                    mis_marker[index_maxlambda]=IS_IN_C;
728                    lambda[index_maxlambda]=-1;
729                    for (j=0;j<n;++j) {
730                        if(mis_marker[j]==IS_AVAILABLE) {
731                            index=&(S->index[S->ptr[j]]);
732                            where_p=(index_t*)bsearch(&i,
733                                            index,
734                                            S->ptr[j + 1]-S->ptr[j],
735                                            sizeof(index_t),
736                                            Paso_comparIndex);
737                            if (where_p!=NULL) {
738                                mis_marker[j]=IS_IN_F;
739                                lambda[j]=-1;
740                                for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
741                                    k=S->index[iptr];
742                                    if(mis_marker[k]==IS_AVAILABLE) {
743                                       lambda[k]++;
744                                    }
745                                }
746                            }
747                            
748                        }
749                    }
750                    
751                }
752            }
753        }
754        
755      }
756      
757      for (i=0;i<n;++i)
758          if(mis_marker[i]==IS_AVAILABLE)
759            mis_marker[i]=IS_IN_F;
760    
761        /*update lambdas*/
762        /*for (i=0;i<n;++i) {
763          if(mis_marker[i]==IS_AVAILABLE) {
764            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);
765            if(maxlambda<lambda[i]) {
766                maxlambda=lambda[i];
767                index_maxlambda=i;
768            }
769          }
770          if(lambda[i]==0) {
771            breakloop=TRUE;
772            break;
773          }
774        }
775        if(breakloop) {
776            break;
777        }
778        
779        for (i=0;i<n;++i) {
780            if(mis_marker[i]==IS_AVAILABLE) {
781                mis_marker[index_maxlambda]=IS_IN_C;
782            }
783            
784            for (j=0;j<n;++j) {
785                if(S_T_[i][j]=IS_STRONG && mis_marker[i]==IS_AVAILABLE) {
786                    mis_marker[j]==IS__IN_F;
787                }
788            }
789        }
790        
791        }
792        */
793    
794        TMPMEMFREE(lambda);
795        
796         /* clean up */
797        if (index_list!=NULL) {
798            #pragma omp parallel for private(i) schedule(static)
799            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
800         }
801    
802        TMPMEMFREE(index_list);
803        Paso_Pattern_free(S);
804        
805    
806        /* swap to TRUE/FALSE in mis_marker */
807         #pragma omp parallel for private(i) schedule(static)
808         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
809        
810    }
811    
812    /*Used in Paso_Pattern_RS_MI*/
813    
814    dim_t how_many(dim_t i,Paso_Pattern * S, bool_t transpose) {
815        dim_t j,n;
816        dim_t total;
817        index_t iptr,*index,*where_p;
818        total=0;
819        
820        n=S->numOutput;
821        
822        if(transpose) {
823            for (j=0;j<n;++j) {
824                index=&(S->index[S->ptr[j]]);
825                where_p=(index_t*)bsearch(&i,
826                                        index,
827                                        S->ptr[j + 1]-S->ptr[j],
828                                        sizeof(index_t),
829                                        Paso_comparIndex);
830                if (where_p!=NULL) {
831                    total++;
832                }
833            }
834        }
835        else {
836            for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
837                    total++;
838            }
839            
840        }
841    return total;
842    }
843    
844    dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {
845        dim_t i;
846        dim_t max=0;
847        dim_t argmax=-1;
848        for (i=0;i<n;++i) {
849            if(max<lambda[i] && lambda[i]!=mask){
850              argmax=i;
851              max=lambda[i];
852            }
853        }
854        return argmax;
855    }
856    
857    
858    /*dim_t how_many(dim_t n,dim_t* S_i, int value1, dim_t* addedSet, int value2) {
859        dim_t j;
860        dim_t total;
861        total=0;
862        for (j=0;j<n;++j) {
863            if(S_i[j]==value1 && addedSet[j]==value2)
864              total++;
865            }
866    return total;
867    }
868    */
869    
870  #undef IS_AVAILABLE  #undef IS_AVAILABLE
871  #undef IS_IN_SET  #undef IS_IN_F
872  #undef IS_REMOVED  #undef IS_IN_C
873    
874    #undef IS_UNDECIDED
875    #undef IS_STRONG
876    #undef IS_WEAK
877    
878    #undef IS_IN_FB
879    #undef IS_IN_FB
880    
881    
882    
883    
884    
885    
886    
887    
888    
889    

Legend:
Removed from v.2548  
changed lines
  Added in v.2760

  ViewVC Help
Powered by ViewVC 1.1.26