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

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

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

revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3283 by gross, Mon Oct 18 22:39:28 2010 UTC
# Line 37  Paso_Pattern* Paso_Pattern_alloc(int typ Line 37  Paso_Pattern* Paso_Pattern_alloc(int typ
37    dim_t i;    dim_t i;
38    Esys_resetError();    Esys_resetError();
39    
   if (type & PATTERN_FORMAT_SYM) {  
     Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: symmetric matrix pattern is not supported yet");  
     return NULL;  
   }  
40    if (ptr!=NULL && index != NULL) {    if (ptr!=NULL && index != NULL) {
41      #pragma omp parallel private(loc_min_index,loc_max_index,i)      #pragma omp parallel private(loc_min_index,loc_max_index,i)
42       {       {
# Line 165  Paso_Pattern* Paso_Pattern_multiply(int Line 161  Paso_Pattern* Paso_Pattern_multiply(int
161    Paso_Pattern*out=NULL;    Paso_Pattern*out=NULL;
162    index_t iptrA,iptrB;    index_t iptrA,iptrB;
163    dim_t i,j,k;    dim_t i,j,k;
164    Paso_IndexList* index_list=NULL;    Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
   
   index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);  
   if (! Esys_checkPtr(index_list)) {  
         #pragma omp parallel for private(i) schedule(static)  
         for(i=0;i<A->numOutput;++i) {  
              index_list[i].extension=NULL;  
              index_list[i].n=0;  
         }  
    }  
165        
166    #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)    #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
167    for(i = 0; i < A->numOutput; i++) {    for(i = 0; i < A->numOutput; i++) {
# Line 182  Paso_Pattern* Paso_Pattern_multiply(int Line 169  Paso_Pattern* Paso_Pattern_multiply(int
169        j = A->index[iptrA];        j = A->index[iptrA];
170        for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {        for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {
171          k = B->index[iptrB];          k = B->index[iptrB];
172          Paso_IndexList_insertIndex(&(index_list[i]),k);      Paso_IndexListArray_insertIndex(index_list,i,k);
173       }       }
174      }      }
175    }    }
176            
177    out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,B->numInput,0);    out=Paso_Pattern_fromIndexListArray(0, index_list,0,B->numInput,0);
178    
179   /* clean up */   /* clean up */
180     if (index_list!=NULL) {   Paso_IndexListArray_free(index_list);
        #pragma omp parallel for private(i) schedule(static)  
         for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
      }  
   TMPMEMFREE(index_list);  
181        
182  return out;  return out;
183  }  }
# Line 212  Paso_Pattern* Paso_Pattern_binop(int typ Line 195  Paso_Pattern* Paso_Pattern_binop(int typ
195    index_t iptrA,iptrB;    index_t iptrA,iptrB;
196    dim_t i,j,k;    dim_t i,j,k;
197    
198    Paso_IndexList* index_list=NULL;    Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
   
  index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);  
    if (! Esys_checkPtr(index_list)) {  
         #pragma omp parallel for private(i) schedule(static)  
         for(i=0;i<A->numOutput;++i) {  
              index_list[i].extension=NULL;  
              index_list[i].n=0;  
         }  
     }  
