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

  ViewVC Help
Powered by ViewVC 1.1.26