/[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 1722 - (show annotations)
Fri Aug 22 04:20:30 2008 UTC (10 years, 10 months ago) by gross
Original Path: trunk/finley/src/IndexList.c
File MIME type: text/plain
File size: 8470 byte(s)
The CSR matrix patters handed ober to ParMetis may not containn the main diagonal entry. 
This revisison is fixing the problem.


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,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