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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 3 months ago) by jfenwick
Original Path: trunk/dudley/src/IndexList.c
File MIME type: text/plain
File size: 7300 byte(s)
Copyright changes
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 3911 * Copyright (c) 2003-2012 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     /**************************************************************/
15    
16 jfenwick 3086 /* Dudley: Converting an element list into a matrix shape */
17 jgs 82
18     /**************************************************************/
19    
20     #include "IndexList.h"
21    
22 ksteube 1312 /* Translate from distributed/local array indices to global indices */
23    
24 jgs 82 /**************************************************************/
25     /* inserts the contributions from the element matrices of elements
26     into the row index col. If symmetric is set, only the upper
27     triangle of the matrix is stored. */
28    
29 jfenwick 3224 void Dudley_IndexList_insertElements(Dudley_IndexList * index_list, Dudley_ElementFile * elements,
30     bool_t reduce_row_order, index_t * row_map,
31     bool_t reduce_col_order, index_t * col_map)
32 jfenwick 3141 {
33 jfenwick 3224 /* 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;
35     dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36     if (elements != NULL)
37     {
38     NN = elements->numNodes;
39     NN_col = (elements->numShapes);
40     NN_row = (elements->numShapes);
41 gross 2748
42 jfenwick 3224 for (color = elements->minColor; color <= elements->maxColor; color++)
43 jfenwick 3141 {
44 jfenwick 3224 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45     for (e = 0; e < elements->numElements; e++)
46     {
47     if (elements->Color[e] == color)
48 jfenwick 3141 {
49 jfenwick 3224 for (kr = 0; kr < NN_row; kr++)
50     {
51     irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52     for (kc = 0; kc < NN_col; kc++)
53 jfenwick 3141 {
54 jfenwick 3224 icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55     Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
56 jfenwick 3141 }
57 jfenwick 3224 }
58 jfenwick 3141 }
59 jfenwick 3224 }
60 gross 2748 }
61 jfenwick 3224 }
62     return;
63 jgs 82 }
64    
65 jfenwick 3224 void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66     Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
67 ksteube 1312 {
68 gross 2748 /* this does not resolve macro elements */
69 jfenwick 3224 index_t color;
70     dim_t e, kr, kc, icol, irow, NN;
71     if (elements != NULL)
72     {
73     NN = elements->numNodes;
74     for (color = elements->minColor; color <= elements->maxColor; color++)
75     {
76     #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77     for (e = 0; e < elements->numElements; e++)
78     {
79     if (elements->Color[e] == color)
80     {
81     for (kr = 0; kr < NN; kr++)
82     {
83     irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84     if ((firstRow <= irow) && (irow < lastRow))
85     {
86     irow -= firstRow;
87     for (kc = 0; kc < NN; kc++)
88     {
89     icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90     Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91     }
92     }
93     }
94     }
95     }
96     }
97 ksteube 1312 }
98     }
99 jfenwick 3224
100     void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101     index_t lastRow, Dudley_ElementFile * elements,
102     index_t * row_map, index_t * col_map)
103 gross 1722 {
104 jfenwick 3224 /* this does not resolve macro elements */
105     index_t color;
106     dim_t e, kr, kc, icol, irow, NN, irow_loc;
107     if (elements != NULL)
108     {
109     NN = elements->numNodes;
110     for (color = elements->minColor; color <= elements->maxColor; color++)
111     {
112     #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113     for (e = 0; e < elements->numElements; e++)
114     {
115     if (elements->Color[e] == color)
116     {
117     for (kr = 0; kr < NN; kr++)
118     {
119     irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
120     if ((firstRow <= irow) && (irow < lastRow))
121     {
122     irow_loc = irow - firstRow;
123     for (kc = 0; kc < NN; kc++)
124     {
125     icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
126     if (icol != irow)
127     Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128     }
129     }
130     }
131     }
132     }
133     }
134 gross 1722 }
135     }
136 ksteube 1312
137 jfenwick 3086 /* inserts row index row into the Dudley_IndexList in if it does not exist */
138 jgs 82
139 jfenwick 3224 void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140     {
141     dim_t i;
142     /* is index in in? */
143     for (i = 0; i < in->n; i++)
144     {
145     if (in->index[i] == index)
146     return;
147     }
148     /* index could not be found */
149     if (in->n == INDEXLIST_LENGTH)
150     {
151     /* if in->index is full check the extension */
152     if (in->extension == NULL)
153     {
154     in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155     if (Dudley_checkPtr(in->extension))
156     return;
157     in->extension->n = 0;
158     in->extension->extension = NULL;
159     }
160     Dudley_IndexList_insertIndex(in->extension, index);
161     }
162     else
163     {
164     /* insert index into in->index */
165     in->index[in->n] = index;
166     in->n++;
167     }
168 jgs 82 }
169    
170 jfenwick 3086 /* counts the number of row indices in the Dudley_IndexList in */
171 jgs 82
172 jfenwick 3224 dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173     {
174     dim_t i;
175     dim_t out = 0;
176     register index_t itmp;
177     if (in == NULL)
178     {
179     return 0;
180 ksteube 1312 }
181 jfenwick 3224 else
182     {
183     for (i = 0; i < in->n; i++)
184     {
185     itmp = in->index[i];
186     if ((itmp >= range_min) && (range_max > itmp))
187     ++out;
188     }
189     return out + Dudley_IndexList_count(in->extension, range_min, range_max);
190     }
191 jgs 82 }
192    
193 jfenwick 3086 /* count the number of row indices in the Dudley_IndexList in */
194 jgs 82
195 jfenwick 3224 void Dudley_IndexList_toArray(Dudley_IndexList * in, index_t * array, index_t range_min, index_t range_max,
196     index_t index_offset)
197     {
198     dim_t i, ptr;
199     register index_t itmp;
200     if (in != NULL)
201     {
202     ptr = 0;
203     for (i = 0; i < in->n; i++)
204     {
205     itmp = in->index[i];
206     if ((itmp >= range_min) && (range_max > itmp))
207     {
208     array[ptr] = itmp + index_offset;
209     ptr++;
210     }
211 ksteube 1312
212 jfenwick 3224 }
213     Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
214 ksteube 1312 }
215 jgs 82 }
216    
217 jfenwick 3086 /* deallocates the Dudley_IndexList in by recursive calls */
218 jgs 82
219 jfenwick 3224 void Dudley_IndexList_free(Dudley_IndexList * in)
220     {
221     if (in != NULL)
222     {
223     Dudley_IndexList_free(in->extension);
224     TMPMEMFREE(in);
225     }
226 jgs 82 }
227    
228 ksteube 1312 /* creates a Paso_pattern from a range of indices */
229 jfenwick 3224 Paso_Pattern *Dudley_IndexList_createPattern(dim_t n0, dim_t n, Dudley_IndexList * index_list, index_t range_min,
230     index_t range_max, index_t index_offset)
231 ksteube 1312 {
232 jfenwick 3224 dim_t *ptr = NULL;
233     register dim_t s, i, itmp;
234     index_t *index = NULL;
235     Paso_Pattern *out = NULL;
236 ksteube 1312
237 jfenwick 3224 ptr = MEMALLOC(n + 1 - n0, index_t);
238     if (!Dudley_checkPtr(ptr))
239     {
240     /* get the number of connections per row */
241     #pragma omp parallel for schedule(static) private(i)
242     for (i = n0; i < n; ++i)
243     {
244     ptr[i - n0] = Dudley_IndexList_count(&index_list[i], range_min, range_max);
245     }
246     /* accumulate ptr */
247     s = 0;
248     for (i = n0; i < n; ++i)
249     {
250     itmp = ptr[i - n0];
251     ptr[i - n0] = s;
252     s += itmp;
253     }
254     ptr[n - n0] = s;
255     /* fill index */
256     index = MEMALLOC(ptr[n - n0], index_t);
257     if (!Dudley_checkPtr(index))
258     {
259     #pragma omp parallel for schedule(static)
260     for (i = n0; i < n; ++i)
261     {
262     Dudley_IndexList_toArray(&index_list[i], &index[ptr[i - n0]], range_min, range_max, index_offset);
263     }
264 gross 3312 out = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n - n0, range_max + index_offset, ptr, index);
265 jfenwick 3224 }
266     }
267     if (!Dudley_noError())
268     {
269     MEMFREE(ptr);
270     MEMFREE(index);
271     Paso_Pattern_free(out);
272     }
273     return out;
274 ksteube 1312 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26