/[escript]/trunk/finley/src/IndexList.c
ViewVC logotype

Contents of /trunk/finley/src/IndexList.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1736 - (show annotations)
Fri Aug 29 02:23:16 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/plain
File size: 8493 byte(s)
This fixes a problem which is typically arising when using reduced order
with MPI and a "small" number of elements per processor. In this case it
can happen that the couple matrix is not using all entries sent to the
processor. The old implementations assumed that the indices will cover
the entire input. This assumption has been removed.


1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* Finley: Converting an element list into a matrix shape */
19
20 /**************************************************************/
21
22 #include "IndexList.h"
23
24 /* Translate from distributed/local array indices to global indices */
25
26 /**************************************************************/
27 /* 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 void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
32 bool_t reduce_row_order, index_t* row_map,
33 bool_t reduce_col_order, index_t* col_map) {
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, *id=NULL;
36 dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
37 if (elements!=NULL) {
38 NN=elements->numNodes;
39 id=TMPMEMALLOC(NN, index_t);
40 if (! Finley_checkPtr(id) ) {
41 for (i=0;i<NN;i++) id[i]=i;
42 if (reduce_col_order) {
43 col_node=elements->ReferenceElement->Type->linearNodes;
44 NN_col=elements->LinearReferenceElement->Type->numNodes;
45 } else {
46 col_node=id;
47 NN_col=elements->ReferenceElement->Type->numNodes;
48 }
49 if (reduce_row_order) {
50 row_node=elements->ReferenceElement->Type->linearNodes;
51 NN_row=elements->LinearReferenceElement->Type->numNodes;
52 } else {
53 row_node=id;
54 NN_row=elements->ReferenceElement->Type->numNodes;
55 }
56 for (color=elements->minColor;color<=elements->maxColor;color++) {
57 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
58 for (e=0;e<elements->numElements;e++) {
59 if (elements->Color[e]==color) {
60 for (kr=0;kr<NN_row;kr++) {
61 irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
62 for (kc=0;kc<NN_col;kc++) {
63 icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
64 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
65 }
66 }
67 }
68 }
69 }
70 TMPMEMFREE(id);
71 }
72 }
73 return;
74 }
75
76 void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
77 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
78 {
79 index_t color;
80 dim_t e,kr,kc,icol,irow, NN;
81 if (elements!=NULL) {
82 NN=elements->numNodes;
83 for (color=elements->minColor;color<=elements->maxColor;color++) {
84 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
85 for (e=0;e<elements->numElements;e++) {
86 if (elements->Color[e]==color) {
87 for (kr=0;kr<NN;kr++) {
88 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
89 if ((firstRow<=irow) && (irow < lastRow)) {
90 irow-=firstRow;
91 for (kc=0;kc<NN;kc++) {
92 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
93 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
94 }
95 }
96 }
97 }
98 }
99 }
100 }
101 }
102 void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
103 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
104 {
105 index_t color;
106 dim_t e,kr,kc,icol,irow, NN,irow_loc;
107 if (elements!=NULL) {
108 NN=elements->numNodes;
109 for (color=elements->minColor;color<=elements->maxColor;color++) {
110 #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
111 for (e=0;e<elements->numElements;e++) {
112 if (elements->Color[e]==color) {
113 for (kr=0;kr<NN;kr++) {
114 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
115 if ((firstRow<=irow) && (irow < lastRow)) {
116 irow_loc=irow-firstRow;
117 for (kc=0;kc<NN;kc++) {
118 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
119 if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
120 }
121 }
122 }
123 }
124 }
125 }
126 }
127 }
128
129 /* inserts row index row into the Finley_IndexList in if it does not exist */
130
131 void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
132 dim_t i;
133 /* is index in in? */
134 for (i=0;i<in->n;i++) {
135 if (in->index[i]==index) return;
136 }
137 /* index could not be found */
138 if (in->n==INDEXLIST_LENGTH) {
139 /* if in->index is full check the extension */
140 if (in->extension==NULL) {
141 in->extension=TMPMEMALLOC(1,Finley_IndexList);
142 if (Finley_checkPtr(in->extension)) return;
143 in->extension->n=0;
144 in->extension->extension=NULL;
145 }
146 Finley_IndexList_insertIndex(in->extension,index);
147 } else {
148 /* insert index into in->index*/
149 in->index[in->n]=index;
150 in->n++;
151 }
152 }
153
154 /* counts the number of row indices in the Finley_IndexList in */
155
156 dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
157 dim_t i;
158 dim_t out=0;
159 register index_t itmp;
160 if (in==NULL) {
161 return 0;
162 } else {
163 for (i=0;i<in->n;i++) {
164 itmp=in->index[i];
165 if ((itmp>=range_min) && (range_max>itmp)) ++out;
166 }
167 return out+Finley_IndexList_count(in->extension, range_min,range_max);
168 }
169 }
170
171 /* count the number of row indices in the Finley_IndexList in */
172
173 void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
174 dim_t i, ptr;
175 register index_t itmp;
176 if (in!=NULL) {
177 ptr=0;
178 for (i=0;i<in->n;i++) {
179 itmp=in->index[i];
180 if ((itmp>=range_min) && (range_max>itmp)) {
181 array[ptr]=itmp+index_offset;
182 ptr++;
183 }
184
185 }
186 Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
187 }
188 }
189
190 /* deallocates the Finley_IndexList in by recursive calls */
191
192 void Finley_IndexList_free(Finley_IndexList* in) {
193 if (in!=NULL) {
194 Finley_IndexList_free(in->extension);
195 TMPMEMFREE(in);
196 }
197 }
198
199 /* creates a Paso_pattern from a range of indices */
200 Paso_Pattern* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
201 {
202 dim_t *ptr=NULL;
203 register dim_t s,i,itmp;
204 index_t *index=NULL;
205 Paso_Pattern* out=NULL;
206
207 ptr=MEMALLOC(n+1-n0,index_t);
208 if (! Finley_checkPtr(ptr) ) {
209 /* get the number of connections per row */
210 #pragma omp parallel for schedule(static) private(i)
211 for(i=n0;i<n;++i) {
212 ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
213 }
214 /* accumulate ptr */
215 s=0;
216 for(i=n0;i<n;++i) {
217 itmp=ptr[i-n0];
218 ptr[i-n0]=s;
219 s+=itmp;
220 }
221 ptr[n-n0]=s;
222 /* fill index */
223 index=MEMALLOC(ptr[n-n0],index_t);
224 if (! Finley_checkPtr(index)) {
225 #pragma omp parallel for schedule(static)
226 for(i=n0;i<n;++i) {
227 Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
228 }
229 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,range_max+index_offset,ptr,index);
230 }
231 }
232 if (! Finley_noError()) {
233 MEMFREE(ptr);
234 MEMFREE(index);
235 Paso_Pattern_free(out);
236 }
237 return out;
238 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26