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

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

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

revision 3282 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3283 by gross, Mon Oct 18 22:39:28 2010 UTC
# Line 12  Line 12 
12  *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /*************************************************************
16    
17  /* Paso: Coarsening strategies                                */   Paso: Coarsening strategies (no MPI)  
18    
19  /**************************************************************/   the given matrix A is splitted in to the form
20    
21    
22           [ A_FF A_FC ]
23     A  =  [           ]
24           [ A_CF A_CC ]
25    
26     such that unknowns/equations in set C are weakly connected via A_CC and strongly connected
27     to at least one unknowns/equation in F via A_CF and A_FC. The unknowns/equations in C and F
28     are marked in the marker_F by FALSE and TRUE respectively.
29    
30     The weak/strong connection is controlled by coarsening_threshold.
31    
32     three strategies are implemented:
33        
34        a) YAIR_SHAPIRA (YS): |a_ij| >= theta * |a_ii|
35        b) Ruge-Stueben (RS): |a_ij| >= theta * max_(k<>i) |a_ik|
36        c) Aggregation :     |a_ij|^2 >= theta**2 * |a_ii||a_jj|
37    
38    where theta = coarsening_threshold/maximum_pattern_degree
39      
40      Remark:
41      
42           - a strong connection in YAIR_SHAPIRA is a strong connection in Aggregation
43      
44    */
45    /**************************************************************
46    
47  /* Author: artak@uq.edu.au                                */    Author: artak@uq.edu.au, l.gross@uq.edu.au                  
48    
49  /**************************************************************/  **************************************************************/
50    
51  #include "Coarsening.h"  #include "Coarsening.h"
52  #include "PasoUtil.h"  #include "PasoUtil.h"
# Line 28  Line 54 
54    
55    
56    
 void Paso_Coarsening_Local(index_t* mis_marker, Paso_SparseMatrix* A_p, double coarsening_threshold, const index_t coarsening_method)  
 {  
   
    const dim_t n=A_p->numRows;  
   /* const dim_t n_block=A_p->row_block_size; */  
    dim_t i;  
     
   
    /* identify independend set of rows/columns */  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i) mis_marker[i]=-1;  
   
    if (coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {  
       Paso_Coarsening_Local_YS(A_p,mis_marker,coarsening_threshold);  
    } else if (coarsening_method == PASO_RUGE_STUEBEN_COARSENING) {  
       Paso_Coarsening_Local_RS(A_p,mis_marker,coarsening_threshold);  
    } else if (coarsening_method == PASO_AGGREGATION_COARSENING) {  
       Paso_Coarsening_Local_Aggregiation(A_p,mis_marker,coarsening_threshold);  
    } else {  
       Paso_Coarsening_Local_Standard_Block(A_p,mis_marker,coarsening_threshold);  
    }  
 }  
   
 /***************************************************************/  
   
 #define IS_AVAILABLE -1  
 #define IS_NOT_AVAILABLE -2  
 #define IS_IN_F -3   /* in F (strong) */  
 #define IS_IN_C -4  /* in C (weak) */  
   
   
 #define IS_UNDECIDED -1                        
 #define IS_STRONG -2  
 #define IS_WEAK -3  
   
   
 #define IS_IN_FA -5  /* test */  
 #define IS_IN_FB -6  /* test */  
   
   
 /* YAIR_SHAPIRA */  
 void Paso_Coarsening_Local_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  
     
    dim_t i,j,bi;  
    /*double sum;*/  
    index_t iptr;  
    /*index_t *index,*where_p;*/  
    bool_t passed=FALSE;  
    dim_t n=A->numRows;  
    dim_t block_size=A->row_block_size;  
    double *diags;  
    double fnorm;  
    diags=MEMALLOC(n,double);  
     
     
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_coup: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i)  
       if(mis_marker[i]==IS_AVAILABLE)  
      mis_marker[i]=IS_IN_C;  
         
       #pragma omp parallel for private(i,j,fnorm,bi) schedule(static)  
       for (i=0;i<n;++i) {  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         j=A->pattern->index[iptr];  
         if(i==j) {  
            if(block_size==1) {  
           diags[i]=A->val[iptr];      
            } else {  
           fnorm=0;  
           for(bi=0;bi<block_size*block_size;++bi)  
           {  
              fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
           }  
           fnorm=sqrt(fnorm);  
           diags[i]=fnorm;  
            }  
             
         }  
      }  
       }  
         
       /*This loop cannot be parallelized, as order matters here.*/  
       for (i=0;i<n;++i) {  
      if (mis_marker[i]==IS_IN_C) {  
         for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
            j=A->pattern->index[iptr];  
            if(block_size==1) {  
           if (j!=i && ABS(A->val[iptr])>=threshold*ABS(diags[i])) {  
              mis_marker[j]=IS_IN_F;  
           }  
            } else {  
           if (j!=i) {  
              fnorm=0;  
              for(bi=0;bi<block_size*block_size;++bi)  
              {  
             fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
              }  
              fnorm=sqrt(fnorm);  
              if (fnorm>=threshold*diags[i]) {  
             mis_marker[j]=IS_IN_F;  
              }  
           }  
            }  
         }  
      }  
       }  
         
       /*This loop cannot be parallelized, as order matters here.*/  
       for (i=0;i<n;i++) {  
      if (mis_marker[i]==IS_IN_F) {  
         passed=TRUE;  
         for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
            j=A->pattern->index[iptr];  
            if (mis_marker[j]==IS_IN_C) {  
           if(block_size==1) {  
              if ((A->val[iptr]/diags[i])>=-threshold) {  
             passed=TRUE;  
              }  
              else {  
             passed=FALSE;  
             break;  
              }  
           } else {  
              fnorm=0;  
              for(bi=0;bi<block_size*block_size;++bi)  
              {  
             fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
              }  
              fnorm=sqrt(fnorm);  
              if ((fnorm/diags[i])<=threshold) {  
             passed=TRUE;  
              }  
              else {  
             passed=FALSE;  
             break;  
              }  
           }  
             
            }  
             
         }  
         if (passed) mis_marker[i]=IS_IN_C;  
      }  
       }  
         
       /* swap to TRUE/FALSE in mis_marker */  
       #pragma omp parallel for private(i) schedule(static)  
       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F) ? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;  
         
       MEMFREE(diags);  
 }  
   
   
 /*  
 * Ruge-Stueben strength of connection mask.  
 *  
 */  
 void Paso_Coarsening_Local_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
 {  
    dim_t i,n,j,bi;  
    index_t iptr;  
    double threshold,max_offdiagonal;  
    double fnorm;  
    dim_t block_size=A->row_block_size;  
     
    Paso_Pattern *out=NULL;  
     
    Paso_IndexList* index_list=NULL;  
     
    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) {  
      index_list[i].extension=NULL;  
      index_list[i].n=0;  
       }  
    }  
     
     
    n=A->numRows;  
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_RS: symmetric matrix pattern is not supported yet");  
       return;  
    }  
    /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j,fnorm,bi) schedule(static)*/  
    for (i=0;i<n;++i) {  
       if(mis_marker[i]==IS_AVAILABLE) {  
      max_offdiagonal = DBL_MIN;  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         if(A->pattern->index[iptr] != i){  
            if(block_size==1) {  
           max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));  
            } else {  
           fnorm=0;  
           for(bi=0;bi<block_size*block_size;++bi)  
           {  
              fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
           }  
           fnorm=sqrt(fnorm);  
           max_offdiagonal = MAX(max_offdiagonal,fnorm);  
            }  
         }  
      }  
       
      threshold = theta*max_offdiagonal;  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         j=A->pattern->index[iptr];  
         if(block_size==1) {  
            if(ABS(A->val[iptr])>=threshold && i!=j) {  
           Paso_IndexList_insertIndex(&(index_list[i]),j);  
           /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/  
            }  
         } else {  
            if (i!=j) {  
           fnorm=0;  
           for(bi=0;bi<block_size*block_size;++bi)  
           {  
              fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
           }  
           fnorm=sqrt(fnorm);  
           if(fnorm>=threshold) {  
              Paso_IndexList_insertIndex(&(index_list[i]),j);  
              /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/  
           }  
            }  
         }  
      }  
       }  
    }  
     
     
    out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);  
     
    /* clean up */  
    if (index_list!=NULL) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
    }  
    TMPMEMFREE(index_list);  
     
    /*Paso_Coarsening_Local_mis(out,mis_marker);*/  
    Paso_Coarsening_Local_greedy(out,mis_marker);  
    Paso_Pattern_free(out);  
 }  
   
 void Paso_Coarsening_Local_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
 {  
    dim_t i,j,n,bi;  
    index_t iptr;  
    double diag,eps_Aii,val;  
    double* diags;  
    double fnorm;  
    dim_t block_size=A->row_block_size;  
     
     
    Paso_Pattern *out=NULL;  
    Paso_IndexList* index_list=NULL;  
     
    n=A->numRows;    
    diags=MEMALLOC(n,double);  
     
    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) {  
      index_list[i].extension=NULL;  
      index_list[i].n=0;  
       }  
    }  
     
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_Aggregiation: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
     
    #pragma omp parallel for private(i,iptr,diag,fnorm,bi) schedule(static)  
    for (i=0;i<n;++i) {  
       diag = 0;  
       for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
      if(A->pattern->index[iptr] == i){  
         if(block_size==1) {  
            diag+=A->val[iptr];  
         } else {  
            fnorm=0;  
            for(bi=0;bi<block_size*block_size;++bi)  
            {  
           fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
            }  
            fnorm=sqrt(fnorm);  
            diag+=fnorm;  
         }  
           
      }  
       }  
       diags[i]=ABS(diag);  
    }  
     
     
    #pragma omp parallel for private(i,iptr,j,val,eps_Aii,fnorm,bi) schedule(static)  
    for (i=0;i<n;++i) {  
       if (mis_marker[i]==IS_AVAILABLE) {  
      eps_Aii = theta*theta*diags[i];  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         j=A->pattern->index[iptr];  
         if(block_size==1) {  
            val=A->val[iptr];  
         } else {  
            fnorm=0;  
            for(bi=0;bi<block_size*block_size;++bi)  
            {  
           fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];  
            }  
            fnorm=sqrt(fnorm);  
            val=fnorm;  
         }  
           
         if((val*val)>=(eps_Aii*diags[j])) {  
            Paso_IndexList_insertIndex(&(index_list[i]),j);  
         }  
      }  
       }  
    }  
     
    out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);  
     
    /* clean up */  
    if (index_list!=NULL) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
    }  
     
    TMPMEMFREE(index_list);  
    MEMFREE(diags);  
     
     
    /*Paso_Coarsening_Local_mis(out,mis_marker);*/  
    Paso_Coarsening_Local_greedy(out,mis_marker);  
    Paso_Pattern_free(out);  
     
 }  
   
 /* Greedy algorithm */  
 void Paso_Coarsening_Local_greedy(Paso_Pattern* pattern, index_t* mis_marker) {  
     
    dim_t i,j;  
    /*double sum;*/  
    index_t iptr;  
    bool_t passed=FALSE;  
    dim_t n=pattern->numOutput;  
     
    if (pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_greedy: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i)  
       if(mis_marker[i]==IS_AVAILABLE)  
      mis_marker[i]=IS_IN_C;  
         
         
       for (i=0;i<n;++i) {  
      if (mis_marker[i]==IS_IN_C) {  
         for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {  
            j=pattern->index[iptr];  
            mis_marker[j]=IS_IN_F;  
         }  
      }  
       }  
         
         
       for (i=0;i<n;i++) {  
      if (mis_marker[i]==IS_IN_F) {  
         passed=TRUE;  
         for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {  
            j=pattern->index[iptr];  
            if (mis_marker[j]==IS_IN_F) {  
           passed=TRUE;  
            }  
            else {  
           passed=FALSE;  
           break;  
            }  
         }  
         if (passed) mis_marker[i]=IS_IN_C;  
      }  
       }  
         
       /* swap to TRUE/FALSE in mis_marker */  
       #pragma omp parallel for private(i) schedule(static)  
       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F) ? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;;  
         
 }  
   
 void Paso_Coarsening_Local_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {  
     
    dim_t i,j;  
    /*double sum;*/  
    index_t iptr;  
    index_t num_colors;  
    index_t* colorOf;  
    register index_t color;  
    bool_t passed=FALSE;  
    dim_t n=pattern->numOutput;  
     
     
    colorOf=MEMALLOC(n,index_t);  
     
    if (pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_greedy: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
    Paso_Pattern_color(pattern,&num_colors,colorOf);  
     
    /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i)  
       if(mis_marker[i]==IS_AVAILABLE)  
      mis_marker[i]=IS_IN_F;  
         
       #pragma omp barrier  
       for (color=0;color<num_colors;++color) {  
      #pragma omp parallel for schedule(static) private(i,iptr,j)  
      for (i=0;i<n;++i) {  
         if (colorOf[i]==color) {    
            if (mis_marker[i]==IS_IN_F) {  
           for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {  
              j=pattern->index[iptr];  
              if (colorOf[j]<color)  
             mis_marker[j]=IS_IN_C;  
           }  
            }  
         }  
      }  
       }  
         
         
       #pragma omp barrier  
       for (color=0;color<num_colors;++color) {  
      #pragma omp parallel for schedule(static) private(i,iptr,j)  
      for (i=0;i<n;i++) {  
         if (colorOf[i]==color) {    
            if (mis_marker[i]==IS_IN_C) {  
           passed=TRUE;  
           for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {  
              j=pattern->index[iptr];  
              if (colorOf[j]<color && passed) {  
             if (mis_marker[j]==IS_IN_C) {  
                passed=TRUE;  
             }  
             else {  
                passed=FALSE;  
                /*break;*/  
             }  
              }  
           }  
           if (passed) mis_marker[i]=IS_IN_F;  
            }  
             
         }  
      }  
       }  
         
       /* swap to TRUE/FALSE in mis_marker */  
       #pragma omp parallel for private(i) schedule(static)  
       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F) ? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;  
         
       MEMFREE(colorOf);  
 }  
   
 /*For testing */  
 void Paso_Coarsening_Local_greedy_diag(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  
     
    dim_t i,j=0,k;  
    double *theta;  
    index_t iptr;  
    dim_t n=A->numRows;  
    double rsum,diag=0;  
    index_t *AvADJ;  
    theta=MEMALLOC(n,double);  
    AvADJ=MEMALLOC(n,index_t);  
     
     
     
     
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_coup: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
     
    #pragma omp parallel for private(i,iptr,j,rsum) schedule(static)  
    for (i=0;i<n;++i) {  
       rsum=0;  
       for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
      j=A->pattern->index[iptr];  
      if(j!=i) {  
         rsum+=ABS(A->val[iptr]);      
      }  
      else {  
         diag=ABS(A->val[iptr]);  
      }  
       }  
       theta[i]=diag/rsum;  
       if(theta[i]>threshold) {  
      mis_marker[i]=IS_IN_F;  
       }  
    }  
     
    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
       k=0;  
         
       for (i=0;i<n;++i) {  
      if(mis_marker[i]==IS_AVAILABLE) {  
         if(k==0) {  
            j=i;  
            k++;  
         }  
         if(theta[j]>theta[i]) {  
            j=i;  
         }  
      }  
       }  
       mis_marker[j]=IS_IN_C;  
         
       for (iptr=A->pattern->ptr[j];iptr<A->pattern->ptr[j+1]; ++iptr) {  
      k=A->pattern->index[iptr];  
      if(mis_marker[k]==IS_AVAILABLE) {  
         AvADJ[k]=1;  
      }  
      else {  
         AvADJ[k]=-1;  
      }  
       
       }  
         
       for (i=0;i<n;++i) {  
      if(AvADJ[i]) {  
         rsum=0;  
         for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
            k=A->pattern->index[iptr];  
            if(k!=i && mis_marker[k]!=IS_IN_C ) {  
           rsum+=ABS(A->val[iptr]);      
            }  
            if(j==i) {  
           diag=ABS(A->val[iptr]);  
            }  
         }  
         theta[i]=diag/rsum;  
         if(theta[i]>threshold) {  
            mis_marker[i]=IS_IN_F;  
         }  
      }  
       }  
         
         
    }  
     
    /* swap to TRUE/FALSE in mis_marker */  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F) ? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;  
     
    MEMFREE(AvADJ);  
    MEMFREE(theta);  
 }  
   
   
 void Paso_Coarsening_Local_YS_plus(Paso_SparseMatrix* A, index_t* mis_marker, double alpha, double taw, double delta) {  
     
    dim_t i,j;  
    /*double sum;*/  
    index_t iptr,*index,*where_p,*diagptr;  
    double *rsum;  
    double sum;  
    dim_t n=A->numRows;  
    Paso_SparseMatrix* A_alpha;  
    diagptr=MEMALLOC(n,index_t);  
    rsum=MEMALLOC(n,double);  
     
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_coup: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
    A_alpha=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern,1,1, FALSE);  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<A->len;++i) {  
       A_alpha->val[i]=A->val[i];  
    }  
     
     
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i)  
       if(mis_marker[i]==IS_AVAILABLE)  
      mis_marker[i]=IS_IN_C;  
         
       /*#pragma omp parallel for private(i,index,where_p) schedule(static)*/  
       for (i=0;i<n;++i) {  
      diagptr[i]=A->pattern->ptr[i];  
      index=&(A->pattern->index[A->pattern->ptr[i]]);  
      where_p=(index_t*)bsearch(&i,  
                    index,  
                    A->pattern->ptr[i + 1]-A->pattern->ptr[i],  
                    sizeof(index_t),  
                    Paso_comparIndex);  
                    if (where_p==NULL) {  
                       Esys_setError(VALUE_ERROR, "Paso_Coarsening_Local_coup: main diagonal element missing.");  
                    } else {  
                       diagptr[i]+=(index_t)(where_p-index);  
                    }  
       }  
         
         
       j=0;  
       for (i=0;i<n;++i) {  
      for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {  
         j=A_alpha->pattern->index[iptr];  
         if(i==j) {  
            A_alpha->val[iptr]=0;  
         }  
         else {  
            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) ) {  
           A_alpha->val[iptr]=0;  
            }  
         }  
           
      }  
       }  
         
         
       for (i=0;i<n;++i) {  
      rsum[i]=0;  
      for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {  
         rsum[i]+=A_alpha->val[iptr];    
      }  
       }  
         
       #pragma omp parallel for private(i,j) schedule(static)  
       for (i=0;i<n;++i) {  
      for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {  
         j=A_alpha->pattern->index[iptr];  
         if (i==j) {  
            A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr]+rsum[i];    
         }  
         else {  
            A_alpha->val[iptr]=A->val[iptr]-A_alpha->val[iptr];  
         }  
           
      }  
       }  
         
       /*This loop cannot be parallelized, as order matters here.*/  
       for (i=0;i<n;++i) {  
      if (mis_marker[i]==IS_IN_C) {  
         for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
            j=A->pattern->index[iptr];  
            if (j!=i && ABS(A->val[iptr])>=alpha*MIN(ABS(A->val[diagptr[i]]),ABS(A->val[diagptr[j]]))) {  
           mis_marker[j]=IS_IN_F;  
            }  
         }  
      }  
       }  
         
         
       /*This loop cannot be parallelized, as order matters here.*/  
       for (i=0;i<n;i++) {  
      if (mis_marker[i]==IS_IN_F) {  
         sum=0;  
         for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {  
            j=A_alpha->pattern->index[iptr];  
            if (mis_marker[j]==IS_IN_C) {  
           sum+=A_alpha->val[iptr];  
            }  
         }  
         if (ABS(sum)>=taw*ABS(A->val[diagptr[i]])) {  
            mis_marker[j]=IS_IN_FA;  
         }  
      }  
       }  
         
       /*This loop cannot be parallelized, as order matters here.*/  
       for (i=0;i<n;i++) {  
      if (mis_marker[i]!=IS_IN_C || mis_marker[i]!=IS_IN_FA) {  
         sum=0;  
         for (iptr=A_alpha->pattern->ptr[i];iptr<A_alpha->pattern->ptr[i+1]; ++iptr) {  
            j=A_alpha->pattern->index[iptr];  
            if (mis_marker[j]==IS_IN_C || mis_marker[j]==IS_IN_FA) {  
           sum+=A_alpha->val[iptr];  
            }  
         }  
         if (ABS(sum)>=delta*ABS(A->val[diagptr[i]])) {  
            mis_marker[j]=IS_IN_FB;  
         }  
         else {  
            mis_marker[j]=IS_IN_C;  
         }  
           
      }  
       }  
         
         
       /* swap to TRUE/FALSE in mis_marker */  
       #pragma omp parallel for private(i) schedule(static)  
       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_FA || mis_marker[i]==IS_IN_FB) ? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;  
         
       MEMFREE(diagptr);  
       MEMFREE(rsum);  
       Paso_SparseMatrix_free(A_alpha);  
 }  
   
   
 void Paso_Coarsening_Local_Standard(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  
 {  
    dim_t i,n,j,k;  
    index_t iptr,jptr;  
    /*index_t *index,*where_p;*/  
    double threshold,max_offdiagonal;  
    dim_t *lambda;   /*mesure of importance */  
    dim_t maxlambda=0;  
    index_t index_maxlambda=0;  
    double time0=0;  
    bool_t verbose=0;  
     
    Paso_Pattern *S=NULL;  
    Paso_Pattern *ST=NULL;  
    Paso_IndexList* index_list=NULL;  
     
    /*dim_t lk;*/  
     
    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) {  
      index_list[i].extension=NULL;  
      index_list[i].n=0;  
       }  
    }  
     
     
    n=A->numRows;  
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_RS: symmetric matrix pattern is not supported yet");  
       return;  
    }  
     
    time0=Esys_timer();  
    k=0;  
    /*S_i={j \in N_i; i strongly coupled to j}*/  
    #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)  
    for (i=0;i<n;++i) {  
       if(mis_marker[i]==IS_AVAILABLE) {  
      max_offdiagonal = DBL_MIN;  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         j=A->pattern->index[iptr];  
         if( j != i){  
            max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));  
         }  
      }  
      threshold = theta*max_offdiagonal;  
      for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
         j=A->pattern->index[iptr];  
         if(ABS(A->val[iptr])>threshold && i!=j) {  
            Paso_IndexList_insertIndex(&(index_list[i]),j);  
         }  
      }  
       }  
    }  
     
    S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);  
    ST=Paso_Coarsening_Local_getTranspose(S);  
     
    time0=Esys_timer()-time0;  
    if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);  
     
    lambda=TMPMEMALLOC(n,dim_t);  
     
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }  
     
    /*S_i={j \in N_i; i strongly coupled to j}*/  
     
    /*  
    #pragma omp parallel for private(i,iptr,lk) schedule(static)  
    for (i=0;i<n;++i) {  
       lk=0;  
       for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
      if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )  
         lk++;  
    }  
    #pragma omp critical  
    k+=lk;  
    if(k==0) {  
       mis_marker[i]=IS_IN_F;  
    }  
    }  
    */  
     
     
    k=0;  
    maxlambda=0;  
     
    time0=Esys_timer();  
     
    for (i=0;i<n;++i) {  
       if(mis_marker[i]==IS_AVAILABLE) {  
      lambda[i]=how_many(k,ST,FALSE);   /* if every row is available then k and i are the same.*/  
      /*lambda[i]=how_many(i,S,TRUE);*/  
      /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/  
      k++;  
      if(maxlambda<lambda[i]) {  
         maxlambda=lambda[i];  
         index_maxlambda=i;  
      }  
       }  
    }  
     
    k=0;  
    time0=Esys_timer()-time0;  
    if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);  
     
    time0=Esys_timer();  
     
    /*Paso_Coarsening_Local_getReport(n,mis_marker);*/  
     
    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {  
         
       if(index_maxlambda<0) {  
      break;  
       }  
         
       i=index_maxlambda;  
       if(mis_marker[i]==IS_AVAILABLE) {  
      mis_marker[index_maxlambda]=IS_IN_C;  
      lambda[index_maxlambda]=IS_NOT_AVAILABLE;  
      for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {  
         j=ST->index[iptr];  
         if(mis_marker[j]==IS_AVAILABLE) {  
            mis_marker[j]=IS_IN_F;  
            lambda[j]=IS_NOT_AVAILABLE;  
            for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {  
           k=S->index[jptr];  
           if(mis_marker[k]==IS_AVAILABLE) {  
              lambda[k]++;  
           }  
            }  
         }  
      }  
      for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {  
         j=S->index[iptr];  
         if(mis_marker[j]==IS_AVAILABLE) {  
            lambda[j]--;  
         }  
      }  
       }  
         
       /* Used when transpose of S is not available */  
       /*  
       for (i=0;i<n;++i) {  
      if(mis_marker[i]==IS_AVAILABLE) {  
         if (i==index_maxlambda) {  
            mis_marker[index_maxlambda]=IS_IN_C;  
            lambda[index_maxlambda]=-1;  
            for (j=0;j<n;++j) {  
           if(mis_marker[j]==IS_AVAILABLE) {  
              index=&(S->index[S->ptr[j]]);  
              where_p=(index_t*)bsearch(&i,  
              index,  
              S->ptr[j + 1]-S->ptr[j],  
              sizeof(index_t),  
              Paso_comparIndex);  
              if (where_p!=NULL) {  
             mis_marker[j]=IS_IN_F;  
             lambda[j]=-1;  
             for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {  
                k=S->index[iptr];  
                if(mis_marker[k]==IS_AVAILABLE) {  
                   lambda[k]++;  
       }  
       }  
       }  
         
       }  
       }  
         
       }  
       }  
       }  
       */  
       index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);  
    }  
     
    time0=Esys_timer()-time0;  
    if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);  
     
    /*Paso_Coarsening_Local_getReport(n,mis_marker);*/  
     
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;++i)  
       if(mis_marker[i]==IS_AVAILABLE) {  
      mis_marker[i]=IS_IN_F;  
       }  
         
       /*Paso_Coarsening_Local_getReport(n,mis_marker);*/  
         
       TMPMEMFREE(lambda);  
     
    /* clean up */  
    if (index_list!=NULL) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
    }  
     
    TMPMEMFREE(index_list);  
    Paso_Pattern_free(S);  
     
    /* swap to TRUE/FALSE in mis_marker */  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);  
     
 }  
   
 void Paso_Coarsening_Local_getReport(dim_t n,index_t* mis_marker) {  
     
    dim_t i,av,nav,f,c;  
     
    av=0;  
    nav=0;  
    f=0;  
    c=0;  
     
    for (i=0;i<n;i++) {  
       if(mis_marker[i]==IS_NOT_AVAILABLE) {  
      nav++;  
       }  
       else if (mis_marker[i]==IS_AVAILABLE) {  
      av++;  
       }  
       else if (mis_marker[i]==IS_IN_F) {  
      f++;  
       }  
       else if (mis_marker[i]==IS_IN_C) {  
      c++;  
       }  
    }  
     
    printf("Available=%d, Not Available = %d, In F = %d, In C = %d \n",av,nav,f,c);  
     
 }  
   
   
