/[escript]/trunk/dudley/src/IndexList.cpp
ViewVC logotype

Diff of /trunk/dudley/src/IndexList.cpp

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

branches/domexper/dudley/src/IndexList.c revision 3152 by jfenwick, Fri Sep 3 05:48:31 2010 UTC trunk/dudley/src/IndexList.c revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
   
14  /**************************************************************/  /**************************************************************/
15    
16  /* Dudley: Converting an element list into a matrix shape     */  /* Dudley: Converting an element list into a matrix shape     */
# Line 27  Line 26 
26     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
27     triangle of the matrix is stored. */     triangle of the matrix is stored. */
28    
29  void Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,  void Dudley_IndexList_insertElements(Dudley_IndexList * index_list, Dudley_ElementFile * elements,
30                                         bool_t reduce_row_order, index_t* row_map,                       bool_t reduce_row_order, index_t * row_map,
31                                         bool_t reduce_col_order, index_t* col_map)                       bool_t reduce_col_order, index_t * col_map)
32  {  {
33    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34    index_t color;      index_t color;
35    Dudley_ReferenceElement*refElement;      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36    dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;      if (elements != NULL)
   if (elements!=NULL)  
   {  
     NN=elements->numNodes;  
     refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);  
     NN_col=(refElement->BasisFunctions->Type->numShapes);  
     NN_col=(refElement->BasisFunctions->Type->numShapes);  
     NN_row=(refElement->BasisFunctions->Type->numShapes) ;  
   
     for (color=elements->minColor;color<=elements->maxColor;color++)  
37      {      {
38      #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)      NN = elements->numNodes;
39      for (e=0;e<elements->numElements;e++)      NN_col = (elements->numShapes);
40        NN_row = (elements->numShapes);
41    
42        for (color = elements->minColor; color <= elements->maxColor; color++)
43      {      {
44          if (elements->Color[e]==color)  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45            for (e = 0; e < elements->numElements; e++)
46            {
47            if (elements->Color[e] == color)
48          {          {
49              for (kr=0;kr<NN_row;kr++)              for (kr = 0; kr < NN_row; kr++)
50                {
51                irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52                for (kc = 0; kc < NN_col; kc++)
53              {              {
54                  irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55                  for (kc=0;kc<NN_col;kc++)                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
                 {  
                 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];  
                 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);  
                 }  
56              }              }
57                }
58          }          }
59            }
60      }      }
61       }      }
62    }      return;
   return;  
63  }  }
64    
65    void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
                                                  Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)  
67  {  {
68  /* this does not resolve macro elements */  /* this does not resolve macro elements */
69      index_t color;      index_t color;
70    dim_t e,kr,kc,icol,irow, NN;      dim_t e, kr, kc, icol, irow, NN;
71    if (elements!=NULL) {      if (elements != NULL)
72      NN=elements->numNodes;      {
73      for (color=elements->minColor;color<=elements->maxColor;color++) {      NN = elements->numNodes;
74             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)      for (color = elements->minColor; color <= elements->maxColor; color++)
75             for (e=0;e<elements->numElements;e++) {      {
76                 if (elements->Color[e]==color) {  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77                     for (kr=0;kr<NN;kr++) {          for (e = 0; e < elements->numElements; e++)
78                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];          {
79                       if ((firstRow<=irow) && (irow < lastRow)) {          if (elements->Color[e] == color)
80                            irow-=firstRow;          {
81                            for (kc=0;kc<NN;kc++) {              for (kr = 0; kr < NN; kr++)
82                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];              {
83                                Dudley_IndexList_insertIndex(&(index_list[irow]),icol);              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84                            }              if ((firstRow <= irow) && (irow < lastRow))
85                        }              {
86                    }                  irow -= firstRow;
87                 }                  for (kc = 0; kc < NN; kc++)
88             }                  {
89      }                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90    }                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91  }                  }
92  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,              }
93                                                                Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)              }
94  {          }
95    /* this does not resolve macro elements */          }
96    index_t color;      }
97    dim_t e,kr,kc,icol,irow, NN,irow_loc;      }
98    if (elements!=NULL) {  }
99      NN=elements->numNodes;  
100      for (color=elements->minColor;color<=elements->maxColor;color++) {  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)                                     index_t lastRow, Dudley_ElementFile * elements,
102             for (e=0;e<elements->numElements;e++) {                                     index_t * row_map, index_t * col_map)
103                 if (elements->Color[e]==color) {  {
104                     for (kr=0;kr<NN;kr++) {      /* this does not resolve macro elements */
105                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];      index_t color;
106                       if ((firstRow<=irow) && (irow < lastRow)) {      dim_t e, kr, kc, icol, irow, NN, irow_loc;
107                            irow_loc=irow-firstRow;      if (elements != NULL)
108                            for (kc=0;kc<NN;kc++) {      {
109                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];      NN = elements->numNodes;
110                                if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);      for (color = elements->minColor; color <= elements->maxColor; color++)
111                            }      {
112                        }  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113                    }          for (e = 0; e < elements->numElements; e++)
114                 }          {
115             }          if (elements->Color[e] == color)
116            {
117                for (kr = 0; kr < NN; kr++)
118                {
119                irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
120                if ((firstRow <= irow) && (irow < lastRow))
121                {
122                    irow_loc = irow - firstRow;
123                    for (kc = 0; kc < NN; kc++)
124                    {
125                    icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
126                    if (icol != irow)
127                        Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128                    }
129                }
130                }
131            }
132            }
133        }
134      }      }
   }  
