# Contents of /trunk/paso/src/Pattern_reduceBandwidth.c

Revision 1981 - (show annotations)
Thu Nov 6 05:27:33 2008 UTC (10 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 9981 byte(s)
```More warning removal.

```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2008 by University of Queensland 5 * 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 14 15 /**************************************************************/ 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 /* Author: gross@access.edu.au */ 28 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;inumOutput;++i) { 47 k=label[i]; 48 for (iptr=pattern->ptr[i];iptrptr[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.dega2.deg) { 71 return 1; 72 } else if (a1.idxa2.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;inumInput;++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]=max_LevelWidth_abort) return FALSE; 115 for (i=firstVertexInLevel[nlvls-1];iptr[k];jptr[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 dim_t i, initial_bandwidth, N=pattern->numOutput, numLabledVertices, numLevels, max_LevelWidth, min_deg,deg, numVerticesInTree=0, bandwidth; 138 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 Paso_setError(VALUE_ERROR,"Paso_Pattern_reduceBandwidth: pattern needs to be for a square matrix."); 143 } else { 144 /* printf("relabeling of %d DOFs started.\n",N); */ 145 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 if (! ( Paso_checkPtr(degAndIdx) || Paso_checkPtr(oldLabel) || Paso_checkPtr(AssignedLevel) || Paso_checkPtr(VerticesInTree) || Paso_checkPtr(firstVertexInLevel) ) ) { 151 /* get the initial bandwidth */ 152 #pragma omp parallel for private(i) 153 for (i=0;iptr[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= 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]; iptr[k+1]-pattern->ptr[k]; 193 if (deg= 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= 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=initial_bandwidth) { 227 /* printf("initial labeling used.\n"); */ 228 #pragma omp parallel for private(i) 229 for (i=0;i