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

temp/finley/src/IndexList.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC branches/domexper/dudley/src/IndexList.c revision 3127 by jfenwick, Wed Sep 1 01:30:28 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    /* 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, *id=NULL;    index_t color;
35    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;    Dudley_ReferenceElement*refElement;
36      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL, isub, numSub;
37    if (elements!=NULL) {    if (elements!=NULL) {
38      NN=elements->numNodes;      NN=elements->numNodes;
39      id=TMPMEMALLOC(NN, index_t);      refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
40      if (! Finley_checkPtr(id) ) {
41         for (i=0;i<NN;i++) id[i]=i;      if (reduce_col_order) {
42         if (reduce_col_order) {            numSub=1;
43            col_node=elements->ReferenceElement->Type->linearNodes;            col_node=refElement->Type->linearNodes;
44            NN_col=elements->LinearReferenceElement->Type->numNodes;            NN_col=(refElement->LinearBasisFunctions->Type->numShapes);
45         } else {      } else {
46            col_node=id;            numSub=refElement->Type->numSubElements;
47            NN_col=elements->ReferenceElement->Type->numNodes;            col_node=refElement->Type->subElementNodes;
48         }            NN_col=(refElement->BasisFunctions->Type->numShapes);
49         if (reduce_row_order) {      }
50            row_node=elements->ReferenceElement->Type->linearNodes;
51            NN_row=elements->LinearReferenceElement->Type->numNodes;      if (reduce_row_order) {
52         } else {            numSub=1;
53            row_node=id;            row_node=refElement->Type->linearNodes;
54            NN_row=elements->ReferenceElement->Type->numNodes;            NN_row=(refElement->LinearBasisFunctions->Type->numShapes);
55         }      } else {
56         for (color=elements->minColor;color<=elements->maxColor;color++) {            numSub=refElement->Type->numSubElements;
57             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)            row_node=refElement->Type->subElementNodes;
58              NN_row=(refElement->BasisFunctions->Type->numShapes) ;
59        }
60
61        for (color=elements->minColor;color<=elements->maxColor;color++) {
62               #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
63             for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
64                 if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
65                     for (kr=0;kr<NN_row;kr++) {                     for (isub=0;isub<numSub; isub++) {
66                       irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                         for (kr=0;kr<NN_row;kr++) {
67                       for (kc=0;kc<NN_col;kc++) {                             irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
68                            icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                             for (kc=0;kc<NN_col;kc++) {
69                            Finley_IndexList_insertIndex(&(index_list[irow]),icol);                                 icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
70                       }                                 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
71                     }                             }
72                           }
73                       }
74                 }                 }
75             }             }
76         }         }
TMPMEMFREE(id);
}
77    }    }
78    return;    return;
79  }  }
80
81  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
82                                                   Finley_ElementFile* elements, index_t* row_map, index_t* col_map)  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
83                                                     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
84  {  {
85    index_t color;  /* this does not resolve macro elements */
86    dim_t e,kr,kc,i,icol,irow, NN;      index_t color;
87      dim_t e,kr,kc,icol,irow, NN;
88    if (elements!=NULL) {    if (elements!=NULL) {
89      NN=elements->numNodes;      NN=elements->numNodes;
90      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
# Line 90  void Finley_IndexList_insertElementsWith Line 97  void Finley_IndexList_insertElementsWith
97                            irow-=firstRow;                            irow-=firstRow;
98                            for (kc=0;kc<NN;kc++) {                            for (kc=0;kc<NN;kc++) {
99                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
100                                Finley_IndexList_insertIndex(&(index_list[irow]),icol);                                Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
101                              }
102                          }
103                      }
104                   }
105               }
106        }
107      }
108    }
109    void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
110                                                                  Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
111    {
112      /* this does not resolve macro elements */
113      index_t color;
114      dim_t e,kr,kc,icol,irow, NN,irow_loc;
115      if (elements!=NULL) {
116        NN=elements->numNodes;
117        for (color=elements->minColor;color<=elements->maxColor;color++) {
118               #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
119               for (e=0;e<elements->numElements;e++) {
120                   if (elements->Color[e]==color) {
121                       for (kr=0;kr<NN;kr++) {
122                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
123                         if ((firstRow<=irow) && (irow < lastRow)) {
124                              irow_loc=irow-firstRow;
125                              for (kc=0;kc<NN;kc++) {
126                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
127                                  if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
128                            }                            }
129                        }                        }
130                    }                    }
# Line 100  void Finley_IndexList_insertElementsWith Line 134  void Finley_IndexList_insertElementsWith
134    }    }
135  }  }
136
137  /* 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 */
138
139  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
140    dim_t i;    dim_t i;
141    /* is index in in? */    /* is index in in? */
142    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
# Line 112  void Finley_IndexList_insertIndex(Finley Line 146  void Finley_IndexList_insertIndex(Finley
146    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
147       /* if in->index is full check the extension */       /* if in->index is full check the extension */
148       if (in->extension==NULL) {       if (in->extension==NULL) {
149          in->extension=TMPMEMALLOC(1,Finley_IndexList);          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
150          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
151          in->extension->n=0;          in->extension->n=0;
152          in->extension->extension=NULL;          in->extension->extension=NULL;
153       }       }
154       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
155    } else {    } else {
156       /* insert index into in->index*/       /* insert index into in->index*/
157       in->index[in->n]=index;       in->index[in->n]=index;
# Line 125  void Finley_IndexList_insertIndex(Finley Line 159  void Finley_IndexList_insertIndex(Finley
159    }    }
160  }  }
161
162  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
163
164  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) {
165    dim_t i;    dim_t i;
166    dim_t out=0;    dim_t out=0;
167    register index_t itmp;    register index_t itmp;
# Line 138  dim_t Finley_IndexList_count(Finley_Inde Line 172  dim_t Finley_IndexList_count(Finley_Inde
172            itmp=in->index[i];            itmp=in->index[i];
173            if ((itmp>=range_min) && (range_max>itmp)) ++out;            if ((itmp>=range_min) && (range_max>itmp)) ++out;
174      }      }
175       return out+Finley_IndexList_count(in->extension, range_min,range_max);       return out+Dudley_IndexList_count(in->extension, range_min,range_max);
176    }    }
177  }  }
178
179  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
180
181  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) {
182    dim_t i, ptr;    dim_t i, ptr;
183    register index_t itmp;    register index_t itmp;
184    if (in!=NULL) {    if (in!=NULL) {
# Line 157  void Finley_IndexList_toArray(Finley_Ind Line 191  void Finley_IndexList_toArray(Finley_Ind
191            }            }
192
193      }      }
194      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);
195    }    }
196  }  }
197
198  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
199
200  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
201    if (in!=NULL) {    if (in!=NULL) {
202      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
203      TMPMEMFREE(in);      TMPMEMFREE(in);
204    }    }
205  }  }
206
207  /* creates a Paso_pattern from a range of indices */  /* creates a Paso_pattern from a range of indices */
208  Paso_Pattern* Finley_IndexList_createPattern(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)
209  {  {
210     dim_t *ptr=NULL;     dim_t *ptr=NULL;
211     register dim_t s,i,itmp;     register dim_t s,i,itmp;
212     index_t *index=NULL;     index_t *index=NULL;
213     Paso_Pattern* out=NULL;     Paso_Pattern* out=NULL;
214
215     ptr=MEMALLOC(n+1,index_t);     ptr=MEMALLOC(n+1-n0,index_t);
216     if (! Finley_checkPtr(ptr) ) {     if (! Dudley_checkPtr(ptr) ) {
217         /* get the number of connections per row */         /* get the number of connections per row */
218         #pragma omp parallel for schedule(static) private(i)         #pragma omp parallel for schedule(static) private(i)
219         for(i=0;i<n;++i) {         for(i=n0;i<n;++i) {
220                ptr[i]=Finley_IndexList_count(&index_list[i],range_min,range_max);                ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
221         }         }
222         /* accumulate ptr */         /* accumulate ptr */
223         s=0;         s=0;
224         for(i=0;i<n;++i) {         for(i=n0;i<n;++i) {
225                 itmp=ptr[i];                 itmp=ptr[i-n0];
226                 ptr[i]=s;                 ptr[i-n0]=s;
227                 s+=itmp;                 s+=itmp;
228         }         }
229         ptr[n]=s;         ptr[n-n0]=s;
230         /* fill index */         /* fill index */
231         index=MEMALLOC(ptr[n],index_t);         index=MEMALLOC(ptr[n-n0],index_t);
232         if (! Finley_checkPtr(index)) {         if (! Dudley_checkPtr(index)) {
233                #pragma omp parallel for schedule(static)                #pragma omp parallel for schedule(static)
234                for(i=0;i<n;++i) {                for(i=n0;i<n;++i) {
235                    Finley_IndexList_toArray(&index_list[i],&index[ptr[i]],range_min,range_max,index_offset);                    Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
236                }                }
237                out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n,ptr,index);                out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
238         }         }
239    }    }
240    if (! Finley_noError()) {    if (! Dudley_noError()) {
241          MEMFREE(ptr);          MEMFREE(ptr);
242          MEMFREE(index);          MEMFREE(index);
243          Paso_Pattern_free(out);          Paso_Pattern_free(out);

Legend:
 Removed from v.1387 changed lines Added in v.3127