57  /*Used in Paso_Coarsening_Local_RS_MI*/  /*Used in Paso_Coarsening_Local_RS_MI*/
58    
59  /*Computes how many connections unknown i have in S.  /*Computes how many connections unknown i have in S.
# Line 1029  dim_t how_many(dim_t i,Paso_Pattern * S, Line 98  dim_t how_many(dim_t i,Paso_Pattern * S,
98        total=S->ptr[i+1]-S->ptr[i];        total=S->ptr[i+1]-S->ptr[i];
99        /*#pragma omp parallel for private(iptr) schedule(static)*/        /*#pragma omp parallel for private(iptr) schedule(static)*/
100        /*for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {        /*for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
101        if(S->index[iptr]!=i && mis_marker[S->index[iptr]]==IS_AVAILABLE)        if(S->index[iptr]!=i && marker_F[S->index[iptr]]==IS_AVAILABLE)
102       total++;       total++;
103     }*/     }*/
104                
# Line 1040  dim_t how_many(dim_t i,Paso_Pattern * S, Line 109  dim_t how_many(dim_t i,Paso_Pattern * S,
109     return total;     return total;
110  }  }
111    
112  dim_t arg_max(dim_t n, dim_t* lambda, dim_t mask) {  
113     dim_t i;  
    dim_t max=-1;  
    dim_t argmax=-1;  
     
    #ifdef _OPENMP  
    dim_t lmax=-1;  
    dim_t li=-1;  
    #pragma omp parallel private(i,lmax,li)  
    {  
       lmax=-1;  
       li=-1;  
       #pragma omp for schedule(static)  
       for (i=0;i<n;++i) {  
      if(lmax<lambda[i] && lambda[i]!=mask){  
         lmax=lambda[i];  
         li=i;  
      }  
       }  
       #pragma omp critical  
       {  
      if (max<=lmax) {  
         if(max==lmax && argmax>li)  
         {  
            argmax=li;  
         }  
         if (max < lmax)  
         {  
            max=lmax;  
            argmax=li;  
         }  
      }  
       }  
    }  
    #else  
    for (i=0;i<n;++i) {  
       if(max<lambda[i] && lambda[i]!=mask){  
      max=lambda[i];  
      argmax=i;  
       }  
    }  
    #endif  
     
    return argmax;  
 }  
