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

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

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

revision 4802 by jfenwick, Thu Feb 6 06:12:20 2014 UTC revision 4803 by caltinay, Wed Mar 26 06:52:28 2014 UTC
# Line 15  Line 15 
15  *****************************************************************************/  *****************************************************************************/
16    
17    
18    /****************************************************************************/
 /************************************************************************************/  
19    
20  /* Paso: Pattern */  /* Paso: Pattern */
21    
22  /************************************************************************************/  /****************************************************************************/
23    
24  /* Author: Lutz Gross, l.gross@uq.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
25    
26  /************************************************************************************/  /****************************************************************************/
27    
28  #include "Paso.h"  #include "Paso.h"
29  #include "Pattern.h"  #include "Pattern.h"
30    
31  /************************************************************************************/  namespace paso {
32    
33  /* allocates a Pattern  */  /* allocates a Pattern  */
34    Pattern* Pattern_alloc(int type, dim_t numOutput, dim_t numInput, index_t* ptr,
35  Paso_Pattern* Paso_Pattern_alloc(int type, dim_t numOutput, dim_t numInput, index_t* ptr, index_t* index) {                         index_t* index)
36    Paso_Pattern *out=NULL;  {
37    index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);      index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38    index_t loc_min_index,loc_max_index,min_index=index_offset,max_index=index_offset-1;      index_t loc_min_index,loc_max_index,min_index=index_offset,max_index=index_offset-1;
39    dim_t i;      dim_t i;
40    Esys_resetError();      Esys_resetError();
41    
42    if (ptr!=NULL && index != NULL) {      if (ptr!=NULL && index != NULL) {
43      #pragma omp parallel private(loc_min_index,loc_max_index,i)  #pragma omp parallel private(loc_min_index,loc_max_index,i)
44       {        {
45          loc_min_index=index_offset;          loc_min_index=index_offset;
46          loc_max_index=index_offset-1;          loc_max_index=index_offset-1;
47          if (type & MATRIX_FORMAT_OFFSET1) {          if (type & MATRIX_FORMAT_OFFSET1) {
48             #pragma omp for schedule(static)  #pragma omp for schedule(static)
49             for (i=0;i<numOutput;++i) {              for (i=0; i < numOutput; ++i) {
50                 if (ptr[i]<ptr[i+1]) {                  if (ptr[i]<ptr[i+1]) {
51                   #ifdef USE_QSORTG  #ifdef USE_QSORTG
52                      qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);                      qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t), comparIndex);
53                   #else  #else
54                      qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);                      qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t), comparIndex);
55                   #endif  #endif
56                   loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);                      loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);
57                   loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);                      loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);
58                 }                  }
59             }              }
60          } else {          } else {
61             #pragma omp for schedule(static)  #pragma omp for schedule(static)
62             for (i=0;i<numOutput;++i) {              for (i=0; i < numOutput; ++i) {
63                 if (ptr[i]<ptr[i+1]) {                  if (ptr[i] < ptr[i+1]) {
64                   #ifdef USE_QSORTG  #ifdef USE_QSORTG
65                      qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);                      qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t), comparIndex);
66                   #else  #else
67                      qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);                      qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t), comparIndex);
68                   #endif  #endif
69                   loc_min_index=MIN(loc_min_index,index[ptr[i]]);                      loc_min_index=MIN(loc_min_index,index[ptr[i]]);
70                   loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);                      loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);
71                 }                  }
72             }              }
73          }          }
74          #pragma omp critical          #pragma omp critical
75          {          {
76             min_index=MIN(loc_min_index,min_index);              min_index=MIN(loc_min_index,min_index);
77             max_index=MAX(loc_max_index,max_index);              max_index=MAX(loc_max_index,max_index);
78            }
79          } // parallel section
80    
81            if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
82                Esys_setError(TYPE_ERROR,"Pattern_alloc: Pattern index out of range.");
83                return NULL;
84          }          }
85      }      }
86      if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {      Pattern* out = new Pattern;
87        Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");      if (!Esys_checkPtr(out)) {
88        return NULL;          out->type = type;
89      }          out->reference_counter = 1;
90    }          out->numOutput = numOutput;
91    out=new Paso_Pattern;          out->numInput = numInput;
92    if (! Esys_checkPtr(out)) {          out->ptr = ptr;
93        out->type=type;          out->index = index;
94        out->reference_counter=1;          out->main_iptr = NULL;
95        out->numOutput=numOutput;          out->coloring = NULL;
96        out->numInput=numInput;          out->numColors = -1;
97        out->ptr=ptr;  
98        out->index=index;          if (out->ptr == NULL) {
99        out->main_iptr = NULL;              out->len = 0;
100        out->coloring = NULL;          } else {
101        out->numColors=-1;              out->len=out->ptr[out->numOutput] - index_offset;
102            }
103        if (out->ptr == NULL) {      }
104            out->len=0;      return out;
       } else {  
           out->len=out->ptr[out->numOutput] - index_offset;  
       }  
   }  
   return out;  
