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

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

  ViewVC Help
Powered by ViewVC 1.1.26