/[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 2827 by artak, Thu Dec 10 02:09:43 2009 UTC revision 2828 by artak, Tue Dec 22 01:24:40 2009 UTC
# Line 313  void Paso_Pattern_greedy(Paso_Pattern* p Line 313  void Paso_Pattern_greedy(Paso_Pattern* p
313            
314  }  }
315    
   
316  void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {  void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
317    
318    dim_t i,j;    dim_t i,j;
# Line 648  void Paso_Pattern_Standard(Paso_SparseMa Line 647  void Paso_Pattern_Standard(Paso_SparseMa
647    Paso_Pattern *ST=NULL;    Paso_Pattern *ST=NULL;
648    Paso_IndexList* index_list=NULL;    Paso_IndexList* index_list=NULL;
649        
   index_t* counter;  
650    /*dim_t lk;*/    /*dim_t lk;*/
651    
652    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);    index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
# Line 685  void Paso_Pattern_Standard(Paso_SparseMa Line 683  void Paso_Pattern_Standard(Paso_SparseMa
683          threshold = theta*max_offdiagonal;          threshold = theta*max_offdiagonal;
684          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) {
685              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
686              if((-A->val[iptr])>=threshold) {              if(ABS(A->val[iptr])>=threshold) {
687                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
688              }              }
689          }          }
690        }        }
691      }      }
692        
     
693    S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);    S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
694    ST=Paso_Pattern_getTranspose(S);    ST=Paso_Pattern_getTranspose(S);
695        
# Line 703  void Paso_Pattern_Standard(Paso_SparseMa Line 700  void Paso_Pattern_Standard(Paso_SparseMa
700    lambda=TMPMEMALLOC(n,dim_t);    lambda=TMPMEMALLOC(n,dim_t);
701        
702    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
703    for (i=0;i<n;++i) { lambda[i]=-1; }    for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }
704        
705    /*S_i={j \in N_i; i strongly coupled to j}*/    /*S_i={j \in N_i; i strongly coupled to j}*/
706    
# Line 731  void Paso_Pattern_Standard(Paso_SparseMa Line 728  void Paso_Pattern_Standard(Paso_SparseMa
728        
729      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
730        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
731          lambda[i]=how_many(k,ST,FALSE);          lambda[i]=how_many(k,ST,FALSE);   /* if every row is available then k and i are the same.*/
732          /*lambda[i]=how_many(i,S,TRUE);*/          /*lambda[i]=how_many(i,S,TRUE);*/
733          /*printf("lambda[%d]=%d, k=%d ",i,lambda[i],k);*/          /*printf("lambda[%d]=%d, k=%d ",i,lambda[i],k);*/
734          k++;          k++;
# Line 749  void Paso_Pattern_Standard(Paso_SparseMa Line 746  void Paso_Pattern_Standard(Paso_SparseMa
746    time0=Paso_timer();    time0=Paso_timer();
747        
748    /*Paso_Pattern_getReport(n,mis_marker);*/    /*Paso_Pattern_getReport(n,mis_marker);*/
     counter=TMPMEMALLOC(n,index_t);  
   
     #pragma omp parallel for private(i) schedule(static)  
     for (i = 0; i < n; ++i) counter[i]=(mis_marker[i]==IS_AVAILABLE);  
       
     Paso_Util_cumsum(n,counter);  
                 
749        
750    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
751            
# Line 766  void Paso_Pattern_Standard(Paso_SparseMa Line 756  void Paso_Pattern_Standard(Paso_SparseMa
756      i=index_maxlambda;      i=index_maxlambda;
757      if(mis_marker[i]==IS_AVAILABLE) {      if(mis_marker[i]==IS_AVAILABLE) {
758          mis_marker[index_maxlambda]=IS_IN_C;          mis_marker[index_maxlambda]=IS_IN_C;
759          lambda[index_maxlambda]=-1;          lambda[index_maxlambda]=IS_NOT_AVAILABLE;
760          for (iptr=ST->ptr[counter[i]];iptr<ST->ptr[counter[i]+1]; ++iptr) {          for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
761              j=ST->index[iptr];              j=ST->index[iptr];
762              if(mis_marker[j]==IS_AVAILABLE) {              if(mis_marker[j]==IS_AVAILABLE) {
763                  mis_marker[j]=IS_IN_F;                  mis_marker[j]=IS_IN_F;
764                  lambda[j]=-1;                  lambda[j]=IS_NOT_AVAILABLE;
765                  for (jptr=S->ptr[counter[j]];jptr<S->ptr[counter[j]+1]; ++jptr) {                  for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
766                          k=S->index[jptr];                          k=S->index[jptr];
767                          if(mis_marker[k]==IS_AVAILABLE) {                          if(mis_marker[k]==IS_AVAILABLE) {
768                             lambda[k]++;                             lambda[k]++;
# Line 815  void Paso_Pattern_Standard(Paso_SparseMa Line 805  void Paso_Pattern_Standard(Paso_SparseMa
805          }          }
806      }      }
807     */     */
808     index_maxlambda=arg_max(n,lambda, -1);     index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);
809    }    }
810        
811    time0=Paso_timer()-time0;    time0=Paso_timer()-time0;
# Line 825  void Paso_Pattern_Standard(Paso_SparseMa Line 815  void Paso_Pattern_Standard(Paso_SparseMa
815        
816    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
817    for (i=0;i<n;++i)    for (i=0;i<n;++i)
818        if(mis_marker[i]==IS_AVAILABLE || mis_marker[i]==IS_NOT_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
819          mis_marker[i]=IS_IN_F;          mis_marker[i]=IS_IN_F;
820         }         }
821                
# Line 848  void Paso_Pattern_Standard(Paso_SparseMa Line 838  void Paso_Pattern_Standard(Paso_SparseMa
838            
839  }  }
840    
 /* Not ready yet for usage, needs intermidiate matrix creation which can be costly. */  
 void Paso_Pattern_Aggressive(Paso_SparseMatrix* A, index_t* mis_marker, double theta) {  
     dim_t i;  
       
     Paso_Pattern_Standard(A, mis_marker, theta);  
       
     #pragma omp parallel for private(i) schedule(static)  
     for (i=0;i<A->numRows;i++) {  
          if(!mis_marker[i]) {  
             mis_marker[i]=IS_AVAILABLE;  
          }  
          else {  
             mis_marker[i]=IS_NOT_AVAILABLE;  
          }  
     }  
       
     Paso_Pattern_Standard(A, mis_marker, theta);  
 }  
   
841  void Paso_Pattern_getReport(dim_t n,index_t* mis_marker) {  void Paso_Pattern_getReport(dim_t n,index_t* mis_marker) {
842            
843      dim_t i,av,nav,f,c;      dim_t i,av,nav,f,c;
# Line 946  dim_t how_many(dim_t i,Paso_Pattern * S, Line 917  dim_t how_many(dim_t i,Paso_Pattern * S,
917                    
918      }      }
919            
920      if (total==0) total=-1;      if (total==0) total=IS_NOT_AVAILABLE;
921            
922  return total;  return total;
923  }  }

Legend:
Removed from v.2827  
changed lines
  Added in v.2828

  ViewVC Help
Powered by ViewVC 1.1.26