105  }  }
106    
107  /* returns a reference to in */  /* returns a reference to in */
108    Pattern* Pattern_getReference(Pattern* in)
109  Paso_Pattern* Paso_Pattern_getReference(Paso_Pattern* in) {  {
110       if (in!=NULL) {      if (in != NULL) {
111          ++(in->reference_counter);          ++(in->reference_counter);
112       }      }
113       return in;      return in;
114  }  }
115        
116  /* deallocates a Pattern: */  /* deallocates a Pattern */
117    void Pattern_free(Pattern* in)
118  void Paso_Pattern_free(Paso_Pattern* in) {  {
119    if (in!=NULL) {      if (in != NULL) {
120       in->reference_counter--;          in->reference_counter--;
121       if (in->reference_counter<=0) {          if (in->reference_counter <= 0) {
122          delete[] in->ptr;              delete[] in->ptr;
123          delete[] in->index;              delete[] in->index;
124      delete[] in->main_iptr;              delete[] in->main_iptr;
125      delete[] in->coloring;              delete[] in->coloring;
126          delete in;              delete in;
127       }          }
128     }      }
129  }  }
130  /* ***********************************************************************************/  /* **************************************************************************/
131    
132  /*  some routines which help to get the matrix pattern from elements: */  /*  some routines which help to get the matrix pattern from elements: */
133    
134  /*  this routine is used by qsort called in Paso_Pattern_alloc */  // this routine is used by qsort called in Pattern_alloc
135    int comparIndex(const void *index1, const void *index2)
136  int Paso_comparIndex(const void *index1,const void *index2){  {
137     index_t Iindex1,Iindex2;      index_t Iindex1,Iindex2;
138     Iindex1=*(index_t*)index1;      Iindex1=*(index_t*)index1;
139     Iindex2=*(index_t*)index2;      Iindex2=*(index_t*)index2;
140     if (Iindex1<Iindex2) {      if (Iindex1 < Iindex2) {
141        return -1;          return -1;
142     } else {      } else {
143        if (Iindex1>Iindex2) {          if (Iindex1 > Iindex2) {
144           return 1;              return 1;
145        } else {          } else {
146           return 0;              return 0;
147        }          }
148     }      }
149  }  }
150    
151  bool Paso_Pattern_isEmpty(Paso_Pattern* in) {  bool Pattern_isEmpty(Pattern* in)
152       if (in != NULL) {  {
153           if ((in->ptr != NULL) && (in->index != NULL)) return FALSE;      if (in != NULL) {
154       }          if ((in->ptr != NULL) && (in->index != NULL))
155       return TRUE;              return FALSE;
156        }
157        return TRUE;
158  }  }
159    
160    /* creates a pattern from a range of indices */
161  /* creates a Paso_pattern from a range of indices */  Pattern* Pattern_fromIndexListArray(dim_t n0,
162  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)          Paso_IndexListArray* index_list_array,
163            index_t range_min, index_t range_max, index_t index_offset)
164  {  {
165     const dim_t n=index_list_array->n;      Pattern* out=NULL;
166     Paso_IndexList* index_list = index_list_array->index_list;      index_t* index=NULL;
167     dim_t *ptr=NULL;      const dim_t n=index_list_array->n;
168     register dim_t s,i,itmp;      dim_t* ptr = new index_t[n+1-n0];
169     index_t *index=NULL;  
170     Paso_Pattern* out=NULL;      if (!Esys_checkPtr(ptr)) {
171            Paso_IndexList* index_list = index_list_array->index_list;
172     ptr=new index_t[n+1-n0];  
173     if (! Esys_checkPtr(ptr) ) {          // get the number of connections per row
174         /* get the number of connections per row */  #pragma omp parallel for schedule(static)
175         #pragma omp parallel for private(i) schedule(static)          for (dim_t i=n0; i < n; ++i) {
176         for(i=n0;i<n;++i) {              ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
177                ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);          }
178         }          // accumulate ptr
179         /* accumulate ptr */          dim_t s=0;
180         s=0;          for (dim_t i=n0; i < n; ++i) {
181         for(i=n0;i<n;++i) {              const dim_t itmp=ptr[i-n0];
182                 itmp=ptr[i-n0];              ptr[i-n0]=s;
183                 ptr[i-n0]=s;              s+=itmp;
184                 s+=itmp;          }
185         }          ptr[n-n0]=s;
186         ptr[n-n0]=s;          // fill index
187         /* fill index */          index=new index_t[ptr[n-n0]];
188         index=new index_t[ptr[n-n0]];          if (!Esys_checkPtr(index)) {
189         if (! Esys_checkPtr(index)) {  #pragma omp parallel for schedule(static)
190                #pragma omp parallel for private(i) schedule(static)              for (dim_t i=n0; i < n; ++i) {
191                for(i=n0;i<n;++i) {                  Paso_IndexList_toArray(&index_list[i], &index[ptr[i-n0]],
192                    Paso_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);                                         range_min, range_max, index_offset);
193                }              }
194                out=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);              out=Pattern_alloc(MATRIX_FORMAT_DEFAULT, n-n0,
195         }                                range_max+index_offset, ptr, index);
196    }          }
197    if (! Esys_noError()) {      }
198        if (!Esys_noError()) {
199          delete[] ptr;          delete[] ptr;
200          delete[] index;          delete[] index;
201          Paso_Pattern_free(out);          Pattern_free(out);
202    }      }
203    return out;      return out;
204  }  }
205    
206  index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)  index_t* Pattern_borrowMainDiagonalPointer(Pattern* A)
207  {  {
208      const dim_t n=A->numOutput;      if (A->main_iptr == NULL) {
209      int fail=0;          const dim_t n=A->numOutput;
210      index_t *index,*where_p, i;          A->main_iptr=new index_t[n];
211                if (!Esys_checkPtr(A->main_iptr)) {
212       if (A->main_iptr == NULL) {              bool fail = false;
213           A->main_iptr=new index_t[n];              // identify the main diagonals
214           if (! Esys_checkPtr(A->main_iptr) ) {  #pragma omp parallel for schedule(static)
215           #pragma omp parallel              for (index_t i = 0; i < n; ++i) {
216               {                  index_t *index=&(A->index[A->ptr[i]]);
217                   /* identify the main diagonals */                  index_t *where_p=reinterpret_cast<index_t*>(bsearch(&i,
218                   #pragma omp for schedule(static) private(i, index, where_p)                              index, (size_t)(A->ptr[i+1] - A->ptr[i]),
219                   for (i = 0; i < n; ++i) {                              sizeof(index_t), comparIndex));
220              index=&(A->index[A->ptr[i]]);                  if (where_p == NULL) {
221              where_p=reinterpret_cast<index_t*>(bsearch(&i,                      fail = true;
222                      index,                  } else {
223                      (size_t) (A->ptr[i + 1]-A->ptr[i]),                      A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);
224                           sizeof(index_t),                  }
225                           Paso_comparIndex));              }
226                                          if (fail) {
227              if (where_p==NULL) {                  delete[] A->main_iptr;
228                  fail=1;                  A->main_iptr=NULL;
229              } else {              }
230                 A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);          }
231              }      }
232                   }      return A->main_iptr;
       
              }  
          if (fail > 0) {  
            delete[] A->main_iptr;  
            A->main_iptr=NULL;  
          }  
   
      }  
      }  
      return A->main_iptr;  