199        
200    #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)    #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
201    for(i = 0; i < B->numOutput; i++){    for(i = 0; i < B->numOutput; i++){
# Line 232  Paso_Pattern* Paso_Pattern_binop(int typ Line 206  Paso_Pattern* Paso_Pattern_binop(int typ
206          j = A->index[iptrA];          j = A->index[iptrA];
207          k = B->index[iptrB];          k = B->index[iptrB];
208          if (j<k) {          if (j<k) {
209             Paso_IndexList_insertIndex(&(index_list[i]),j);         Paso_IndexListArray_insertIndex(index_list,i,j);
210             iptrA++;             iptrA++;
211          } else if (j>k) {          } else if (j>k) {
212              Paso_IndexList_insertIndex(&(index_list[i]),k);         Paso_IndexListArray_insertIndex(index_list,i,k);
213              iptrB++;              iptrB++;
214          } else if (j==k) {          } else if (j==k) {
215              Paso_IndexList_insertIndex(&(index_list[i]),j);         Paso_IndexListArray_insertIndex(index_list,i,j);
216              iptrB++;              iptrB++;
217              iptrA++;              iptrA++;
218          }          }
219      }      }
220      while(iptrA < A->ptr[i+1]) {      while(iptrA < A->ptr[i+1]) {
221          j = A->index[iptrA];          j = A->index[iptrA];
222          Paso_IndexList_insertIndex(&(index_list[i]),j);      Paso_IndexListArray_insertIndex(index_list,i,j);
223          iptrA++;          iptrA++;
224      }      }
225      while(iptrB < B->ptr[i+1]) {      while(iptrB < B->ptr[i+1]) {
226          k = B->index[iptrB];          k = B->index[iptrB];
227          Paso_IndexList_insertIndex(&(index_list[i]),k);      Paso_IndexListArray_insertIndex(index_list,i,k);
228          iptrB++;          iptrB++;
229      }      }
230    }    }
231    
232    out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,A->numInput,0);    out=Paso_Pattern_fromIndexListArray(0, index_list,0,A->numInput,0);
233    
234    
235   /* clean up */   /* clean up */
236     if (index_list!=NULL) {   Paso_IndexListArray_free(index_list);
         #pragma omp parallel for private(i) schedule(static)  
         for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);  
      }  
   TMPMEMFREE(index_list);  
237    
238    return out;    return out;
239  }  }
240    
 /* inserts row index row into the Paso_IndexList in if it does not exist */  
   
 void Paso_IndexList_insertIndex(Paso_IndexList* in, index_t index) {  
   dim_t i;  
   /* is index in in? */  
   for (i=0;i<in->n;i++) {  
     if (in->index[i]==index)  return;  
   }  
   /* index could not be found */  
   if (in->n==INDEXLIST_LENGTH) {  
      /* if in->index is full check the extension */  
      if (in->extension==NULL) {  
         in->extension=TMPMEMALLOC(1,Paso_IndexList);  
         if (Esys_checkPtr(in->extension)) return;  
         in->extension->n=0;  
         in->extension->extension=NULL;  
      }  
      Paso_IndexList_insertIndex(in->extension,index);  
   } else {  
      /* insert index into in->index*/  
      in->index[in->n]=index;  
      in->n++;  
   }  
 }  
   
 /* counts the number of row indices in the Paso_IndexList in */  
   
 dim_t Paso_IndexList_count(Paso_IndexList* in, index_t range_min,index_t range_max) {  
   dim_t i;  
   dim_t out=0;  
   register index_t itmp;  
   if (in==NULL) {  
      return 0;  
   } else {  
     for (i=0;i<in->n;i++) {  
           itmp=in->index[i];  
           if ((itmp>=range_min) && (range_max>itmp)) ++out;  
     }  
      return out+Paso_IndexList_count(in->extension, range_min,range_max);  
   }  
 }  
   
 /* count the number of row indices in the Paso_IndexList in */  
   
 void Paso_IndexList_toArray(Paso_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {  
   dim_t i, ptr;  
   register index_t itmp;  
   if (in!=NULL) {  
     ptr=0;  
     for (i=0;i<in->n;i++) {  
           itmp=in->index[i];  
           if ((itmp>=range_min) && (range_max>itmp)) {  
              array[ptr]=itmp+index_offset;  
              ptr++;  
           }  
   
     }  
     Paso_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);  
   }  
 }  
   
 /* deallocates the Paso_IndexList in by recursive calls */  
   
 void Paso_IndexList_free(Paso_IndexList* in) {  
   if (in!=NULL) {  
     Paso_IndexList_free(in->extension);  
     TMPMEMFREE(in);  
   }  
 }  
   
