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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 10002 byte(s)
Merging dudley and scons updates from branches

1 ksteube 1315
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1315
14 ksteube 1811
15 ksteube 1315 /**************************************************************/
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 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
28 ksteube 1315
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 jfenwick 1981 dim_t i, initial_bandwidth, N=pattern->numOutput, numLabledVertices, numLevels, max_LevelWidth, min_deg,deg, numVerticesInTree=0, bandwidth;
138 ksteube 1315 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 jfenwick 3259 Esys_setError(VALUE_ERROR,"Paso_Pattern_reduceBandwidth: pattern needs to be for a square matrix.");
143 caltinay 2873 } else if (N > 0) {
144 gross 1841 /* printf("relabeling of %d DOFs started.\n",N); */
145 ksteube 1315 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 jfenwick 3259 if (! ( Esys_checkPtr(degAndIdx) || Esys_checkPtr(oldLabel) || Esys_checkPtr(AssignedLevel) || Esys_checkPtr(VerticesInTree) || Esys_checkPtr(firstVertexInLevel) ) ) {
151 ksteube 1315 /* 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 gross 1841 /* printf("initial bandwidth = %d\n",initial_bandwidth); */
156 ksteube 1315 /* 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 gross 1841 /* printf("bandwidth after DOF relabeling= %d\n",bandwidth); */
226 ksteube 1315 if (bandwidth>=initial_bandwidth) {
227 gross 1841 /* printf("initial labeling used.\n"); */
228 ksteube 1315 #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