# Contents of /branches/doubleplusgood/paso/src/Pattern_reduceBandwidth.cpp

Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 9 months ago) by jfenwick
File size: 10068 byte(s)

 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2013 by University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development since 2012 by School of Earth Sciences 13 * 14 *****************************************************************************/ 15 16 17 /************************************************************************************/ 18 19 /* Paso: pattern: a simple algorithm to create a new labeling */ 20 /* of input/output with a minimum bandwidth. The algorithm */ 21 /* starts from a vertex with a minimum number of neighbours */ 22 /* (connections in the pattern). From this root it searches */ 23 /* recursively the neighbours. In the last level a new root */ 24 /* is picked (again the vertex with the minimum number of */ 25 /* neighbours) and the process is repeated. The algorithm is */ 26 /* repeated until the maximum number of vertices in a level */ 27 /* cannot be reduced. */ 28 29 /************************************************************************************/ 30 31 /* Author: Lutz Gross, l.gross@uq.edu.au */ 32 33 /************************************************************************************/ 34 35 36 #include "Paso.h" 37 #include "Pattern.h" 38 #include "Common.h" 39 40 /* calculate initial bandwidth for a given labeling */ 41 dim_t Paso_Pattern_getBandwidth(Paso_Pattern* pattern, index_t* label) { 42 register index_t k; 43 index_t iptr; 44 dim_t bandwidth = 0, local_bandwidth=0,i ; 45 46 #pragma omp parallel private(local_bandwidth) 47 { 48 local_bandwidth=0; 49 #pragma omp for private(i,iptr,k) 50 for (i=0;inumOutput;++i) { 51 k=label[i]; 52 for (iptr=pattern->ptr[i];iptrptr[i+1];++iptr) { 53 local_bandwidth = MAX(local_bandwidth, ABS(k - label[pattern->index[iptr]])); 54 } 55 } 56 #pragma omp critical 57 { 58 bandwidth=MAX(local_bandwidth,bandwidth); 59 } 60 } 61 return bandwidth; 62 } 63 /* structures used to create an ordering by increasing degree */ 64 typedef struct Paso_DegreeAndIdx { 65 dim_t deg; 66 index_t idx; 67 } Paso_DegreeAndIdx; 68 69 int Paso_comparDegreeAndIdx(const void *arg1,const void *arg2){ 70 Paso_DegreeAndIdx a1=*(Paso_DegreeAndIdx*)arg1; 71 Paso_DegreeAndIdx a2=*(Paso_DegreeAndIdx*)arg2; 72 if (a1.dega2.deg) { 75 return 1; 76 } else if (a1.idxa2.idx) { 79 return 1; 80 } else { 81 return 0; 82 } 83 } 84 85 /* Paso_Pattern_dropTree drops a tree in pattern from root */ 86 87 /* root - on input the starting point of the tree. */ 88 /* AssignedLevel- array of length pattern->numOutput indicating the level assigned to vertex */ 89 /* VerticesInTree - on output contains the vertices used in tree (array of length pattern->numOutput) */ 90 /* numLevels- on output the number of levels used. */ 91 /* firstVertexInLevel - on output firstVertexInLevel[i] points to the first vertex of level i in VerticesInTree. */ 92 /* firstVertexInLevel[i+1]-firstVertexInLevel[i] is the number of vertices in level i. */ 93 /* (array of length pattern->numOutput+1) */ 94 /* max_LevelWidth_abort- input param which triggers early return if LevelWidth becomes >= max_LevelWidth_abort */ 95 bool_t Paso_Pattern_dropTree(index_t root, 96 Paso_Pattern *pattern, 97 index_t *AssignedLevel, 98 index_t *VerticesInTree, 99 dim_t* numLevels, 100 index_t* firstVertexInLevel, 101 dim_t max_LevelWidth_abort, 102 dim_t N) 103 { 104 dim_t nlvls,i; 105 index_t level_top,k ,itest,j; 106 107 #pragma omp parallel for private(i) 108 for (i=0;inumInput;++i) AssignedLevel[i]=-1; 109 110 nlvls=0; 111 AssignedLevel[root] = 0; 112 VerticesInTree[0] = root; 113 firstVertexInLevel[0]=0; 114 level_top=firstVertexInLevel[0]+1; 115 while (firstVertexInLevel[nlvls]=max_LevelWidth_abort) return FALSE; 119 for (i=firstVertexInLevel[nlvls-1];iptr[k];jptr[k+1];++j) { 122 itest = pattern->index[j]; 123 if (AssignedLevel[itest]<0) { 124 #ifdef BOUNDS_CHECK 125 if (itest < 0 || itest >= N) { printf("BOUNDS_CHECK %s %d itest=%d\n", __FILE__, __LINE__, itest); exit(1); } 126 if (level_top < 0 || level_top >= N) { printf("BOUNDS_CHECK %s %d level_top=%d\n", __FILE__, __LINE__, level_top); exit(1); } 127 #endif 128 AssignedLevel[itest] = nlvls; 129 VerticesInTree[level_top] = itest; 130 level_top++; 131 } 132 } 133 } 134 } 135 *numLevels=nlvls; 136 return TRUE; 137 } 138 139 /* the driver */ 140 void Paso_Pattern_reduceBandwidth(Paso_Pattern* pattern,index_t* oldToNew) { 141 dim_t i, initial_bandwidth, N=pattern->numOutput, numLabledVertices, numLevels, max_LevelWidth, min_deg,deg, numVerticesInTree=0, bandwidth; 142 Paso_DegreeAndIdx* degAndIdx=NULL; 143 index_t root, *AssignedLevel=NULL, *VerticesInTree=NULL, *firstVertexInLevel=NULL,k, *oldLabel=NULL; 144 /* check input */ 145 if (N != pattern->numInput) { 146 Esys_setError(VALUE_ERROR,"Paso_Pattern_reduceBandwidth: pattern needs to be for a square matrix."); 147 } else if (N > 0) { 148 /* printf("relabeling of %d DOFs started.\n",N); */ 149 degAndIdx=new Paso_DegreeAndIdx[N]; 150 oldLabel=new bool_t[N]; 151 AssignedLevel=new index_t[N]; 152 VerticesInTree=new index_t[N]; 153 firstVertexInLevel=new index_t[N+1]; 154 if (! ( Esys_checkPtr(degAndIdx) || Esys_checkPtr(oldLabel) || Esys_checkPtr(AssignedLevel) || Esys_checkPtr(VerticesInTree) || Esys_checkPtr(firstVertexInLevel) ) ) { 155 /* get the initial bandwidth */ 156 #pragma omp parallel for private(i) 157 for (i=0;iptr[i+1]-pattern->ptr[i]; 166 } 167 /* create an ordering with increasing degree */ 168 #ifdef USE_QSORTG 169 qsortG(degAndIdx,(size_t)N,sizeof(Paso_DegreeAndIdx),Paso_comparDegreeAndIdx); 170 #else 171 qsort(degAndIdx,(size_t)N,sizeof(Paso_DegreeAndIdx),Paso_comparDegreeAndIdx); 172 #endif 173 174 root=degAndIdx[0].idx; 175 numLabledVertices=0; 176 177 while (root>=0) { 178 max_LevelWidth=N+1; 179 180 while (Paso_Pattern_dropTree(root,pattern,AssignedLevel,VerticesInTree, 181 &numLevels,firstVertexInLevel,max_LevelWidth,N) ) { 182 183 /* find new maximum level width */ 184 max_LevelWidth=0; 185 for (i=0;i= N+1) { printf("BOUNDS_CHECK %s %d i=%d N=%d\n", __FILE__, __LINE__, i, N); exit(1); } 188 #endif 189 max_LevelWidth=MAX(max_LevelWidth,firstVertexInLevel[i+1]-firstVertexInLevel[i]); 190 } 191 /* find a vertex in the last level which has minimum degree */ 192 min_deg=N+1; 193 root=-1; 194 for (i = firstVertexInLevel[numLevels-1]; iptr[k+1]-pattern->ptr[k]; 197 if (deg= 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); } 207 #endif 208 oldLabel[numLabledVertices+i]=VerticesInTree[i]; 209 } 210 } 211 /* now the vertices in the current tree */ 212 for (i=0;i= N) { printf("BOUNDS_CHECK %s %d i=%d numLabeledVertices=%d root=%d N=%d\n", __FILE__, __LINE__, i, numLabledVertices, root, N); exit(1); } 215 #endif 216 oldToNew[oldLabel[numLabledVertices+i]]=numLabledVertices+i; 217 } 218 numLabledVertices+=numVerticesInTree; 219 /* new search for a vertex which is not labeled yet */ 220 root=-1; 221 for (i=0;i=initial_bandwidth) { 231 /* printf("initial labeling used.\n"); */ 232 #pragma omp parallel for private(i) 233 for (i=0;i