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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4257 - (hide annotations)
Wed Feb 27 03:42:40 2013 UTC (6 years, 7 months ago) by jfenwick
Original Path: branches/doubleplusgood/dudley/src/IndexList.c
File MIME type: text/plain
File size: 7485 byte(s)
Some simple experiments for c++ Finley

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26