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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 ksteube 1315
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4657 * Copyright (c) 2003-2014 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 ksteube 1315
17 ksteube 1811
18 jfenwick 3981 /************************************************************************************/
19 ksteube 1315
20 caltinay 3642 /* 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 ksteube 1315
30 jfenwick 3981 /************************************************************************************/
31 ksteube 1315
32 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
33 ksteube 1315
34 jfenwick 3981 /************************************************************************************/
35 ksteube 1315
36    
37     #include "Paso.h"
38     #include "Pattern.h"
39     #include "Common.h"
40    
41 caltinay 4803 namespace paso {
42    
43 caltinay 3642 /* calculate initial bandwidth for a given labeling */
44 caltinay 4803 dim_t Pattern_getBandwidth(Pattern* pattern, index_t* label) {
45 ksteube 1315 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 caltinay 4803 /* Pattern_dropTree drops a tree in pattern from root */
89 ksteube 1315
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 caltinay 3642 /* numLevels- on output the number of levels used. */
94 ksteube 1315 /* 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 caltinay 4803 bool Pattern_dropTree(index_t root,
99     Pattern *pattern,
100 ksteube 1315 index_t *AssignedLevel,
101     index_t *VerticesInTree,
102     dim_t* numLevels,
103     index_t* firstVertexInLevel,
104     dim_t max_LevelWidth_abort,
105 caltinay 4803 dim_t N)
106 ksteube 1315 {
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 caltinay 4803 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 ksteube 1315 #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 caltinay 4803 void Pattern_reduceBandwidth(Pattern* pattern,index_t* oldToNew) {
144 jfenwick 1981 dim_t i, initial_bandwidth, N=pattern->numOutput, numLabledVertices, numLevels, max_LevelWidth, min_deg,deg, numVerticesInTree=0, bandwidth;
145 ksteube 1315 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 caltinay 4803 Esys_setError(VALUE_ERROR,"Pattern_reduceBandwidth: pattern needs to be for a square matrix.");
150 caltinay 2873 } else if (N > 0) {
151 gross 1841 /* printf("relabeling of %d DOFs started.\n",N); */
152 jfenwick 4291 degAndIdx=new Paso_DegreeAndIdx[N];
153 jfenwick 4521 oldLabel=new index_t[N];
154 jfenwick 4291 AssignedLevel=new index_t[N];
155     VerticesInTree=new index_t[N];
156     firstVertexInLevel=new index_t[N+1];
157 jfenwick 3259 if (! ( Esys_checkPtr(degAndIdx) || Esys_checkPtr(oldLabel) || Esys_checkPtr(AssignedLevel) || Esys_checkPtr(VerticesInTree) || Esys_checkPtr(firstVertexInLevel) ) ) {
158 ksteube 1315 /* get the initial bandwidth */
159     #pragma omp parallel for private(i)
160     for (i=0;i<N;++i) oldToNew[i]=i;
161 caltinay 4803 initial_bandwidth=Pattern_getBandwidth(pattern,oldToNew);
162 gross 1841 /* printf("initial bandwidth = %d\n",initial_bandwidth); */
163 ksteube 1315 /* 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 caltinay 4803 while (Pattern_dropTree(root,pattern,AssignedLevel,VerticesInTree,
184 ksteube 1315 &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 caltinay 4803 if (i < 0 || i >= N+1) { printf("BOUNDS_CHECK %s %d i=%d N=%d\n", __FILE__, __LINE__, i, N); exit(1); }
191 ksteube 1315 #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 caltinay 4803 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 ksteube 1315 #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 caltinay 4803 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 ksteube 1315 #endif
219     oldToNew[oldLabel[numLabledVertices+i]]=numLabledVertices+i;
220     }
221     numLabledVertices+=numVerticesInTree;
222 caltinay 3642 /* new search for a vertex which is not labeled yet */
223 ksteube 1315 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 caltinay 4803 bandwidth=Pattern_getBandwidth(pattern,oldToNew);
232 gross 1841 /* printf("bandwidth after DOF relabeling= %d\n",bandwidth); */
233 ksteube 1315 if (bandwidth>=initial_bandwidth) {
234 gross 1841 /* printf("initial labeling used.\n"); */
235 ksteube 1315 #pragma omp parallel for private(i)
236     for (i=0;i<N;++i) oldToNew[i]=i;
237     }
238     }
239 jfenwick 4291 delete[] degAndIdx;
240     delete[] oldLabel;
241     delete[] AssignedLevel;
242     delete[] VerticesInTree;
243     delete[] firstVertexInLevel;
244 ksteube 1315 }
245     }
246 caltinay 3642
247 caltinay 4803 } // 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