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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1736 - (hide annotations)
Fri Aug 29 02:23:16 2008 UTC (11 years, 7 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 jgs 82
2 ksteube 1312 /* $Id$ */
3 jgs 82
4 ksteube 1312 /*******************************************************
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 jgs 82
16     /**************************************************************/
17    
18 ksteube 1312 /* Finley: Converting an element list into a matrix shape */
19 jgs 82
20     /**************************************************************/
21    
22     #include "IndexList.h"
23    
24 ksteube 1312 /* Translate from distributed/local array indices to global indices */
25    
26 jgs 82 /**************************************************************/
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 ksteube 971 void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
32 ksteube 1312 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 gross 1028 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 jgs 82 if (elements!=NULL) {
38 ksteube 1312 NN=elements->numNodes;
39 gross 1028 id=TMPMEMALLOC(NN, index_t);
40     if (! Finley_checkPtr(id) ) {
41 ksteube 1312 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 jgs 82 }
73     return;
74     }
75    
76 ksteube 1312 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 phornby 1628 dim_t e,kr,kc,icol,irow, NN;
81 ksteube 1312 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 gross 1722 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 ksteube 1312
129 jgs 82 /* inserts row index row into the Finley_IndexList in if it does not exist */
130    
131 jgs 123 void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
132     dim_t i;
133 jgs 82 /* is index in in? */
134     for (i=0;i<in->n;i++) {
135 jgs 102 if (in->index[i]==index) return;
136 jgs 82 }
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 jgs 102 in->extension=TMPMEMALLOC(1,Finley_IndexList);
142 jgs 82 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 ksteube 1312 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 jgs 82 if (in==NULL) {
161     return 0;
162     } else {
163 ksteube 1312 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 jgs 82 }
169     }
170    
171     /* count the number of row indices in the Finley_IndexList in */
172    
173 ksteube 1312 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 jgs 82 if (in!=NULL) {
177 ksteube 1312 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 jgs 82 }
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 ksteube 1312 /* creates a Paso_pattern from a range of indices */
200 gross 1552 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 ksteube 1312 {
202     dim_t *ptr=NULL;
203     register dim_t s,i,itmp;
204     index_t *index=NULL;
205     Paso_Pattern* out=NULL;
206    
207 gross 1552 ptr=MEMALLOC(n+1-n0,index_t);
208 ksteube 1312 if (! Finley_checkPtr(ptr) ) {
209     /* get the number of connections per row */
210     #pragma omp parallel for schedule(static) private(i)
211 gross 1552 for(i=n0;i<n;++i) {
212     ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
213 ksteube 1312 }
214     /* accumulate ptr */
215     s=0;
216 gross 1552 for(i=n0;i<n;++i) {
217     itmp=ptr[i-n0];
218     ptr[i-n0]=s;
219 ksteube 1312 s+=itmp;
220     }
221 gross 1552 ptr[n-n0]=s;
222 ksteube 1312 /* fill index */
223 gross 1552 index=MEMALLOC(ptr[n-n0],index_t);
224 ksteube 1312 if (! Finley_checkPtr(index)) {
225     #pragma omp parallel for schedule(static)
226 gross 1552 for(i=n0;i<n;++i) {
227     Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
228 ksteube 1312 }
229 gross 1736 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,range_max+index_offset,ptr,index);
230 ksteube 1312 }
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