/[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 2760 by artak, Thu Nov 19 05:22:45 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             #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.2381  
changed lines
  Added in v.2760

  ViewVC Help
Powered by ViewVC 1.1.26