# Diff of /branches/domexper/dudley/src/IndexList.c

trunk/finley/src/IndexList.c revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC branches/domexper/dudley/src/IndexList.c revision 3141 by jfenwick, Thu Sep 2 23:52:41 2010 UTC
# Line 1  Line 1
1
/* \$Id\$ */

2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 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
17  /* Finley: Converting an element list into a matrix shape     */  /* Dudley: Converting an element list into a matrix shape     */
18
19  /**************************************************************/  /**************************************************************/
20
# Line 28  Line 27
27     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
28     triangle of the matrix is stored. */     triangle of the matrix is stored. */
29
30  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,  void Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,
31                                         bool_t reduce_row_order, index_t* row_map,                                         bool_t reduce_row_order, index_t* row_map,
32                                         bool_t reduce_col_order, index_t* col_map) {                                         bool_t reduce_col_order, index_t* col_map)
33    {
34    /* 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 */
35    index_t color, *id=NULL;    index_t color;
36    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;    Dudley_ReferenceElement*refElement;
37      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL;
38      if (elements!=NULL)
39      {
40        NN=elements->numNodes;
41        refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
42        if (reduce_col_order)
43        {
44        col_node=refElement->Type->linearNodes;
45        NN_col=(refElement->LinearBasisFunctions->Type->numShapes);
46        }
47        else
48        {
49        col_node=refElement->Type->subElementNodes;
50        NN_col=(refElement->BasisFunctions->Type->numShapes);
51        }
52        if (reduce_row_order)
53        {
54        row_node=refElement->Type->linearNodes;
55        NN_row=(refElement->LinearBasisFunctions->Type->numShapes);
56        } else {
57        row_node=refElement->Type->subElementNodes;
58        NN_row=(refElement->BasisFunctions->Type->numShapes) ;
59        }
60
61        for (color=elements->minColor;color<=elements->maxColor;color++)
62        {
63        #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
64        for (e=0;e<elements->numElements;e++)
65        {
66            if (elements->Color[e]==color)
67            {
68                for (kr=0;kr<NN_row;kr++)
69                {
70                    irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,0,NN_row)],e,NN)]];
71                    for (kc=0;kc<NN_col;kc++)
72                    {
73                    icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,0,NN_col)],e,NN)]];
74                    Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
75                    }
76                }
77            }
78            }
79        }
80      }
81      return;
82    }
83
84
85    void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
86                                                     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
87    {
88    /* this does not resolve macro elements */
89        index_t color;
90      dim_t e,kr,kc,icol,irow, NN;
91    if (elements!=NULL) {    if (elements!=NULL) {
92      NN=elements->numNodes;      NN=elements->numNodes;
93      id=TMPMEMALLOC(NN, index_t);      for (color=elements->minColor;color<=elements->maxColor;color++) {
if (! Finley_checkPtr(id) ) {
for (i=0;i<NN;i++) id[i]=i;
if (reduce_col_order) {
col_node=elements->ReferenceElement->Type->linearNodes;
NN_col=elements->LinearReferenceElement->Type->numNodes;
} else {
col_node=id;
NN_col=elements->ReferenceElement->Type->numNodes;
}
if (reduce_row_order) {
row_node=elements->ReferenceElement->Type->linearNodes;
NN_row=elements->LinearReferenceElement->Type->numNodes;
} else {
row_node=id;
NN_row=elements->ReferenceElement->Type->numNodes;
}
for (color=elements->minColor;color<=elements->maxColor;color++) {
94             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
95             for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
96                 if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
97                     for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
98                       irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
99                       for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
100                            icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow-=firstRow;
101                            Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
102                       }                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
103                     }                                Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
104                              }
105                          }
106                      }
107                 }                 }
108             }             }
}
TMPMEMFREE(id);
109      }      }
110    }    }
return;
111  }  }
112    void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
113  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,                                                                Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
114  {  {
115      /* this does not resolve macro elements */
116    index_t color;    index_t color;
117    dim_t e,kr,kc,icol,irow, NN;    dim_t e,kr,kc,icol,irow, NN,irow_loc;
118    if (elements!=NULL) {    if (elements!=NULL) {
119      NN=elements->numNodes;      NN=elements->numNodes;
120      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
121             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
122             for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
123                 if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
124                     for (kr=0;kr<NN;kr++) {                     for (kr=0;kr<NN;kr++) {
125                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
126                       if ((firstRow<=irow) && (irow < lastRow)) {                       if ((firstRow<=irow) && (irow < lastRow)) {
127                            irow-=firstRow;                            irow_loc=irow-firstRow;
128                            for (kc=0;kc<NN;kc++) {                            for (kc=0;kc<NN;kc++) {
129                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
130                                Finley_IndexList_insertIndex(&(index_list[irow]),icol);                                if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
131                            }                            }
132                        }                        }
133                    }                    }
# Line 100  void Finley_IndexList_insertElementsWith Line 137  void Finley_IndexList_insertElementsWith
137    }    }
138  }  }
139
140  /* inserts row index row into the Finley_IndexList in if it does not exist */  /* inserts row index row into the Dudley_IndexList in if it does not exist */
141
142  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
143    dim_t i;    dim_t i;
144    /* is index in in? */    /* is index in in? */
145    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
# Line 112  void Finley_IndexList_insertIndex(Finley Line 149  void Finley_IndexList_insertIndex(Finley
149    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
150       /* if in->index is full check the extension */       /* if in->index is full check the extension */
151       if (in->extension==NULL) {       if (in->extension==NULL) {
152          in->extension=TMPMEMALLOC(1,Finley_IndexList);          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
153          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
154          in->extension->n=0;          in->extension->n=0;
155          in->extension->extension=NULL;          in->extension->extension=NULL;
156       }       }
157       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
158    } else {    } else {
159       /* insert index into in->index*/       /* insert index into in->index*/
160       in->index[in->n]=index;       in->index[in->n]=index;
# Line 125  void Finley_IndexList_insertIndex(Finley Line 162  void Finley_IndexList_insertIndex(Finley
162    }    }
163  }  }
164
165  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
166
167  dim_t Finley_IndexList_count(Finley_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) {
168    dim_t i;    dim_t i;
169    dim_t out=0;    dim_t out=0;
170    register index_t itmp;    register index_t itmp;
# Line 138  dim_t Finley_IndexList_count(Finley_Inde Line 175  dim_t Finley_IndexList_count(Finley_Inde
175            itmp=in->index[i];            itmp=in->index[i];
176            if ((itmp>=range_min) && (range_max>itmp)) ++out;            if ((itmp>=range_min) && (range_max>itmp)) ++out;
177      }      }
178       return out+Finley_IndexList_count(in->extension, range_min,range_max);       return out+Dudley_IndexList_count(in->extension, range_min,range_max);
179    }    }
180  }  }
181
182  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
183
184  void Finley_IndexList_toArray(Finley_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, index_t index_offset) {
185    dim_t i, ptr;    dim_t i, ptr;
186    register index_t itmp;    register index_t itmp;
187    if (in!=NULL) {    if (in!=NULL) {
# Line 157  void Finley_IndexList_toArray(Finley_Ind Line 194  void Finley_IndexList_toArray(Finley_Ind
194            }            }
195
196      }      }
197      Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);      Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
198    }    }
199  }  }
200
201  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
202
203  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
204    if (in!=NULL) {    if (in!=NULL) {
205      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
206      TMPMEMFREE(in);      TMPMEMFREE(in);
207    }    }
208  }  }
209
210  /* creates a Paso_pattern from a range of indices */  /* creates a Paso_pattern from a range of indices */
211  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,index_t range_max,index_t index_offset)
212  {  {
213     dim_t *ptr=NULL;     dim_t *ptr=NULL;
214     register dim_t s,i,itmp;     register dim_t s,i,itmp;
# Line 179  Paso_Pattern* Finley_IndexList_createPat Line 216  Paso_Pattern* Finley_IndexList_createPat
216     Paso_Pattern* out=NULL;     Paso_Pattern* out=NULL;
217
218     ptr=MEMALLOC(n+1-n0,index_t);     ptr=MEMALLOC(n+1-n0,index_t);
219     if (! Finley_checkPtr(ptr) ) {     if (! Dudley_checkPtr(ptr) ) {
220         /* get the number of connections per row */         /* get the number of connections per row */
221         #pragma omp parallel for schedule(static) private(i)         #pragma omp parallel for schedule(static) private(i)
222         for(i=n0;i<n;++i) {         for(i=n0;i<n;++i) {
223                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);                ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
224         }         }
225         /* accumulate ptr */         /* accumulate ptr */
226         s=0;         s=0;
# Line 195  Paso_Pattern* Finley_IndexList_createPat Line 232  Paso_Pattern* Finley_IndexList_createPat
232         ptr[n-n0]=s;         ptr[n-n0]=s;
233         /* fill index */         /* fill index */
234         index=MEMALLOC(ptr[n-n0],index_t);         index=MEMALLOC(ptr[n-n0],index_t);
235         if (! Finley_checkPtr(index)) {         if (! Dudley_checkPtr(index)) {
236                #pragma omp parallel for schedule(static)                #pragma omp parallel for schedule(static)
237                for(i=n0;i<n;++i) {                for(i=n0;i<n;++i) {
238                    Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);                    Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
239                }                }
240                out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,ptr,index);                out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
241         }         }
242    }    }
243    if (! Finley_noError()) {    if (! Dudley_noError()) {
244          MEMFREE(ptr);          MEMFREE(ptr);
245          MEMFREE(index);          MEMFREE(index);
246          Paso_Pattern_free(out);          Paso_Pattern_free(out);

Legend:
 Removed from v.1628 changed lines Added in v.3141