/[escript]/branches/domexper/dudley/src/IndexList.c
ViewVC logotype

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

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

trunk/finley/src/IndexList.c revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC branches/domexper/dudley/src/IndexList.c revision 3086 by jfenwick, Thu Aug 5 05:07:58 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  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
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) * (refElement->Type->numSides) ;
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) * (refElement->Type->numSides) ;
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) * (refElement->Type->numSides) ;
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) * (refElement->Type->numSides) ;
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        index_t color;
87    dim_t e,kr,kc,icol,irow, NN;    dim_t e,kr,kc,icol,irow, NN;
88    if (elements!=NULL) {    if (elements!=NULL) {
89      NN=elements->numNodes;      NN=elements->numNodes;
# 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 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)
209  {  {
210     dim_t *ptr=NULL;     dim_t *ptr=NULL;
211     register dim_t s,i,itmp;     register dim_t s,i,itmp;
# Line 179  Paso_Pattern* Finley_IndexList_createPat Line 213  Paso_Pattern* Finley_IndexList_createPat
213     Paso_Pattern* out=NULL;     Paso_Pattern* out=NULL;
214    
215     ptr=MEMALLOC(n+1-n0,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=n0;i<n;++i) {         for(i=n0;i<n;++i) {
220                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);
221         }         }
222         /* accumulate ptr */         /* accumulate ptr */
223         s=0;         s=0;
# Line 195  Paso_Pattern* Finley_IndexList_createPat Line 229  Paso_Pattern* Finley_IndexList_createPat
229         ptr[n-n0]=s;         ptr[n-n0]=s;
230         /* fill index */         /* fill index */
231         index=MEMALLOC(ptr[n-n0],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=n0;i<n;++i) {                for(i=n0;i<n;++i) {
235                    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);
236                }                }
237                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);
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.1628  
changed lines
  Added in v.3086

  ViewVC Help
Powered by ViewVC 1.1.26