/[escript]/trunk/dudley/src/IndexList.cpp
ViewVC logotype

Annotation of /trunk/dudley/src/IndexList.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3144 - (hide annotations)
Fri Sep 3 00:49:02 2010 UTC (9 years ago) by jfenwick
Original Path: branches/domexper/dudley/src/IndexList.c
File MIME type: text/plain
File size: 8388 byte(s)
row_node, col_node (Assemble params), node_selection 
removed
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 jfenwick 3086 /* Dudley: Converting an element list into a matrix shape */
18 jgs 82
19     /**************************************************************/
20    
21     #include "IndexList.h"
22    
23 ksteube 1312 /* Translate from distributed/local array indices to global indices */
24    
25 jgs 82 /**************************************************************/
26     /* inserts the contributions from the element matrices of elements
27     into the row index col. If symmetric is set, only the upper
28     triangle of the matrix is stored. */
29    
30 jfenwick 3086 void Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,
31 ksteube 1312 bool_t reduce_row_order, index_t* row_map,
32 jfenwick 3141 bool_t reduce_col_order, index_t* col_map)
33     {
34 ksteube 1312 /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
35 gross 2748 index_t color;
36 jfenwick 3086 Dudley_ReferenceElement*refElement;
37 jfenwick 3141 dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL;
38     if (elements!=NULL)
39     {
40 ksteube 1312 NN=elements->numNodes;
41 jfenwick 3086 refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
42 jfenwick 3141 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 gross 2748 } else {
57 jfenwick 3141 row_node=refElement->Type->subElementNodes;
58     NN_row=(refElement->BasisFunctions->Type->numShapes) ;
59 gross 2748 }
60    
61 jfenwick 3141 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 gross 2748 }
79 jfenwick 3144 }
80 jgs 82 }
81     return;
82     }
83    
84 gross 2748
85 jfenwick 3086 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 ksteube 1312 {
88 gross 2748 /* this does not resolve macro elements */
89     index_t color;
90 phornby 1628 dim_t e,kr,kc,icol,irow, NN;
91 ksteube 1312 if (elements!=NULL) {
92     NN=elements->numNodes;
93     for (color=elements->minColor;color<=elements->maxColor;color++) {
94     #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
95     for (e=0;e<elements->numElements;e++) {
96     if (elements->Color[e]==color) {
97     for (kr=0;kr<NN;kr++) {
98     irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
99     if ((firstRow<=irow) && (irow < lastRow)) {
100     irow-=firstRow;
101     for (kc=0;kc<NN;kc++) {
102     icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
103 jfenwick 3086 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
104 ksteube 1312 }
105     }
106     }
107     }
108     }
109     }
110     }
111     }
112 jfenwick 3086 void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
113     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
114 gross 1722 {
115 gross 2748 /* this does not resolve macro elements */
116 gross 1722 index_t color;
117     dim_t e,kr,kc,icol,irow, NN,irow_loc;
118     if (elements!=NULL) {
119     NN=elements->numNodes;
120     for (color=elements->minColor;color<=elements->maxColor;color++) {
121     #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
122     for (e=0;e<elements->numElements;e++) {
123     if (elements->Color[e]==color) {
124     for (kr=0;kr<NN;kr++) {
125     irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
126     if ((firstRow<=irow) && (irow < lastRow)) {
127     irow_loc=irow-firstRow;
128     for (kc=0;kc<NN;kc++) {
129     icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
130 jfenwick 3086 if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
131 gross 1722 }
132     }
133     }
134     }
135     }
136     }
137     }
138     }
139 ksteube 1312
140 jfenwick 3086 /* inserts row index row into the Dudley_IndexList in if it does not exist */
141 jgs 82
142 jfenwick 3086 void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
143 jgs 123 dim_t i;
144 jgs 82 /* is index in in? */
145     for (i=0;i<in->n;i++) {
146 jgs 102 if (in->index[i]==index) return;
147 jgs 82 }
148     /* index could not be found */
149     if (in->n==INDEXLIST_LENGTH) {
150     /* if in->index is full check the extension */
151     if (in->extension==NULL) {
152 jfenwick 3086 in->extension=TMPMEMALLOC(1,Dudley_IndexList);
153     if (Dudley_checkPtr(in->extension)) return;
154 jgs 82 in->extension->n=0;
155     in->extension->extension=NULL;
156     }
157 jfenwick 3086 Dudley_IndexList_insertIndex(in->extension,index);
158 jgs 82 } else {
159     /* insert index into in->index*/
160     in->index[in->n]=index;
161     in->n++;
162     }
163     }
164    
165 jfenwick 3086 /* counts the number of row indices in the Dudley_IndexList in */
166 jgs 82
167 jfenwick 3086 dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
168 ksteube 1312 dim_t i;
169     dim_t out=0;
170     register index_t itmp;
171 jgs 82 if (in==NULL) {
172     return 0;
173     } else {
174 ksteube 1312 for (i=0;i<in->n;i++) {
175     itmp=in->index[i];
176     if ((itmp>=range_min) && (range_max>itmp)) ++out;
177     }
178 jfenwick 3086 return out+Dudley_IndexList_count(in->extension, range_min,range_max);
179 jgs 82 }
180     }
181    
182 jfenwick 3086 /* count the number of row indices in the Dudley_IndexList in */
183 jgs 82
184 jfenwick 3086 void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
185 ksteube 1312 dim_t i, ptr;
186     register index_t itmp;
187 jgs 82 if (in!=NULL) {
188 ksteube 1312 ptr=0;
189     for (i=0;i<in->n;i++) {
190     itmp=in->index[i];
191     if ((itmp>=range_min) && (range_max>itmp)) {
192     array[ptr]=itmp+index_offset;
193     ptr++;
194     }
195    
196     }
197 jfenwick 3086 Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
198 jgs 82 }
199     }
200    
201 jfenwick 3086 /* deallocates the Dudley_IndexList in by recursive calls */
202 jgs 82
203 jfenwick 3086 void Dudley_IndexList_free(Dudley_IndexList* in) {
204 jgs 82 if (in!=NULL) {
205 jfenwick 3086 Dudley_IndexList_free(in->extension);
206 jgs 82 TMPMEMFREE(in);
207     }
208     }
209    
210 ksteube 1312 /* creates a Paso_pattern from a range of indices */
211 jfenwick 3086 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 ksteube 1312 {
213     dim_t *ptr=NULL;
214     register dim_t s,i,itmp;
215     index_t *index=NULL;
216     Paso_Pattern* out=NULL;
217    
218 gross 1552 ptr=MEMALLOC(n+1-n0,index_t);
219 jfenwick 3086 if (! Dudley_checkPtr(ptr) ) {
220 ksteube 1312 /* get the number of connections per row */
221     #pragma omp parallel for schedule(static) private(i)
222 gross 1552 for(i=n0;i<n;++i) {
223 jfenwick 3086 ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
224 ksteube 1312 }
225     /* accumulate ptr */
226     s=0;
227 gross 1552 for(i=n0;i<n;++i) {
228     itmp=ptr[i-n0];
229     ptr[i-n0]=s;
230 ksteube 1312 s+=itmp;
231     }
232 gross 1552 ptr[n-n0]=s;
233 ksteube 1312 /* fill index */
234 gross 1552 index=MEMALLOC(ptr[n-n0],index_t);
235 jfenwick 3086 if (! Dudley_checkPtr(index)) {
236 ksteube 1312 #pragma omp parallel for schedule(static)
237 gross 1552 for(i=n0;i<n;++i) {
238 jfenwick 3086 Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
239 ksteube 1312 }
240 gross 2551 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
241 ksteube 1312 }
242     }
243 jfenwick 3086 if (! Dudley_noError()) {
244 ksteube 1312 MEMFREE(ptr);
245     MEMFREE(index);
246     Paso_Pattern_free(out);
247     }
248     return out;
249     }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26