← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/IndexList.c

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/IndexList.c
 
1
 
 
2
/*******************************************************
 
3
*
 
4
* Copyright (c) 2003-2010 by University of Queensland
 
5
* 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
 
 
14
/**************************************************************/
 
15
 
 
16
/* Dudley: Converting an element list into a matrix shape     */
 
17
 
 
18
/**************************************************************/
 
19
 
 
20
#include "IndexList.h"
 
21
 
 
22
/* Translate from distributed/local array indices to global indices */
 
23
 
 
24
/**************************************************************/
 
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
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
{
 
33
    /* 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
 
 
42
        for (color = elements->minColor; color <= elements->maxColor; color++)
 
43
        {
 
44
#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
                {
 
49
                    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
                        {
 
54
                            icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
 
55
                            Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
 
56
                        }
 
57
                    }
 
58
                }
 
59
            }
 
60
        }
 
61
    }
 
62
    return;
 
63
}
 
64
 
 
65
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
{
 
68
/* this does not resolve macro elements */
 
69
    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
    }
 
98
}
 
99
 
 
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
{
 
104
    /* 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
    }
 
135
}
 
136
 
 
137
/* inserts row index row into the Dudley_IndexList in if it does not exist */
 
138
 
 
139
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
}
 
169
 
 
170
/* counts the number of row indices in the Dudley_IndexList in */
 
171
 
 
172
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
    }
 
181
    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
}
 
192
 
 
193
/* count the number of row indices in the Dudley_IndexList in */
 
194
 
 
195
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
 
 
212
        }
 
213
        Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
 
214
    }
 
215
}
 
216
 
 
217
/* deallocates the Dudley_IndexList in by recursive calls */
 
218
 
 
219
void Dudley_IndexList_free(Dudley_IndexList * in)
 
220
{
 
221
    if (in != NULL)
 
222
    {
 
223
        Dudley_IndexList_free(in->extension);
 
224
        TMPMEMFREE(in);
 
225
    }
 
226
}
 
227
 
 
228
/* creates a Paso_pattern from a range of indices */
 
229
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
{
 
232
    dim_t *ptr = NULL;
 
233
    register dim_t s, i, itmp;
 
234
    index_t *index = NULL;
 
235
    Paso_Pattern *out = NULL;
 
236
 
 
237
    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
            out = Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT, n - n0, range_max + index_offset, ptr, index);
 
265
        }
 
266
    }
 
267
    if (!Dudley_noError())
 
268
    {
 
269
        MEMFREE(ptr);
 
270
        MEMFREE(index);
 
271
        Paso_Pattern_free(out);
 
272
    }
 
273
    return out;
 
274
}