/[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 2830 by artak, Tue Dec 22 01:24:40 2009 UTC revision 2831 by artak, Thu Jan 7 05:08:19 2010 UTC
# Line 47  Line 47 
47    
48    
49  #define IS_IN_FA -5  /* test */  #define IS_IN_FA -5  /* test */
50  #define IS_IN_FB -6  /* test */  #define IS_IN_FB -6  /* test */
51    
52    
53    void Paso_Pattern_Read(char *fileName,dim_t n,index_t* mis_marker) {
54        
55        dim_t i;
56        int scan_ret;
57        
58        FILE *fileHandle_p = NULL;
59        
60        fileHandle_p = fopen( fileName, "r" );
61        if( fileHandle_p == NULL )
62        {
63            Paso_setError(IO_ERROR, "Paso_Pattern_Read: Cannot open file for reading.");
64            return;
65        }
66        
67        for (i=0;i<n;++i) {    
68            scan_ret=fscanf( fileHandle_p, "%d\n", &mis_marker[i]);
69            if (scan_ret!=1) {
70                Paso_setError(IO_ERROR, "Paso_Pattern_Read: Cannot read line from file.");
71            return;
72            }
73        }
74        
75        fclose(fileHandle_p);
76        
77        /* swap to TRUE/FALSE in mis_marker */
78         #pragma omp parallel for private(i) schedule(static)
79         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=1);
80    }
81    
82    
83    
84  void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
85    
# Line 157  void Paso_Pattern_RS(Paso_SparseMatrix* Line 189  void Paso_Pattern_RS(Paso_SparseMatrix*
189          max_offdiagonal = DBL_MIN;          max_offdiagonal = DBL_MIN;
190          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) {
191              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
192                    max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);                    max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
193              }              }
194          }          }
195                    
196          threshold = theta*max_offdiagonal;          threshold = theta*max_offdiagonal;
197          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) {
198              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
199              if((-A->val[iptr])>=threshold) {              if(ABS(A->val[iptr])>=threshold && i!=j) {
200                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
201                  /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/                  /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
202              }              }
# Line 675  void Paso_Pattern_Standard(Paso_SparseMa Line 707  void Paso_Pattern_Standard(Paso_SparseMa
707          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) {
708              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
709              if( j != i){              if( j != i){
                 if(A->val[iptr]<0) {  
710                    max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));                    max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
                 }  
711              }              }
712          }          }
713          threshold = theta*max_offdiagonal;          threshold = theta*max_offdiagonal;
714          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) {
715              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
716              if(ABS(A->val[iptr])>=threshold) {              if(ABS(A->val[iptr])>=threshold && i!=j) {
717                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
718              }              }
719          }          }
# Line 692  void Paso_Pattern_Standard(Paso_SparseMa Line 722  void Paso_Pattern_Standard(Paso_SparseMa
722        
723    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);
724    ST=Paso_Pattern_getTranspose(S);    ST=Paso_Pattern_getTranspose(S);
     
725            
726    time0=Paso_timer()-time0;    time0=Paso_timer()-time0;
727    if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);    if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
# Line 730  void Paso_Pattern_Standard(Paso_SparseMa Line 759  void Paso_Pattern_Standard(Paso_SparseMa
759        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
760          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.*/
761          /*lambda[i]=how_many(i,S,TRUE);*/          /*lambda[i]=how_many(i,S,TRUE);*/
762          /*printf("lambda[%d]=%d, k=%d ",i,lambda[i],k);*/          /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/
763          k++;          k++;
764          if(maxlambda<lambda[i]) {          if(maxlambda<lambda[i]) {
765              maxlambda=lambda[i];              maxlambda=lambda[i];
# Line 770  void Paso_Pattern_Standard(Paso_SparseMa Line 799  void Paso_Pattern_Standard(Paso_SparseMa
799                  }                  }
800              }              }
801          }          }
802            for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
803                j=S->index[iptr];
804                if(mis_marker[j]==IS_AVAILABLE) {
805                               lambda[j]--;
806                }
807            }
808                
809      }      }
810            
811     /* Used when transpose of S is not available */     /* Used when transpose of S is not available */
Legend:
Removed from v.2830  
changed lines
  Added in v.2831
  ViewVC Help
Powered by ViewVC 1.1.26