241  /* creates a Paso_pattern from a range of indices */  /* creates a Paso_pattern from a range of indices */
242  Paso_Pattern* Paso_IndexList_createPattern(dim_t n0, dim_t n,Paso_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)  Paso_Pattern* Paso_Pattern_fromIndexListArray(dim_t n0, Paso_IndexListArray* index_list_array,index_t range_min,index_t range_max,index_t index_offset)
243  {  {
244       const dim_t n=index_list_array->n;
245       Paso_IndexList* index_list = index_list_array->index_list;
246     dim_t *ptr=NULL;     dim_t *ptr=NULL;
247     register dim_t s,i,itmp;     register dim_t s,i,itmp;
248     index_t *index=NULL;     index_t *index=NULL;
# Line 348  Paso_Pattern* Paso_IndexList_createPatte Line 250  Paso_Pattern* Paso_IndexList_createPatte
250    
251     ptr=MEMALLOC(n+1-n0,index_t);     ptr=MEMALLOC(n+1-n0,index_t);
252     if (! Esys_checkPtr(ptr) ) {     if (! Esys_checkPtr(ptr) ) {
253         /* get the number of connections per row */         /* get the number of connections per row */
254         #pragma omp parallel for private(i) schedule(static)         #pragma omp parallel for private(i) schedule(static)
255         for(i=n0;i<n;++i) {         for(i=n0;i<n;++i) {
256                ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);                ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
# Line 383  index_t* Paso_Pattern_borrowMainDiagonal Line 285  index_t* Paso_Pattern_borrowMainDiagonal
285  {  {
286      const dim_t n=A->numOutput;      const dim_t n=A->numOutput;
287      int fail=0;      int fail=0;
288      index_t i,iptr,iptr_main;      index_t *index,*where_p, i;
289            
290       if (A->main_iptr == NULL) {       if (A->main_iptr == NULL) {
291           A->main_iptr=MEMALLOC(n,index_t);           A->main_iptr=MEMALLOC(n,index_t);
# Line 391  index_t* Paso_Pattern_borrowMainDiagonal Line 293  index_t* Paso_Pattern_borrowMainDiagonal
293           #pragma omp parallel           #pragma omp parallel
294               {               {
295                   /* identify the main diagonals */                   /* identify the main diagonals */
296                   #pragma omp for schedule(static) private(i,iptr,iptr_main)                   #pragma omp for schedule(static) private(i, index, where_p)
297                   for (i = 0; i < n; ++i) {                   for (i = 0; i < n; ++i) {
298                          iptr_main=A->ptr[0]-1;              index=&(A->index[A->ptr[i]]);
299                          for (iptr=A->ptr[i];iptr<A->ptr[i+1]; iptr++) {              where_p=bsearch(&i,
300                                if (A->index[iptr]==i) {                      index,
301                                     iptr_main=iptr;                      (size_t) (A->ptr[i + 1]-A->ptr[i]),
302                                     break;                           sizeof(index_t),
303                                }                           Paso_comparIndex);
304                          }                            
305                          A->main_iptr[i]=iptr_main;              if (where_p==NULL) {
306                          if (iptr_main==A->ptr[0]-1) fail=1;                  fail=1;
307                    }              } else {
308                   A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);
309                }
310                     }
311            
312               }               }
313           if (fail > 0) {           if (fail > 0) {
# Line 414  index_t* Paso_Pattern_borrowMainDiagonal Line 319  index_t* Paso_Pattern_borrowMainDiagonal
319       }       }
320       return A->main_iptr;       return A->main_iptr;
321  }  }
322                  
323    
324  dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)  dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)
325  {  {
# Line 435  index_t* Paso_Pattern_borrowColoringPoin Line 341  index_t* Paso_Pattern_borrowColoringPoin
341        }        }
342     }     }
343     return A->coloring;     return A->coloring;
 }  
344    }
345    dim_t Paso_Pattern_maxDeg(Paso_Pattern* A)
346    {
347       dim_t deg=0, loc_deg=0, i;
348       const dim_t n=A->numInput;
349    
350       #pragma omp parallel private(i, loc_deg)
351       {
352             loc_deg=0;
353         #pragma omp for schedule(static)
354         for (i = 0; i < n; ++i) {
355            loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);
356         }
357         #pragma omp critical
358         {
359            deg=MAX(deg, loc_deg);
360             }
361       }
362       return deg;
363    }

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

  ViewVC Help
Powered by ViewVC 1.1.26