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

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

  ViewVC Help
Powered by ViewVC 1.1.26