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

trunk/finley/src/IndexList.c revision 1736 by gross, Fri Aug 29 02:23:16 2008 UTC trunk/dudley/src/IndexList.c revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1
1
/* \$Id\$ */

2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2012 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
11   *  *
12   *******************************************************/  *******************************************************/
13
14  /**************************************************************/  /**************************************************************/
15
16  /* Finley: Converting an element list into a matrix shape     */  /* Dudley: Converting an element list into a matrix shape     */
17
18  /**************************************************************/  /**************************************************************/
19
# Line 28  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 Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_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    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */  {
33    index_t color, *id=NULL;      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;      index_t color;
35    if (elements!=NULL) {      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36      NN=elements->numNodes;      if (elements != NULL)
37      id=TMPMEMALLOC(NN, index_t);      {
38      if (! Finley_checkPtr(id) ) {      NN = elements->numNodes;
39         for (i=0;i<NN;i++) id[i]=i;      NN_col = (elements->numShapes);
40         if (reduce_col_order) {      NN_row = (elements->numShapes);
41            col_node=elements->ReferenceElement->Type->linearNodes;
42            NN_col=elements->LinearReferenceElement->Type->numNodes;      for (color = elements->minColor; color <= elements->maxColor; color++)
43         } else {      {
44            col_node=id;  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45            NN_col=elements->ReferenceElement->Type->numNodes;          for (e = 0; e < elements->numElements; e++)
46         }          {
47         if (reduce_row_order) {          if (elements->Color[e] == color)
48            row_node=elements->ReferenceElement->Type->linearNodes;          {
49            NN_row=elements->LinearReferenceElement->Type->numNodes;              for (kr = 0; kr < NN_row; kr++)
50         } else {              {
51            row_node=id;              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52            NN_row=elements->ReferenceElement->Type->numNodes;              for (kc = 0; kc < NN_col; kc++)
53         }              {
54         for (color=elements->minColor;color<=elements->maxColor;color++) {                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
56             for (e=0;e<elements->numElements;e++) {              }
57                 if (elements->Color[e]==color) {              }
58                     for (kr=0;kr<NN_row;kr++) {          }
59                       irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];          }
60                       for (kc=0;kc<NN_col;kc++) {      }
61                            icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];      }
62                            Finley_IndexList_insertIndex(&(index_list[irow]),icol);      return;
63                       }  }
64                     }
65                 }  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66             }                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
67         }  {
68         TMPMEMFREE(id);  /* this does not resolve macro elements */
69      }      index_t color;
70    }      dim_t e, kr, kc, icol, irow, NN;
71    return;      if (elements != NULL)
72  }      {
73        NN = elements->numNodes;
74  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,      for (color = elements->minColor; color <= elements->maxColor; color++)
75                                                   Finley_ElementFile* elements, index_t* row_map, index_t* col_map)      {
76  {  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77    index_t color;          for (e = 0; e < elements->numElements; e++)
78    dim_t e,kr,kc,icol,irow, NN;          {
79    if (elements!=NULL) {          if (elements->Color[e] == color)
80      NN=elements->numNodes;          {
81      for (color=elements->minColor;color<=elements->maxColor;color++) {              for (kr = 0; kr < NN; kr++)
82             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)              {
83             for (e=0;e<elements->numElements;e++) {              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84                 if (elements->Color[e]==color) {              if ((firstRow <= irow) && (irow < lastRow))
85                     for (kr=0;kr<NN;kr++) {              {
86                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];                  irow -= firstRow;
87                       if ((firstRow<=irow) && (irow < lastRow)) {                  for (kc = 0; kc < NN; kc++)
88                            irow-=firstRow;                  {
89                            for (kc=0;kc<NN;kc++) {                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91                                Finley_IndexList_insertIndex(&(index_list[irow]),icol);                  }
92                            }              }
93                        }              }
94                    }          }
95                 }          }
96             }      }
97      }      }
98    }  }
99  }
100  void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101                                                                Finley_ElementFile* elements, index_t* row_map, index_t* col_map)                                     index_t lastRow, Dudley_ElementFile * elements,
102  {                                     index_t * row_map, index_t * col_map)
103    index_t color;  {
104    dim_t e,kr,kc,icol,irow, NN,irow_loc;      /* this does not resolve macro elements */
105    if (elements!=NULL) {      index_t color;
106      NN=elements->numNodes;      dim_t e, kr, kc, icol, irow, NN, irow_loc;
107      for (color=elements->minColor;color<=elements->maxColor;color++) {      if (elements != NULL)
108             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)      {
109             for (e=0;e<elements->numElements;e++) {      NN = elements->numNodes;
110                 if (elements->Color[e]==color) {      for (color = elements->minColor; color <= elements->maxColor; color++)
111                     for (kr=0;kr<NN;kr++) {      {
112                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113                       if ((firstRow<=irow) && (irow < lastRow)) {          for (e = 0; e < elements->numElements; e++)
114                            irow_loc=irow-firstRow;          {
115                            for (kc=0;kc<NN;kc++) {          if (elements->Color[e] == color)
116                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];          {
117                                if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);              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  /* inserts row index row into the Finley_IndexList in if it does not exist */                      Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128                    }
129  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {              }
130    dim_t i;              }
131    /* is index in in? */          }
132    for (i=0;i<in->n;i++) {          }
133      if (in->index[i]==index)  return;      }
134    }      }
135    /* index could not be found */  }
136    if (in->n==INDEXLIST_LENGTH) {
137       /* if in->index is full check the extension */  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138       if (in->extension==NULL) {
139          in->extension=TMPMEMALLOC(1,Finley_IndexList);  void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140          if (Finley_checkPtr(in->extension)) return;  {
141          in->extension->n=0;      dim_t i;
142          in->extension->extension=NULL;      /* is index in in? */
143       }      for (i = 0; i < in->n; i++)
144       Finley_IndexList_insertIndex(in->extension,index);      {
145    } else {      if (in->index[i] == index)
146       /* insert index into in->index*/          return;
147       in->index[in->n]=index;      }
148       in->n++;      /* index could not be found */
149    }      if (in->n == INDEXLIST_LENGTH)
150  }      {
151        /* if in->index is full check the extension */
152  /* counts the number of row indices in the Finley_IndexList in */      if (in->extension == NULL)
153        {
154  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {          in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155    dim_t i;          if (Dudley_checkPtr(in->extension))
156    dim_t out=0;          return;
157    register index_t itmp;          in->extension->n = 0;
158    if (in==NULL) {          in->extension->extension = NULL;
159       return 0;      }
160    } else {      Dudley_IndexList_insertIndex(in->extension, index);
161      for (i=0;i<in->n;i++) {      }
162            itmp=in->index[i];      else
163            if ((itmp>=range_min) && (range_max>itmp)) ++out;      {
164      }      /* insert index into in->index */
165       return out+Finley_IndexList_count(in->extension, range_min,range_max);      in->index[in->n] = index;
166    }      in->n++;
167  }      }
168    }
169  /* count the number of row indices in the Finley_IndexList in */
170    /* counts the number of row indices in the Dudley_IndexList in */
171  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
172    dim_t i, ptr;  dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173    register index_t itmp;  {
174    if (in!=NULL) {      dim_t i;
175      ptr=0;      dim_t out = 0;
176      for (i=0;i<in->n;i++) {      register index_t itmp;
177            itmp=in->index[i];      if (in == NULL)
178            if ((itmp>=range_min) && (range_max>itmp)) {      {
179               array[ptr]=itmp+index_offset;      return 0;
180               ptr++;      }
181            }      else
182        {
183      }      for (i = 0; i < in->n; i++)
184      Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);      {
185    }          itmp = in->index[i];
186  }          if ((itmp >= range_min) && (range_max > itmp))
187            ++out;
188  /* deallocates the Finley_IndexList in by recursive calls */      }
189        return out + Dudley_IndexList_count(in->extension, range_min, range_max);
190  void Finley_IndexList_free(Finley_IndexList* in) {      }
191    if (in!=NULL) {  }
192      Finley_IndexList_free(in->extension);
193      TMPMEMFREE(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,
196                      index_t index_offset)
197    {
198        dim_t i, ptr;
199        register index_t itmp;
200        if (in != NULL)
201        {
202        ptr = 0;
203        for (i = 0; i < in->n; i++)
204        {
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        }
215    }
216
217    /* deallocates the Dudley_IndexList in by recursive calls */
218
219    void Dudley_IndexList_free(Dudley_IndexList * in)
220    {
221        if (in != NULL)
222        {
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* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_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 (! Finley_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]=Finley_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 (! Finley_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                    Finley_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,1,1,n-n0,range_max+index_offset,ptr,index);  #pragma omp parallel for schedule(static)
260         }          for (i = n0; i < n; ++i)
261    }          {
262    if (! Finley_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.1736 changed lines Added in v.3911