/[escript]/trunk/paso/src/Pattern_reduceBandwidth.cpp
ViewVC logotype

Contents of /trunk/paso/src/Pattern_reduceBandwidth.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4803 - (show annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 10 months ago) by caltinay
File size: 10210 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

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

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/Pattern_reduceBandwidth.cpp:3531-3826 /branches/lapack2681/paso/src/Pattern_reduceBandwidth.cpp:2682-2741 /branches/pasowrap/paso/src/Pattern_reduceBandwidth.cpp:3661-3674 /branches/py3_attempt2/paso/src/Pattern_reduceBandwidth.cpp:3871-3891 /branches/restext/paso/src/Pattern_reduceBandwidth.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/Pattern_reduceBandwidth.cpp:3669-3791 /branches/stage3.0/paso/src/Pattern_reduceBandwidth.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/Pattern_reduceBandwidth.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/Pattern_reduceBandwidth.cpp:3517-3974 /release/3.0/paso/src/Pattern_reduceBandwidth.cpp:2591-2601 /trunk/paso/src/Pattern_reduceBandwidth.cpp:4257-4344 /trunk/ripley/test/python/paso/src/Pattern_reduceBandwidth.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26