233  }  }
234                                
235    dim_t Pattern_getNumColors(Pattern* A)
 dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)  
236  {  {
237     Paso_Pattern_borrowColoringPointer(A);  /* make sure numColors is defined */      // make sure numColors is defined
238     return A->numColors;      Pattern_borrowColoringPointer(A);
239        return A->numColors;
240  }  }
241  index_t* Paso_Pattern_borrowColoringPointer(Paso_Pattern* A)  
242    index_t* Pattern_borrowColoringPointer(Pattern* A)
243  {  {
244     dim_t n=A->numInput;      // is coloring available?
245     /* is coloring available ? */      if (A->coloring == NULL) {
246     if (A->coloring == NULL) {          const dim_t n = A->numInput;
247                  A->coloring = new index_t[n];
248        A->coloring=new index_t[n];          if (!Esys_checkPtr(A->coloring)) {
249        if ( ! Esys_checkPtr(A->coloring)) {              Pattern_color(A, &(A->numColors), A->coloring);
250       Paso_Pattern_color(A,&(A->numColors),A->coloring);              if (!Esys_noError()) {
251       if (! Esys_noError()) {                  delete[] A->coloring;
252          delete[] A->coloring;                  A->coloring = NULL;
253       }              }
254        }          }
255     }      }
256     return A->coloring;      return A->coloring;
257  }  }
258  dim_t Paso_Pattern_maxDeg(Paso_Pattern* A)  
259    dim_t Pattern_maxDeg(Pattern* A)
260  {  {
261     dim_t deg=0, loc_deg=0, i;      dim_t deg = 0;
262     const dim_t n=A->numInput;      const dim_t n=A->numInput;
263    
264     #pragma omp parallel private(i, loc_deg)  #pragma omp parallel
265     {      {
266           loc_deg=0;          dim_t loc_deg=0;
267       #pragma omp for schedule(static)  #pragma omp for schedule(static)
268       for (i = 0; i < n; ++i) {          for (dim_t i = 0; i < n; ++i) {
269          loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);              loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);
270       }          }
271       #pragma omp critical  #pragma omp critical
272       {          {
273          deg=MAX(deg, loc_deg);              deg = MAX(deg, loc_deg);
274           }          }
275     }      }
276     return deg;      return deg;
277  }  }
278    
279    } // namespace paso
280    

Legend:
Removed from v.4802  
changed lines
  Added in v.4803

  ViewVC Help
Powered by ViewVC 1.1.26