/[escript]/branches/arrexp_2137_win/paso/src/Pattern_reduceBandwidth.c
ViewVC logotype

Contents of /branches/arrexp_2137_win/paso/src/Pattern_reduceBandwidth.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2202 - (show annotations)
Fri Jan 9 01:28:32 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 9981 byte(s)
Branching the array experiments from version 2137.
The idea is to make the changes required for the c++ changes to compile 
on windows without bringing in the later python changes.


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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 /* Paso: pattern: a simple algorithm to create a new labeling */
18 /* of input/output with a minimum bandwidth. The algorithm */
19 /* starts from a vertex with a minumum number of naigbours (connections in the pattern). */
20 /* From this root it searches recursively the naighbours. In the last level */
21 /* a new root is picked (again the vertex with the minumum number of naigbours) and the */
22 /* process is repeated. The algorithm is repeated until the maximum number of vertices */
23 /* in a level cannot be reduced. */
24
25 /**************************************************************/
26
27 /* Author: gross@access.edu.au */
28
29 /**************************************************************/
30
31
32 #include "Paso.h"
33 #include "Pattern.h"
34 #include "Common.h"
35
36 /* calculate initial badwdisth for a given labeling */
37 dim_t Paso_Pattern_getBandwidth(Paso_Pattern* pattern, index_t* label) {
38 register index_t k;
39 index_t iptr;
40 dim_t bandwidth = 0, local_bandwidth=0,i ;
41
42 #pragma omp parallel private(local_bandwidth)
43 {
44 local_bandwidth=0;
45 #pragma omp for private(i,iptr,k)
46 for (i=0;i<pattern->numOutput;++i) {
47 k=label[i];
48 for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1];++iptr) {
49 local_bandwidth = MAX(local_bandwidth, ABS(k - label[pattern->index[iptr]]));
50 }
51 }
52 #pragma omp critical
53 {
54 bandwidth=MAX(local_bandwidth,bandwidth);
55 }
56 }
57 return bandwidth;
58 }
59 /* structures used to create an ordering by increasing degree */
60 typedef struct Paso_DegreeAndIdx {
61 dim_t deg;
62 index_t idx;
63 } Paso_DegreeAndIdx;
64
65 int Paso_comparDegreeAndIdx(const void *arg1,const void *arg2){
66 Paso_DegreeAndIdx a1=*(Paso_DegreeAndIdx*)arg1;
67 Paso_DegreeAndIdx a2=*(Paso_DegreeAndIdx*)arg2;
68 if (a1.deg<a2.deg) {
69 return -1;
70 } else if (a1.deg>a2.deg) {
71 return 1;
72 } else if (a1.idx<a2.idx) {
73 return -1;
74 } else if (a1.idx>a2.idx) {
75 return 1;
76 } else {
77 return 0;
78 }
79 }
80
81 /* Paso_Pattern_dropTree drops a tree in pattern from root */
82
83 /* root - on input the starting point of the tree. */
84 /* AssignedLevel- array of length pattern->numOutput indicating the level assigned to vertex */
85 /* VerticesInTree - on output contains the vertices used in tree (array of length pattern->numOutput) */
86 /* numLevels- on output the number of level used. */
87 /* firstVertexInLevel - on output firstVertexInLevel[i] points to the first vertex of level i in VerticesInTree. */
88 /* firstVertexInLevel[i+1]-firstVertexInLevel[i] is the number of vertices in level i. */
89 /* (array of length pattern->numOutput+1) */
90 /* max_LevelWidth_abort- input param which triggers early return if LevelWidth becomes >= max_LevelWidth_abort */
91 bool_t Paso_Pattern_dropTree(index_t root,
92 Paso_Pattern *pattern,
93 index_t *AssignedLevel,
94 index_t *VerticesInTree,
95 dim_t* numLevels,
96 index_t* firstVertexInLevel,
97 dim_t max_LevelWidth_abort,
98 dim_t N)
99 {
100 dim_t nlvls,i;
101 index_t level_top,k ,itest,j;
102
103 #pragma omp parallel for private(i)
104 for (i=0;i<pattern->numInput;++i) AssignedLevel[i]=-1;
105
106 nlvls=0;
107 AssignedLevel[root] = 0;
108 VerticesInTree[0] = root;
109 firstVertexInLevel[0]=0;
110 level_top=firstVertexInLevel[0]+1;
111 while (firstVertexInLevel[nlvls]<level_top) {
112 nlvls++;
113 firstVertexInLevel[nlvls]=level_top;
114 if (firstVertexInLevel[nlvls]-firstVertexInLevel[nlvls-1]>=max_LevelWidth_abort) return FALSE;
115 for (i=firstVertexInLevel[nlvls-1];i<firstVertexInLevel[nlvls];++i) {
116 k = VerticesInTree[i];
117 for (j=pattern->ptr[k];j<pattern->ptr[k+1];++j) {
118 itest = pattern->index[j];
119 if (AssignedLevel[itest]<0) {
120 #ifdef BOUNDS_CHECK
121 if (itest < 0 || itest >= N) { printf("BOUNDS_CHECK %s %d itest=%d\n", __FILE__, __LINE__, itest); exit(1); }
122 if (level_top < 0 || level_top >= N) { printf("BOUNDS_CHECK %s %d level_top=%d\n", __FILE__, __LINE__, level_top); exit(1); }
123 #endif
124 AssignedLevel[itest] = nlvls;
125 VerticesInTree[level_top] = itest;
126 level_top++;
127 }
128 }
129 }
130 }
131 *numLevels=nlvls;
132 return TRUE;
133 }
134
135 /* the driver */
136 void Paso_Pattern_reduceBandwidth(Paso_Pattern* pattern,index_t* oldToNew) {
137 dim_t i, initial_bandwidth, N=pattern->numOutput, numLabledVertices, numLevels, max_LevelWidth, min_deg,deg, numVerticesInTree=0, bandwidth;
138 Paso_DegreeAndIdx* degAndIdx=NULL;
139 index_t root, *AssignedLevel=NULL, *VerticesInTree=NULL, *firstVertexInLevel=NULL,k, *oldLabel=NULL;
140 /* check input */
141 if (N != pattern->numInput) {
142 Paso_setError(VALUE_ERROR,"Paso_Pattern_reduceBandwidth: pattern needs to be for a square matrix.");
143 } else {
144 /* printf("relabeling of %d DOFs started.\n",N); */
145 degAndIdx=TMPMEMALLOC(N,Paso_DegreeAndIdx);
146 oldLabel=TMPMEMALLOC(N,bool_t);
147 AssignedLevel=TMPMEMALLOC(N,index_t);
148 VerticesInTree=TMPMEMALLOC(N,index_t);
149 firstVertexInLevel=TMPMEMALLOC(N+1,index_t);
150 if (! ( Paso_checkPtr(degAndIdx) || Paso_checkPtr(oldLabel) || Paso_checkPtr(AssignedLevel) || Paso_checkPtr(VerticesInTree) || Paso_checkPtr(firstVertexInLevel) ) ) {
151 /* get the initial bandwidth */
152 #pragma omp parallel for private(i)
153 for (i=0;i<N;++i) oldToNew[i]=i;
154 initial_bandwidth=Paso_Pattern_getBandwidth(pattern,oldToNew);
155 /* printf("initial bandwidth = %d\n",initial_bandwidth); */
156 /* get the initial bandwidth */
157 #pragma omp parallel for private(i)
158 for (i=0;i<N;++i) {
159 oldToNew[i]=-1;
160 degAndIdx[i].idx=i;
161 degAndIdx[i].deg=pattern->ptr[i+1]-pattern->ptr[i];
162 }
163 /* create an ordering with increasing degree */
164 #ifdef USE_QSORTG
165 qsortG(degAndIdx,(size_t)N,sizeof(Paso_DegreeAndIdx),Paso_comparDegreeAndIdx);
166 #else
167 qsort(degAndIdx,(size_t)N,sizeof(Paso_DegreeAndIdx),Paso_comparDegreeAndIdx);
168 #endif
169
170 root=degAndIdx[0].idx;
171 numLabledVertices=0;
172
173 while (root>=0) {
174 max_LevelWidth=N+1;
175
176 while (Paso_Pattern_dropTree(root,pattern,AssignedLevel,VerticesInTree,
177 &numLevels,firstVertexInLevel,max_LevelWidth,N) ) {
178
179 /* find new maximum level width */
180 max_LevelWidth=0;
181 for (i=0;i<numLevels;++i) {
182 #ifdef BOUNDS_CHECK
183 if (i < 0 || i >= N+1) { printf("BOUNDS_CHECK %s %d i=%d N=%d\n", __FILE__, __LINE__, i, N); exit(1); }
184 #endif
185 max_LevelWidth=MAX(max_LevelWidth,firstVertexInLevel[i+1]-firstVertexInLevel[i]);
186 }
187 /* find a vertex in the last level which has minimum degree */
188 min_deg=N+1;
189 root=-1;
190 for (i = firstVertexInLevel[numLevels-1]; i<firstVertexInLevel[numLevels];++i) {
191 k=VerticesInTree[i];
192 deg=pattern->ptr[k+1]-pattern->ptr[k];
193 if (deg<min_deg) {
194 min_deg=deg;
195 root=k;
196 }
197 }
198 /* save the vertices in the current tree */
199 numVerticesInTree=firstVertexInLevel[numLevels];
200 for (i=0;i<firstVertexInLevel[numLevels];++i) {
201 #ifdef BOUNDS_CHECK
202 if (numLabledVertices+i < 0 || numLabledVertices+i >= N) { printf("BOUNDS_CHECK %s %d i=%d numLabledVertices=%d root=%d N=%d firstVertexInLevel[numLevels]=%d\n", __FILE__, __LINE__, i, numLabledVertices, root, N, firstVertexInLevel[numLevels]); exit(1); }
203 #endif
204 oldLabel[numLabledVertices+i]=VerticesInTree[i];
205 }
206 }
207 /* now the vertices in the current tree */
208 for (i=0;i<numVerticesInTree;++i) {
209 #ifdef BOUNDS_CHECK
210 if (numLabledVertices+i < 0 || numLabledVertices+i >= N) { printf("BOUNDS_CHECK %s %d i=%d numLabledVertices=%d root=%d N=%d\n", __FILE__, __LINE__, i, numLabledVertices, root, N); exit(1); }
211 #endif
212 oldToNew[oldLabel[numLabledVertices+i]]=numLabledVertices+i;
213 }
214 numLabledVertices+=numVerticesInTree;
215 /* new search for a vertex which is not labled yet */
216 root=-1;
217 for (i=0;i<N;++i) {
218 if (oldToNew[degAndIdx[i].idx] < 0) {
219 root=degAndIdx[i].idx;
220 break;
221 }
222 }
223 } /* end of while root loop */
224 bandwidth=Paso_Pattern_getBandwidth(pattern,oldToNew);
225 /* printf("bandwidth after DOF relabeling= %d\n",bandwidth); */
226 if (bandwidth>=initial_bandwidth) {
227 /* printf("initial labeling used.\n"); */
228 #pragma omp parallel for private(i)
229 for (i=0;i<N;++i) oldToNew[i]=i;
230 }
231 }
232 TMPMEMFREE(degAndIdx);
233 TMPMEMFREE(oldLabel);
234 TMPMEMFREE(AssignedLevel);
235 TMPMEMFREE(VerticesInTree);
236 TMPMEMFREE(firstVertexInLevel);
237 }
238 }

  ViewVC Help
Powered by ViewVC 1.1.26