114    
115    
116  Paso_Pattern* Paso_Coarsening_Local_getTranspose(Paso_Pattern* P){  Paso_Pattern* Paso_Coarsening_Local_getTranspose(Paso_Pattern* P){
117        
118     Paso_Pattern *outpattern=NULL;     Paso_Pattern *outpattern=NULL;
119        
    Paso_IndexList* index_list=NULL;  
120     dim_t C=P->numInput;     dim_t C=P->numInput;
121     dim_t F=P->numOutput-C;     dim_t F=P->numOutput-C;
122     dim_t n=C+F;     dim_t n=C+F;
123     dim_t i,j;     dim_t i,j;
124     index_t iptr;     index_t iptr;
125         Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(C);
    index_list=TMPMEMALLOC(C,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<C;++i) {  
      index_list[i].extension=NULL;  
      index_list[i].n=0;  
       }  
    }  
126        
127        
128     /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/     /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
129     for (i=0;i<n;++i) {     for (i=0;i<n;++i) {
130        for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {        for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {
131       j=P->index[iptr];       j=P->index[iptr];
132       Paso_IndexList_insertIndex(&(index_list[j]),i);       Paso_IndexListArray_insertIndex(index_list,i,j);
133        }        }
134     }     }
135        
136     outpattern=Paso_IndexList_createPattern(0, C,index_list,0,n,0);     outpattern=Paso_Pattern_fromIndexListArray(0, index_list,0,n,0);
137        
138     /* clean up */     /* clean up */
139     if (index_list!=NULL) {     Paso_IndexListArray_free(index_list);
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);  
    }  
    TMPMEMFREE(index_list);  
