/[escript]/branches/domexper/dudley/src/IndexList.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/IndexList.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3144 - (show annotations)
Fri Sep 3 00:49:02 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 8388 byte(s)
row_node, col_node (Assemble params), node_selection 
removed
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
17 /* Dudley: Converting an element list into a matrix shape */
18
19 /**************************************************************/
20
21 #include "IndexList.h"
22
23 /* Translate from distributed/local array indices to global indices */
24
25 /**************************************************************/
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 void Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,
31 bool_t reduce_row_order, index_t* row_map,
32 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 */
35 index_t color;
36 Dudley_ReferenceElement*refElement;
37 dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL;
38 if (elements!=NULL)
39 {
40 NN=elements->numNodes;
41 refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
42 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 } else {
57 row_node=refElement->Type->subElementNodes;
58 NN_row=(refElement->BasisFunctions->Type->numShapes) ;
59 }
60
61 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 }
79 }
80 }
81 return;
82 }
83
84
85 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 {
88 /* this does not resolve macro elements */
89 index_t color;
90 dim_t e,kr,kc,icol,irow, NN;
91 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 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
104 }
105 }
106 }
107 }
108 }
109 }
110 }
111 }
112 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 {
115 /* this does not resolve macro elements */
116 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 if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
131 }
132 }
133 }
134 }
135 }
136 }
137 }
138 }
139
140 /* inserts row index row into the Dudley_IndexList in if it does not exist */
141
142 void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
143 dim_t i;
144 /* is index in in? */
145 for (i=0;i<in->n;i++) {
146 if (in->index[i]==index) return;
147 }
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 in->extension=TMPMEMALLOC(1,Dudley_IndexList);
153 if (Dudley_checkPtr(in->extension)) return;
154 in->extension->n=0;
155 in->extension->extension=NULL;
156 }
157 Dudley_IndexList_insertIndex(in->extension,index);
158 } else {
159 /* insert index into in->index*/
160 in->index[in->n]=index;
161 in->n++;
162 }
163 }
164
165 /* counts the number of row indices in the Dudley_IndexList in */
166
167 dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
168 dim_t i;
169 dim_t out=0;
170 register index_t itmp;
171 if (in==NULL) {
172 return 0;
173 } else {
174 for (i=0;i<in->n;i++) {
175 itmp=in->index[i];
176 if ((itmp>=range_min) && (range_max>itmp)) ++out;
177 }
178 return out+Dudley_IndexList_count(in->extension, range_min,range_max);
179 }
180 }
181
182 /* count the number of row indices in the Dudley_IndexList in */
183
184 void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
185 dim_t i, ptr;
186 register index_t itmp;
187 if (in!=NULL) {
188 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 Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
198 }
199 }
200
201 /* deallocates the Dudley_IndexList in by recursive calls */
202
203 void Dudley_IndexList_free(Dudley_IndexList* in) {
204 if (in!=NULL) {
205 Dudley_IndexList_free(in->extension);
206 TMPMEMFREE(in);
207 }
208 }
209
210 /* creates a Paso_pattern from a range of indices */
211 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 {
213 dim_t *ptr=NULL;
214 register dim_t s,i,itmp;
215 index_t *index=NULL;
216 Paso_Pattern* out=NULL;
217
218 ptr=MEMALLOC(n+1-n0,index_t);
219 if (! Dudley_checkPtr(ptr) ) {
220 /* get the number of connections per row */
221 #pragma omp parallel for schedule(static) private(i)
222 for(i=n0;i<n;++i) {
223 ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
224 }
225 /* accumulate ptr */
226 s=0;
227 for(i=n0;i<n;++i) {
228 itmp=ptr[i-n0];
229 ptr[i-n0]=s;
230 s+=itmp;
231 }
232 ptr[n-n0]=s;
233 /* fill index */
234 index=MEMALLOC(ptr[n-n0],index_t);
235 if (! Dudley_checkPtr(index)) {
236 #pragma omp parallel for schedule(static)
237 for(i=n0;i<n;++i) {
238 Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
239 }
240 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
241 }
242 }
243 if (! Dudley_noError()) {
244 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