/[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 2551 - (show annotations)
Thu Jul 23 09:19:15 2009 UTC (9 years, 11 months ago) by gross
Original Path: trunk/finley/src/IndexList.c
File MIME type: text/plain
File size: 8454 byte(s)
a problem with the sparse matrix unrolling fixed.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 /* Finley: 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 Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
31 bool_t reduce_row_order, index_t* row_map,
32 bool_t reduce_col_order, index_t* col_map) {
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, *id=NULL;
35 dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
36 if (elements!=NULL) {
37 NN=elements->numNodes;
38 id=TMPMEMALLOC(NN, index_t);
39 if (! Finley_checkPtr(id) ) {
40 for (i=0;i<NN;i++) id[i]=i;
41 if (reduce_col_order) {
42 col_node=elements->ReferenceElement->Type->linearNodes;
43 NN_col=elements->LinearReferenceElement->Type->numNodes;
44 } else {
45 col_node=id;
46 NN_col=elements->ReferenceElement->Type->numNodes;
47 }
48 if (reduce_row_order) {
49 row_node=elements->ReferenceElement->Type->linearNodes;
50 NN_row=elements->LinearReferenceElement->Type->numNodes;
51 } else {
52 row_node=id;
53 NN_row=elements->ReferenceElement->Type->numNodes;
54 }
55 for (color=elements->minColor;color<=elements->maxColor;color++) {
56 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
57 for (e=0;e<elements->numElements;e++) {
58 if (elements->Color[e]==color) {
59 for (kr=0;kr<NN_row;kr++) {
60 irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
61 for (kc=0;kc<NN_col;kc++) {
62 icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
63 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
64 }
65 }
66 }
67 }
68 }
69 TMPMEMFREE(id);
70 }
71 }
72 return;
73 }
74
75 void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
76 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
77 {
78 index_t color;
79 dim_t e,kr,kc,icol,irow, NN;
80 if (elements!=NULL) {
81 NN=elements->numNodes;
82 for (color=elements->minColor;color<=elements->maxColor;color++) {
83 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
84 for (e=0;e<elements->numElements;e++) {
85 if (elements->Color[e]==color) {
86 for (kr=0;kr<NN;kr++) {
87 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
88 if ((firstRow<=irow) && (irow < lastRow)) {
89 irow-=firstRow;
90 for (kc=0;kc<NN;kc++) {
91 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
92 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
93 }
94 }
95 }
96 }
97 }
98 }
99 }
100 }
101 void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
102 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
103 {
104 index_t color;
105 dim_t e,kr,kc,icol,irow, NN,irow_loc;
106 if (elements!=NULL) {
107 NN=elements->numNodes;
108 for (color=elements->minColor;color<=elements->maxColor;color++) {
109 #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
110 for (e=0;e<elements->numElements;e++) {
111 if (elements->Color[e]==color) {
112 for (kr=0;kr<NN;kr++) {
113 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
114 if ((firstRow<=irow) && (irow < lastRow)) {
115 irow_loc=irow-firstRow;
116 for (kc=0;kc<NN;kc++) {
117 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
118 if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
119 }
120 }
121 }
122 }
123 }
124 }
125 }
126 }
127
128 /* inserts row index row into the Finley_IndexList in if it does not exist */
129
130 void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
131 dim_t i;
132 /* is index in in? */
133 for (i=0;i<in->n;i++) {
134 if (in->index[i]==index) return;
135 }
136 /* index could not be found */
137 if (in->n==INDEXLIST_LENGTH) {
138 /* if in->index is full check the extension */
139 if (in->extension==NULL) {
140 in->extension=TMPMEMALLOC(1,Finley_IndexList);
141 if (Finley_checkPtr(in->extension)) return;
142 in->extension->n=0;
143 in->extension->extension=NULL;
144 }
145 Finley_IndexList_insertIndex(in->extension,index);
146 } else {
147 /* insert index into in->index*/
148 in->index[in->n]=index;
149 in->n++;
150 }
151 }
152
153 /* counts the number of row indices in the Finley_IndexList in */
154
155 dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
156 dim_t i;
157 dim_t out=0;
158 register index_t itmp;
159 if (in==NULL) {
160 return 0;
161 } else {
162 for (i=0;i<in->n;i++) {
163 itmp=in->index[i];
164 if ((itmp>=range_min) && (range_max>itmp)) ++out;
165 }
166 return out+Finley_IndexList_count(in->extension, range_min,range_max);
167 }
168 }
169
170 /* count the number of row indices in the Finley_IndexList in */
171
172 void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
173 dim_t i, ptr;
174 register index_t itmp;
175 if (in!=NULL) {
176 ptr=0;
177 for (i=0;i<in->n;i++) {
178 itmp=in->index[i];
179 if ((itmp>=range_min) && (range_max>itmp)) {
180 array[ptr]=itmp+index_offset;
181 ptr++;
182 }
183
184 }
185 Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
186 }
187 }
188
189 /* deallocates the Finley_IndexList in by recursive calls */
190
191 void Finley_IndexList_free(Finley_IndexList* in) {
192 if (in!=NULL) {
193 Finley_IndexList_free(in->extension);
194 TMPMEMFREE(in);
195 }
196 }
197
198 /* creates a Paso_pattern from a range of indices */
199 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)
200 {
201 dim_t *ptr=NULL;
202 register dim_t s,i,itmp;
203 index_t *index=NULL;
204 Paso_Pattern* out=NULL;
205
206 ptr=MEMALLOC(n+1-n0,index_t);
207 if (! Finley_checkPtr(ptr) ) {
208 /* get the number of connections per row */
209 #pragma omp parallel for schedule(static) private(i)
210 for(i=n0;i<n;++i) {
211 ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
212 }
213 /* accumulate ptr */
214 s=0;
215 for(i=n0;i<n;++i) {
216 itmp=ptr[i-n0];
217 ptr[i-n0]=s;
218 s+=itmp;
219 }
220 ptr[n-n0]=s;
221 /* fill index */
222 index=MEMALLOC(ptr[n-n0],index_t);
223 if (! Finley_checkPtr(index)) {
224 #pragma omp parallel for schedule(static)
225 for(i=n0;i<n;++i) {
226 Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
227 }
228 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
229 }
230 }
231 if (! Finley_noError()) {
232 MEMFREE(ptr);
233 MEMFREE(index);
234 Paso_Pattern_free(out);
235 }
236 return out;
237 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26