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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1722 - (hide annotations)
Fri Aug 22 04:20:30 2008 UTC (11 years, 2 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 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 1552 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,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