/[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 2764 by artak, Thu Nov 19 23:48:05 2009 UTC revision 2765 by artak, Fri Nov 20 01:49:19 2009 UTC
# Line 639  void Paso_Pattern_YS_plus(Paso_SparseMat Line 639  void Paso_Pattern_YS_plus(Paso_SparseMat
639  void Paso_Pattern_RS_MI(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_RS_MI(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
640  {  {
641    dim_t i,n,j,k;    dim_t i,n,j,k;
642    index_t iptr,*index,*where_p;    index_t iptr,jptr;
643      /*index_t *index,*where_p;*/
644    double threshold,max_offdiagonal;    double threshold,max_offdiagonal;
645    dim_t *lambda;   /*mesure of importance */    dim_t *lambda;   /*mesure of importance */
646    /*bool_t breakloop=FALSE;*/    /*bool_t breakloop=FALSE;*/
647    dim_t maxlambda=0;    dim_t maxlambda=0;
648    index_t index_maxlambda=0;    index_t index_maxlambda=0;
649      double time0=0;
650      bool_t verbose=0;
651        
652    Paso_Pattern *S=NULL;    Paso_Pattern *S=NULL;
653      Paso_Pattern *ST=NULL;
654    Paso_IndexList* index_list=NULL;    Paso_IndexList* index_list=NULL;
655    
656   index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);   index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
# Line 665  void Paso_Pattern_RS_MI(Paso_SparseMatri Line 669  void Paso_Pattern_RS_MI(Paso_SparseMatri
669      return;      return;
670    }    }
671        
672        time0=Paso_timer();
673      
674      /*S_i={j \in N_i; i strongly coupled to j}*/      /*S_i={j \in N_i; i strongly coupled to j}*/
675      /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/      /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/
676      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
# Line 687  void Paso_Pattern_RS_MI(Paso_SparseMatri Line 693  void Paso_Pattern_RS_MI(Paso_SparseMatri
693          }          }
694         }         }
695        }        }
       
696        
697      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);
698      ST=Paso_Pattern_getTranspose(S);
699      
700      time0=Paso_timer()-time0;
701      if (verbose) fprintf(stderr,"timing: RS filtering and pattern creation: %e\n",time0);
702        
703    lambda=TMPMEMALLOC(n,dim_t);    lambda=TMPMEMALLOC(n,dim_t);
704        
705    for (i=0;i<n;++i) {    #pragma omp parallel for private(i) schedule(static)
706       lambda[i]=-1;    for (i=0;i<n;++i) { lambda[i]=-1; }
    }  
