/[escript]/branches/doubleplusgood/paso/src/Pattern_reduceBandwidth.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/Pattern_reduceBandwidth.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 9 months ago) by jfenwick
File size: 10068 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.26