140        
141     return outpattern;     return outpattern;
142  }  }
# Line 1131  Paso_Pattern* Paso_Coarsening_Local_getT Line 144  Paso_Pattern* Paso_Coarsening_Local_getT
144    
145  /************** BLOCK COARSENENING *********************/  /************** BLOCK COARSENENING *********************/
146    
147  void Paso_Coarsening_Local_Standard_Block(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Coarsening_Local_Standard_Block(Paso_SparseMatrix* A, index_t* marker_F, double theta)
148  {  {
149     dim_t i,n,j,k;     const dim_t n=A->numRows;
150      
151       dim_t i,j,k;
152     index_t iptr,jptr;     index_t iptr,jptr;
153     /*index_t *index,*where_p;*/     /*index_t *index,*where_p;*/
154     double threshold,max_offdiagonal;     double threshold,max_offdiagonal;
# Line 1149  void Paso_Coarsening_Local_Standard_Bloc Line 164  void Paso_Coarsening_Local_Standard_Bloc
164        
165     Paso_Pattern *S=NULL;     Paso_Pattern *S=NULL;
166     Paso_Pattern *ST=NULL;     Paso_Pattern *ST=NULL;
167     Paso_IndexList* index_list=NULL;     Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(n);
     
    /*dim_t lk;*/  
     
    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) {  
      index_list[i].extension=NULL;  
      index_list[i].n=0;  
       }  
    }  
     
     
    n=A->numRows;  
    if (A->pattern->type & PATTERN_FORMAT_SYM) {  
       Esys_setError(TYPE_ERROR,"Paso_Coarsening_Local_RS: symmetric matrix pattern is not supported yet");  
       return;  
    }  
168        
169     time0=Esys_timer();     time0=Esys_timer();
170     k=0;     k=0;
171     /*Paso_Coarsening_Local_getReport(n,mis_marker);*/     /*Paso_Coarsening_Local_getReport(n,marker_F);*/
172     /*printf("Blocks %d %d\n",n_block,A->len);*/     /*printf("Blocks %d %d\n",n_block,A->len);*/
173        
174     /*S_i={j \in N_i; i strongly coupled to j}*/     /*S_i={j \in N_i; i strongly coupled to j}*/
175     #pragma omp parallel for private(i,bi,fnorm,iptr,max_offdiagonal,threshold,j) schedule(static)     #pragma omp parallel for private(i,bi,fnorm,iptr,max_offdiagonal,threshold,j) schedule(static)
176     for (i=0;i<n;++i) {     for (i=0;i<n;++i) {
177        if(mis_marker[i]==IS_AVAILABLE) {        if(marker_F[i]==IS_AVAILABLE) {
178       max_offdiagonal = DBL_MIN;       max_offdiagonal = DBL_MIN;
179       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) {
180          j=A->pattern->index[iptr];          j=A->pattern->index[iptr];
# Line 1202  void Paso_Coarsening_Local_Standard_Bloc Line 199  void Paso_Coarsening_Local_Standard_Bloc
199             }             }
200             fnorm=sqrt(fnorm);             fnorm=sqrt(fnorm);
201             if(fnorm>=threshold) {             if(fnorm>=threshold) {
202            Paso_IndexList_insertIndex(&(index_list[i]),j);            Paso_IndexListArray_insertIndex(index_list,i,j);
203             }             }
204          }          }
205       }       }
206        }        }
207     }     }
208        
209     S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);     S=Paso_Pattern_fromIndexListArray(0,index_list,0,A->pattern->numInput,0);
210     ST=Paso_Coarsening_Local_getTranspose(S);     ST=Paso_Coarsening_Local_getTranspose(S);
211        
212     /*printf("Patterns len %d %d\n",S->len,ST->len);*/     /*printf("Patterns len %d %d\n",S->len,ST->len);*/
# Line 1235  void Paso_Coarsening_Local_Standard_Bloc Line 232  void Paso_Coarsening_Local_Standard_Bloc
232     #pragma omp critical     #pragma omp critical
233     k+=lk;     k+=lk;
234     if(k==0) {     if(k==0) {
235        mis_marker[i]=IS_IN_F;        marker_F[i]=TRUE;
236     }     }
237     }     }
238     */     */
# Line 1247  void Paso_Coarsening_Local_Standard_Bloc Line 244  void Paso_Coarsening_Local_Standard_Bloc
244     time0=Esys_timer();     time0=Esys_timer();
245        
246     for (i=0;i<n;++i) {     for (i=0;i<n;++i) {
247        if(mis_marker[i]==IS_AVAILABLE) {        if(marker_F[i]==IS_AVAILABLE) {
248       lambda[i]=how_many(k,ST,FALSE);   /* if every row is available then k and i are the same.*/       lambda[i]=how_many(k,ST,FALSE);   /* if every row is available then k and i are the same.*/
249       /*lambda[i]=how_many(i,S,TRUE);*/       /*lambda[i]=how_many(i,S,TRUE);*/
250       /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/       /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/
# Line 1265  void Paso_Coarsening_Local_Standard_Bloc Line 262  void Paso_Coarsening_Local_Standard_Bloc
262        
263     time0=Esys_timer();     time0=Esys_timer();
264        
265     /*Paso_Coarsening_Local_getReport(n,mis_marker);*/     /*Paso_Coarsening_Local_getReport(n,marker_F);*/
266        
267     while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {     while (Paso_Util_isAny(n,marker_F,IS_AVAILABLE)) {
268                
269        if(index_maxlambda<0) {        if(index_maxlambda<0) {
270       break;       break;
271        }        }
272                
273        i=index_maxlambda;        i=index_maxlambda;
274        if(mis_marker[i]==IS_AVAILABLE) {        if(marker_F[i]==IS_AVAILABLE) {
275       mis_marker[index_maxlambda]=IS_IN_C;       marker_F[index_maxlambda]=FALSE;
276       lambda[index_maxlambda]=IS_NOT_AVAILABLE;       lambda[index_maxlambda]=IS_NOT_AVAILABLE;
277       for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {       for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
278          j=ST->index[iptr];          j=ST->index[iptr];
279          if(mis_marker[j]==IS_AVAILABLE) {          if(marker_F[j]==IS_AVAILABLE) {
280             mis_marker[j]=IS_IN_F;             marker_F[j]=TRUE;
281             lambda[j]=IS_NOT_AVAILABLE;             lambda[j]=IS_NOT_AVAILABLE;
282             for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {             for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
283            k=S->index[jptr];            k=S->index[jptr];
284            if(mis_marker[k]==IS_AVAILABLE) {            if(marker_F[k]==IS_AVAILABLE) {
285               lambda[k]++;               lambda[k]++;
286            }            }
287             }             }
# Line 1292  void Paso_Coarsening_Local_Standard_Bloc Line 289  void Paso_Coarsening_Local_Standard_Bloc
289       }       }
290       for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {       for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
291          j=S->index[iptr];          j=S->index[iptr];
292          if(mis_marker[j]==IS_AVAILABLE) {          if(marker_F[j]==IS_AVAILABLE) {
293             lambda[j]--;             lambda[j]--;
294          }          }
295       }       }
# Line 1301  void Paso_Coarsening_Local_Standard_Bloc Line 298  void Paso_Coarsening_Local_Standard_Bloc
298        /* Used when transpose of S is not available */        /* Used when transpose of S is not available */
299        /*        /*
300        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
301       if(mis_marker[i]==IS_AVAILABLE) {       if(marker_F[i]==IS_AVAILABLE) {
302          if (i==index_maxlambda) {          if (i==index_maxlambda) {
303             mis_marker[index_maxlambda]=IS_IN_C;             marker_F[index_maxlambda]=FALSE;
304             lambda[index_maxlambda]=-1;             lambda[index_maxlambda]=-1;
305             for (j=0;j<n;++j) {             for (j=0;j<n;++j) {
306            if(mis_marker[j]==IS_AVAILABLE) {            if(marker_F[j]==IS_AVAILABLE) {
307               index=&(S->index[S->ptr[j]]);               index=&(S->index[S->ptr[j]]);
308               where_p=(index_t*)bsearch(&i,               where_p=(index_t*)bsearch(&i,
309               index,               index,
# Line 1314  void Paso_Coarsening_Local_Standard_Bloc Line 311  void Paso_Coarsening_Local_Standard_Bloc
311               sizeof(index_t),               sizeof(index_t),
312               Paso_comparIndex);               Paso_comparIndex);
313               if (where_p!=NULL) {               if (where_p!=NULL) {
314              mis_marker[j]=IS_IN_F;              marker_F[j]=TRUE;
315              lambda[j]=-1;              lambda[j]=-1;
316              for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {              for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
317                 k=S->index[iptr];                 k=S->index[iptr];
318                 if(mis_marker[k]==IS_AVAILABLE) {                 if(marker_F[k]==IS_AVAILABLE) {
319                    lambda[k]++;                    lambda[k]++;
320        }        }
321        }        }
# Line 1337  void Paso_Coarsening_Local_Standard_Bloc Line 334  void Paso_Coarsening_Local_Standard_Bloc
334     time0=Esys_timer()-time0;     time0=Esys_timer()-time0;
335     if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);     if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
336        
337     /*Paso_Coarsening_Local_getReport(n,mis_marker);*/     /*Paso_Coarsening_Local_getReport(n,marker_F);*/
338        
339     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
340     for (i=0;i<n;++i)     for (i=0;i<n;++i)
341        if(mis_marker[i]==IS_AVAILABLE) {        if(marker_F[i]==IS_AVAILABLE) {
342       mis_marker[i]=IS_IN_F;       marker_F[i]=TRUE;
343        }        }
344                
345        /*Paso_Coarsening_Local_getReport(n,mis_marker);*/        /*Paso_Coarsening_Local_getReport(n,marker_F);*/
346                
347        TMPMEMFREE(lambda);        TMPMEMFREE(lambda);
348        
349     /* clean up */     /* clean up */
350     if (index_list!=NULL) {     Paso_IndexListArray_free(index_list);
       #pragma omp parallel for private(i) schedule(static)  
       for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
    }  
     
    TMPMEMFREE(index_list);  
351     Paso_Pattern_free(S);     Paso_Pattern_free(S);
352        
353     /* swap to TRUE/FALSE in mis_marker */     /* swap to TRUE/FALSE in marker_F */
354     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
355     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F)? PASO_COARSENING_IN_F : PASO_COARSENING_IN_C;     for (i=0;i<n;i++) marker_F[i]=(marker_F[i]==TRUE)? TRUE : FALSE;
356        
357  }  }
358    
# Line 1368  void Paso_Coarsening_Local_Standard_Bloc Line 360  void Paso_Coarsening_Local_Standard_Bloc
360    
361  #undef IS_AVAILABLE  #undef IS_AVAILABLE
362  #undef IS_NOT_AVAILABLE  #undef IS_NOT_AVAILABLE
363  #undef IS_IN_F  #undef TRUE
364  #undef IS_IN_C  #undef FALSE
365    
366  #undef IS_UNDECIDED  #undef IS_UNDECIDED
367  #undef IS_STRONG  #undef IS_STRONG
368  #undef IS_WEAK  #undef IS_WEAK
369    
370  #undef IS_IN_FB  #undef TRUEB
371  #undef IS_IN_FB  #undef TRUEB
372    

Legend:
Removed from v.3282  
changed lines
  Added in v.3283

  ViewVC Help
Powered by ViewVC 1.1.26