707        
708    /*S_i={j \in N_i; i strongly coupled to j}*/    /*S_i={j \in N_i; i strongly coupled to j}*/
709    
710    k=0;    k=0;
711    maxlambda=0;    maxlambda=0;
712        
713      
714      time0=Paso_timer();
715      
716      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
717        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
718          lambda[i]=how_many(i,S,TRUE);          lambda[i]=how_many(i,ST,FALSE);
719            /*lambda[i]=how_many(i,S,TRUE);*/
720          /*printf("lambda[%d]=%d, ",i,lambda[i]);*/          /*printf("lambda[%d]=%d, ",i,lambda[i]);*/
721          if(maxlambda<lambda[i]) {          if(maxlambda<lambda[i]) {
722              maxlambda=lambda[i];              maxlambda=lambda[i];
# Line 712  void Paso_Pattern_RS_MI(Paso_SparseMatri Line 724  void Paso_Pattern_RS_MI(Paso_SparseMatri
724          }          }
725        }        }
726      }      }
727    
728      time0=Paso_timer()-time0;
729      if (verbose) fprintf(stderr,"timing: Lambdas computations at the begining: %e\n",time0);
730      
731      time0=Paso_timer();
732        
733    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {    while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
734            
# Line 720  void Paso_Pattern_RS_MI(Paso_SparseMatri Line 737  void Paso_Pattern_RS_MI(Paso_SparseMatri
737          break;          break;
738      }      }
739    
740    
741    
742        for (i=0;i<n;++i) {
743            if(mis_marker[i]==IS_AVAILABLE) {
744                if (i==index_maxlambda) {
745                    mis_marker[index_maxlambda]=IS_IN_C;
746                    lambda[index_maxlambda]=-1;
747                    for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
748                        j=ST->index[iptr];
749                        if(mis_marker[j]==IS_AVAILABLE) {
750                            mis_marker[j]=IS_IN_F;
751                            lambda[j]=-1;
752                            for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
753                                    k=S->index[jptr];
754                                    if(mis_marker[k]==IS_AVAILABLE) {
755                                       lambda[k]++;
756                                    }
757                            }
758                        }
759                    }
760                }
761            }
762        }
763    
764    
765        
766       /* Used when transpose of S is not available */
767       /*
768      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
769          if(mis_marker[i]==IS_AVAILABLE) {          if(mis_marker[i]==IS_AVAILABLE) {
770              if (i==index_maxlambda) {              if (i==index_maxlambda) {
# Line 750  void Paso_Pattern_RS_MI(Paso_SparseMatri Line 795  void Paso_Pattern_RS_MI(Paso_SparseMatri
795              }              }
796          }          }
797      }      }
798           */
799    }    }
800        
801      time0=Paso_timer()-time0;
802      if (verbose) fprintf(stderr,"timing: Loop : %e\n",time0);
803      
804      #pragma omp parallel for private(i) schedule(static)
805    for (i=0;i<n;++i)    for (i=0;i<n;++i)
806        if(mis_marker[i]==IS_AVAILABLE)        if(mis_marker[i]==IS_AVAILABLE)
807          mis_marker[i]=IS_IN_F;          mis_marker[i]=IS_IN_F;
808    
     /*update lambdas*/  
     /*for (i=0;i<n;++i) {  
       if(mis_marker[i]==IS_AVAILABLE) {  
         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);  
         if(maxlambda<lambda[i]) {  
             maxlambda=lambda[i];  
             index_maxlambda=i;  
         }  
       }  
       if(lambda[i]==0) {  
         breakloop=TRUE;  
         break;  
       }  
     }  
     if(breakloop) {  
         break;  
     }  
       
     for (i=0;i<n;++i) {  
         if(mis_marker[i]==IS_AVAILABLE) {  
             mis_marker[index_maxlambda]=IS_IN_C;  
         }  
           
         for (j=0;j<n;++j) {  
             if(S_T_[i][j]=IS_STRONG && mis_marker[i]==IS_AVAILABLE) {  
                 mis_marker[j]==IS__IN_F;  
             }  
         }  
     }  
       
     }  
     */  
809    
810      TMPMEMFREE(lambda);      TMPMEMFREE(lambda);
811            
# Line 854  dim_t arg_max(dim_t n, dim_t* lambda, di Line 871  dim_t arg_max(dim_t n, dim_t* lambda, di
871  }  }
872    
873    
874    Paso_Pattern* Paso_Pattern_getTranspose(Paso_Pattern* P){
875    
876      Paso_Pattern *outpattern=NULL;
877      
878      Paso_IndexList* index_list=NULL;
879      dim_t C=P->numInput;
880      dim_t F=P->numOutput-C;
881      dim_t n=C+F;
882      dim_t i,j;
883      index_t iptr;
884    
885      index_list=TMPMEMALLOC(C,Paso_IndexList);
886       if (! Paso_checkPtr(index_list)) {
887            #pragma omp parallel for private(i) schedule(static)
888            for(i=0;i<C;++i) {
889                 index_list[i].extension=NULL;
890                 index_list[i].n=0;
891            }
892        }
893      
894    
895        /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
896        for (i=0;i<n;++i) {
897              for (iptr=P->ptr[i];iptr<P->ptr[i+1]; ++iptr) {
898                 j=P->index[iptr];
899                 Paso_IndexList_insertIndex(&(index_list[j]),i);
900            }
901        }
902      
903        outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0);
904        
905         /* clean up */
906       if (index_list!=NULL) {
907            #pragma omp parallel for private(i) schedule(static)
908            for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
909         }
910        TMPMEMFREE(index_list);
911    
912        return outpattern;
913    }
914    
915    
916  /*dim_t how_many(dim_t n,dim_t* S_i, int value1, dim_t* addedSet, int value2) {  /*dim_t how_many(dim_t n,dim_t* S_i, int value1, dim_t* addedSet, int value2) {
917      dim_t j;      dim_t j;
918      dim_t total;      dim_t total;

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

  ViewVC Help
Powered by ViewVC 1.1.26