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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3145 - (show annotations)
Fri Sep 3 01:25:52 2010 UTC (9 years ago) by jfenwick
Original Path: branches/domexper/dudley/src/IndexList.c
File MIME type: text/plain
File size: 8130 byte(s)
Removed linearNodes and subElementNodes from ReferenceElementInfo
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;
38 if (elements!=NULL)
39 {
40 NN=elements->numNodes;
41 refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
42 if (reduce_col_order)
43 {
44 NN_col=(refElement->LinearBasisFunctions->Type->numShapes);
45 }
46 else
47 {
48 NN_col=(refElement->BasisFunctions->Type->numShapes);
49 }
50 if (reduce_row_order)
51 {
52 NN_row=(refElement->LinearBasisFunctions->Type->numShapes);
53 } else {
54 NN_row=(refElement->BasisFunctions->Type->numShapes) ;
55 }
56
57 for (color=elements->minColor;color<=elements->maxColor;color++)
58 {
59 #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
60 for (e=0;e<elements->numElements;e++)
61 {
62 if (elements->Color[e]==color)
63 {
64 for (kr=0;kr<NN_row;kr++)
65 {
66 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
67 for (kc=0;kc<NN_col;kc++)
68 {
69 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
70 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
71 }
72 }
73 }
74 }
75 }
76 }
77 return;
78 }
79
80
81 void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
82 Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
83 {
84 /* this does not resolve macro elements */
85 index_t color;
86 dim_t e,kr,kc,icol,irow, NN;
87 if (elements!=NULL) {
88 NN=elements->numNodes;
89 for (color=elements->minColor;color<=elements->maxColor;color++) {
90 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
91 for (e=0;e<elements->numElements;e++) {
92 if (elements->Color[e]==color) {
93 for (kr=0;kr<NN;kr++) {
94 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
95 if ((firstRow<=irow) && (irow < lastRow)) {
96 irow-=firstRow;
97 for (kc=0;kc<NN;kc++) {
98 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
99 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
100 }
101 }
102 }
103 }
104 }
105 }
106 }
107 }
108 void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
109 Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
110 {
111 /* this does not resolve macro elements */
112 index_t color;
113 dim_t e,kr,kc,icol,irow, NN,irow_loc;
114 if (elements!=NULL) {
115 NN=elements->numNodes;
116 for (color=elements->minColor;color<=elements->maxColor;color++) {
117 #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
118 for (e=0;e<elements->numElements;e++) {
119 if (elements->Color[e]==color) {
120 for (kr=0;kr<NN;kr++) {
121 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
122 if ((firstRow<=irow) && (irow < lastRow)) {
123 irow_loc=irow-firstRow;
124 for (kc=0;kc<NN;kc++) {
125 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
126 if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
127 }
128 }
129 }
130 }
131 }
132 }
133 }
134 }
135
136 /* inserts row index row into the Dudley_IndexList in if it does not exist */
137
138 void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
139 dim_t i;
140 /* is index in in? */
141 for (i=0;i<in->n;i++) {
142 if (in->index[i]==index) return;
143 }
144 /* index could not be found */
145 if (in->n==INDEXLIST_LENGTH) {
146 /* if in->index is full check the extension */
147 if (in->extension==NULL) {
148 in->extension=TMPMEMALLOC(1,Dudley_IndexList);
149 if (Dudley_checkPtr(in->extension)) return;
150 in->extension->n=0;
151 in->extension->extension=NULL;
152 }
153 Dudley_IndexList_insertIndex(in->extension,index);
154 } else {
155 /* insert index into in->index*/
156 in->index[in->n]=index;
157 in->n++;
158 }
159 }
160
161 /* counts the number of row indices in the Dudley_IndexList in */
162
163 dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
164 dim_t i;
165 dim_t out=0;
166 register index_t itmp;
167 if (in==NULL) {
168 return 0;
169 } else {
170 for (i=0;i<in->n;i++) {
171 itmp=in->index[i];
172 if ((itmp>=range_min) && (range_max>itmp)) ++out;
173 }
174 return out+Dudley_IndexList_count(in->extension, range_min,range_max);
175 }
176 }
177
178 /* count the number of row indices in the Dudley_IndexList in */
179
180 void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
181 dim_t i, ptr;
182 register index_t itmp;
183 if (in!=NULL) {
184 ptr=0;
185 for (i=0;i<in->n;i++) {
186 itmp=in->index[i];
187 if ((itmp>=range_min) && (range_max>itmp)) {
188 array[ptr]=itmp+index_offset;
189 ptr++;
190 }
191
192 }
193 Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
194 }
195 }
196
197 /* deallocates the Dudley_IndexList in by recursive calls */
198
199 void Dudley_IndexList_free(Dudley_IndexList* in) {
200 if (in!=NULL) {
201 Dudley_IndexList_free(in->extension);
202 TMPMEMFREE(in);
203 }
204 }
205
206 /* creates a Paso_pattern from a range of indices */
207 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)
208 {
209 dim_t *ptr=NULL;
210 register dim_t s,i,itmp;
211 index_t *index=NULL;
212 Paso_Pattern* out=NULL;
213
214 ptr=MEMALLOC(n+1-n0,index_t);
215 if (! Dudley_checkPtr(ptr) ) {
216 /* get the number of connections per row */
217 #pragma omp parallel for schedule(static) private(i)
218 for(i=n0;i<n;++i) {
219 ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
220 }
221 /* accumulate ptr */
222 s=0;
223 for(i=n0;i<n;++i) {
224 itmp=ptr[i-n0];
225 ptr[i-n0]=s;
226 s+=itmp;
227 }
228 ptr[n-n0]=s;
229 /* fill index */
230 index=MEMALLOC(ptr[n-n0],index_t);
231 if (! Dudley_checkPtr(index)) {
232 #pragma omp parallel for schedule(static)
233 for(i=n0;i<n;++i) {
234 Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
235 }
236 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
237 }
238 }
239 if (! Dudley_noError()) {
240 MEMFREE(ptr);
241 MEMFREE(index);
242 Paso_Pattern_free(out);
243 }
244 return out;
245 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26