135  }  }
136    
137  /* inserts row index row into the Dudley_IndexList in if it does not exist */  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138    
139  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140    dim_t i;  {
141    /* is index in in? */      dim_t i;
142    for (i=0;i<in->n;i++) {      /* is index in in? */
143      if (in->index[i]==index)  return;      for (i = 0; i < in->n; i++)
144    }      {
145    /* index could not be found */      if (in->index[i] == index)
146    if (in->n==INDEXLIST_LENGTH) {          return;
147       /* if in->index is full check the extension */      }
148       if (in->extension==NULL) {      /* index could not be found */
149          in->extension=TMPMEMALLOC(1,Dudley_IndexList);      if (in->n == INDEXLIST_LENGTH)
150          if (Dudley_checkPtr(in->extension)) return;      {
151          in->extension->n=0;      /* if in->index is full check the extension */
152          in->extension->extension=NULL;      if (in->extension == NULL)
153       }      {
154       Dudley_IndexList_insertIndex(in->extension,index);          in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155    } else {          if (Dudley_checkPtr(in->extension))
156       /* insert index into in->index*/          return;
157       in->index[in->n]=index;          in->extension->n = 0;
158       in->n++;          in->extension->extension = NULL;
159    }      }
160        Dudley_IndexList_insertIndex(in->extension, index);
161        }
162        else
163        {
164        /* insert index into in->index */
165        in->index[in->n] = index;
166        in->n++;
167        }
168  }  }
169    
170  /* counts the number of row indices in the Dudley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
171    
172  dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {  dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173    dim_t i;  {
174    dim_t out=0;      dim_t i;
175    register index_t itmp;      dim_t out = 0;
176    if (in==NULL) {      register index_t itmp;
177       return 0;      if (in == NULL)
178    } else {      {
179      for (i=0;i<in->n;i++) {      return 0;
180            itmp=in->index[i];      }
181            if ((itmp>=range_min) && (range_max>itmp)) ++out;      else
182        {
183        for (i = 0; i < in->n; i++)
184        {
185            itmp = in->index[i];
186            if ((itmp >= range_min) && (range_max > itmp))
187            ++out;
188        }
189        return out + Dudley_IndexList_count(in->extension, range_min, range_max);
190      }      }
      return out+Dudley_IndexList_count(in->extension, range_min,range_max);  
   }  
191  }  }
192    
193  /* count the number of row indices in the Dudley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
194    
195  void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {  void Dudley_IndexList_toArray(Dudley_IndexList * in, index_t * array, index_t range_min, index_t range_max,
196    dim_t i, ptr;                    index_t index_offset)
197    register index_t itmp;  {
198    if (in!=NULL) {      dim_t i, ptr;
199      ptr=0;      register index_t itmp;
200      for (i=0;i<in->n;i++) {      if (in != NULL)
201            itmp=in->index[i];      {
202            if ((itmp>=range_min) && (range_max>itmp)) {      ptr = 0;
203               array[ptr]=itmp+index_offset;      for (i = 0; i < in->n; i++)
204               ptr++;      {
205            }          itmp = in->index[i];
206            if ((itmp >= range_min) && (range_max > itmp))
207            {
208            array[ptr] = itmp + index_offset;
209            ptr++;
210            }
211    
212        }
213        Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
214      }      }
     Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);  
   }  
215  }  }
216    
217  /* deallocates the Dudley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
218    
219  void Dudley_IndexList_free(Dudley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList * in)
220    if (in!=NULL) {  {
221      Dudley_IndexList_free(in->extension);      if (in != NULL)
222      TMPMEMFREE(in);      {
223    }      Dudley_IndexList_free(in->extension);
224        TMPMEMFREE(in);
225        }
226  }  }
227    
228  /* creates a Paso_pattern from a range of indices */  /* creates a Paso_pattern from a range of indices */
229  Paso_Pattern* Dudley_IndexList_createPattern(dim_t n0, dim_t n,Dudley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)  Paso_Pattern *Dudley_IndexList_createPattern(dim_t n0, dim_t n, Dudley_IndexList * index_list, index_t range_min,
230                             index_t range_max, index_t index_offset)
231  {  {
232     dim_t *ptr=NULL;      dim_t *ptr = NULL;
233     register dim_t s,i,itmp;      register dim_t s, i, itmp;
234     index_t *index=NULL;      index_t *index = NULL;
235     Paso_Pattern* out=NULL;      Paso_Pattern *out = NULL;
236    
237     ptr=MEMALLOC(n+1-n0,index_t);      ptr = MEMALLOC(n + 1 - n0, index_t);
238     if (! Dudley_checkPtr(ptr) ) {      if (!Dudley_checkPtr(ptr))
239         /* get the number of connections per row */      {
240         #pragma omp parallel for schedule(static) private(i)      /* get the number of connections per row */
241         for(i=n0;i<n;++i) {  #pragma omp parallel for schedule(static) private(i)
242                ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);      for (i = n0; i < n; ++i)
243         }      {
244         /* accumulate ptr */          ptr[i - n0] = Dudley_IndexList_count(&index_list[i], range_min, range_max);
245         s=0;      }
246         for(i=n0;i<n;++i) {      /* accumulate ptr */
247                 itmp=ptr[i-n0];      s = 0;
248                 ptr[i-n0]=s;      for (i = n0; i < n; ++i)
249                 s+=itmp;      {
250         }          itmp = ptr[i - n0];
251         ptr[n-n0]=s;          ptr[i - n0] = s;
252         /* fill index */          s += itmp;
253         index=MEMALLOC(ptr[n-n0],index_t);      }
254         if (! Dudley_checkPtr(index)) {      ptr[n - n0] = s;
255                #pragma omp parallel for schedule(static)      /* fill index */
256                for(i=n0;i<n;++i) {      index = MEMALLOC(ptr[n - n0], index_t);
257                    Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);      if (!Dudley_checkPtr(index))
258                }      {
259                out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);  #pragma omp parallel for schedule(static)
260         }          for (i = n0; i < n; ++i)
261    }          {
262    if (! Dudley_noError()) {          Dudley_IndexList_toArray(&index_list[i], &index[ptr[i - n0]], range_min, range_max, index_offset);
263          MEMFREE(ptr);          }
264          MEMFREE(index);          out = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n - n0, range_max + index_offset, ptr, index);
265          Paso_Pattern_free(out);      }
266    }      }
267    return out;      if (!Dudley_noError())
268        {
269        MEMFREE(ptr);
270        MEMFREE(index);
271        Paso_Pattern_free(out);
272        }
273        return out;
274  }  }

Legend:
Removed from v.3152  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26