Line data Source code
1 :
2 : #include "common/common.h"
3 : #include "common/profiler/profiler.h"
4 :
5 : #include <algorithm>
6 : #include <vector>
7 : #include <queue>
8 :
9 : #include "IOrderingAlgorithm.h"
10 : #include "lib_algebra/algebra_common/permutation_util.h"
11 : #include "native_cuthill_mckee.h"
12 :
13 : namespace ug{
14 :
15 0 : bool CheckPermutationBijective(const std::vector<size_t> &perm)
16 : {
17 : std::vector<size_t> invPerm;
18 0 : invPerm.resize(perm.size());
19 : bool bId = true;
20 0 : for(size_t i=0; i<perm.size(); i++) invPerm[i] = (size_t) (-1);
21 :
22 0 : for(size_t i=0; i<perm.size(); i++)
23 : {
24 0 : UG_COND_THROW(invPerm[perm[i]] != (size_t) (-1), "not a bijective permutation "
25 : "(double mapping to index " << perm[i] << " by indices " << invPerm[perm[i]] << " and " << i << ")!");
26 0 : bId = bId && perm[i] == i;
27 0 : invPerm[perm[i]] = i;
28 : }
29 0 : return bId;
30 0 : }
31 :
32 : /// help class to provide compare operator for indices based on their degree
33 : /**
34 : * This class is used to provide an ordering for indices. The ordering relation
35 : * is based on the connectivity-degree, i.e. on the number of connections the
36 : * index has. The indices with less connections are ordered first.
37 : */
38 : struct CompareDegree {
39 : /// constructor, passing field with connections for each index
40 : CompareDegree(const std::vector<std::vector<size_t> >& vInfo) : m_vCon(vInfo) {}
41 :
42 : /// comparison operator
43 : bool operator() (size_t i,size_t j)
44 : {
45 : UG_ASSERT(i < m_vCon.size(), "Invalid index.");
46 : UG_ASSERT(j < m_vCon.size(), "Invalid index.");
47 0 : return (m_vCon[i].size() < m_vCon[j].size());
48 : }
49 :
50 : private:
51 : /// storage field of connections of each index
52 : const std::vector<std::vector<size_t> >& m_vCon;
53 : };
54 :
55 :
56 : /// greatest common divisor
57 : static size_t gcd(size_t a, size_t b)
58 : {
59 0 : while (b)
60 : {
61 0 : size_t r = a % b;
62 : a = b;
63 : b = r;
64 : }
65 : return a;
66 : }
67 :
68 :
69 0 : static size_t findBlockSize(const std::vector<std::vector<size_t> >& vvConnection)
70 : {
71 : size_t nDoF = vvConnection.size();
72 :
73 : // find first non-empty connection
74 : size_t cd = 0;
75 0 : while (cd < nDoF && !vvConnection[cd].size())
76 0 : ++cd;
77 0 : if (cd == nDoF)
78 : return nDoF;
79 :
80 : size_t blockSize = 1;
81 0 : for (size_t i = cd+1; i < nDoF; ++i)
82 : {
83 0 : if (!vvConnection[i].size())
84 : {
85 0 : ++blockSize;
86 0 : continue;
87 : }
88 :
89 : cd = gcd(blockSize, cd);
90 : blockSize = 1;
91 : }
92 :
93 : return gcd(blockSize, cd);
94 : }
95 :
96 :
97 : // computes ordering using Cuthill-McKee algorithm
98 0 : void ComputeCuthillMcKeeOrder(std::vector<size_t>& vNewIndex,
99 : std::vector<std::vector<size_t> >& vvConnection,
100 : bool bReverse,
101 : bool bPreserveConsec)
102 : {
103 : PROFILE_FUNC();
104 : // list of sorted (will be filled later)
105 : std::vector<size_t> vNewOrder;
106 :
107 : const std::size_t nDoF = vvConnection.size();
108 :
109 : // create flag list to remember already handled indices
110 0 : std::vector<bool> vHandled(nDoF, false);
111 :
112 : // Sort neighbours by degree (i.e. by number of neighbours those have)
113 : CompareDegree myCompDegree(vvConnection);
114 0 : for(size_t i = 0; i < nDoF; ++i)
115 : {
116 : // indices with no adjacent indices are marked as handled (and skipped)
117 0 : if(vvConnection[i].size() == 0){
118 : vHandled[i] = true;
119 : }
120 : // sort adjacent index by degree
121 : // edit (mbreit, 10-11-2015): using stable_sort here because implementations
122 : // of std::sort seem to vary among different platforms with regard to the
123 : // sorting outcome of entries where myCompDegree is "=", causing differences
124 : // in convergence behavior in ILU-T
125 : else {
126 0 : std::stable_sort( vvConnection[i].begin(),
127 : vvConnection[i].end(), myCompDegree);
128 : }
129 : }
130 :
131 : // also sort vvConnection itself, this is extremely useful if there are many
132 : // identity rows, as finding the "start" index again and again will take
133 : // REALLY much time in that case
134 0 : std::vector<size_t> sorting(nDoF);
135 : size_t szSort = sorting.size();
136 0 : for (size_t i = 0; i < szSort; ++i)
137 0 : sorting[i] = i;
138 0 : std::stable_sort(sorting.begin(), sorting.end(), myCompDegree);
139 :
140 : // start with first index
141 : size_t firstNonHandled = 0;
142 :
143 : // handle all indices
144 : while(true)
145 : {
146 : // find first unhandled index
147 : size_t i_notHandled = firstNonHandled;
148 0 : for(; i_notHandled < szSort; ++i_notHandled){
149 0 : if(!vHandled[sorting[i_notHandled]]) {firstNonHandled = i_notHandled; break;}
150 : }
151 :
152 : // check if one unhandled vertex left
153 0 : if(i_notHandled == szSort) break;
154 :
155 : /* // This is no longer necessary as firstUnhandled automatically
156 : has smallest degree due to sorting.
157 : // Find node with smallest degree for all remaining indices
158 : size_t start = firstNonHandled;
159 : for(size_t i = start+1; i < vHandled.size(); ++i)
160 : {
161 : if(!vHandled[i] && vvConnection[i].size() < vvConnection[start].size())
162 : start = i;
163 : }
164 : */
165 0 : size_t start = sorting[firstNonHandled];
166 :
167 : // Create queue of adjacent vertices
168 : std::queue<size_t> qAdjacent;
169 : qAdjacent.push(start);
170 :
171 : // add adjacent vertices to mapping
172 0 : while(!qAdjacent.empty())
173 : {
174 : // get next index
175 0 : const size_t front = qAdjacent.front();
176 :
177 : // if not handled
178 0 : if(!vHandled[front])
179 : {
180 : // Add to mapping
181 0 : vNewOrder.push_back(front);
182 : vHandled[front] = true;
183 :
184 : // add adjacent to queue
185 0 : for(size_t i = 0; i < vvConnection[front].size(); ++i)
186 : {
187 0 : const size_t ind = vvConnection[front][i];
188 :
189 0 : if(!vHandled[ind])
190 : qAdjacent.push(ind);
191 : }
192 : }
193 :
194 : // pop index
195 : qAdjacent.pop();
196 : }
197 0 : }
198 :
199 : // Create list of mapping
200 0 : vNewIndex.clear(); vNewIndex.resize(nDoF, (size_t) -1);
201 :
202 0 : if (bPreserveConsec)
203 : {
204 : // write new indices into out array
205 : size_t cnt = 0;
206 0 : if(bReverse)
207 : {
208 0 : for(size_t newInd = 0; newInd < nDoF; ++newInd)
209 : {
210 : // only distribute indices that were already connected before
211 0 : if(vvConnection[newInd].size() == 0) continue;
212 :
213 : // get old index
214 : UG_ASSERT(cnt < vNewOrder.size(), "cnt: " << cnt << ", ordered: " << vNewOrder.size())
215 0 : const size_t oldInd = vNewOrder[vNewOrder.size() - 1 - cnt]; ++cnt;
216 : UG_ASSERT(oldInd < nDoF, "oldInd: "<<oldInd<<", size: "<<nDoF)
217 :
218 : // set new index to order
219 0 : vNewIndex[oldInd] = newInd;
220 : }
221 : }
222 : else
223 : {
224 0 : for(size_t newInd = 0; newInd < nDoF; ++newInd)
225 : {
226 : // skip non-sorted indices
227 0 : if(vvConnection[newInd].size() == 0) continue;
228 :
229 : // get old index
230 : UG_ASSERT(cnt < vNewOrder.size(), "cnt: " << cnt << ", ordered: " << vNewOrder.size())
231 0 : const size_t oldInd = vNewOrder[cnt++];
232 : UG_ASSERT(oldInd < nDoF, "oldInd: "<<oldInd<<", size: "<<nDoF)
233 :
234 : // set new index to order
235 0 : vNewIndex[oldInd] = newInd;
236 : }
237 : }
238 :
239 : // check if all ordered indices have been written
240 0 : if(cnt != vNewOrder.size())
241 0 : UG_THROW("OrderCuthillMcKee: Not all indices sorted that must be sorted: "
242 : << cnt <<" written, but should write: " << vNewOrder.size());
243 :
244 : // fill non-sorted indices (preserving consecutive indexing)
245 : // find minimal block size (most probably the number of functions in underlying algebra)
246 0 : size_t blockSize = findBlockSize(vvConnection);
247 :
248 : // keep blocks, that are due to really unconnected DoFs
249 : // (and not to unblocked algebra), where they were before
250 0 : for (size_t i = 0; i < nDoF; i += blockSize)
251 : {
252 0 : if (vNewIndex[i] == (size_t) -1)
253 0 : vNewIndex[i] = i;
254 : }
255 :
256 : // fill all blocks
257 0 : for (size_t i = 0; i < nDoF; i += blockSize)
258 0 : for (size_t j = 1; j < blockSize; ++j)
259 0 : vNewIndex[i+j] = vNewIndex[i] + j;
260 : }
261 : else
262 : {
263 : size_t newOrdSz = vNewOrder.size();
264 0 : if (bReverse)
265 : {
266 0 : for (size_t i = 0; i < newOrdSz; ++i)
267 : {
268 0 : size_t oldInd = vNewOrder[newOrdSz - 1 - i];
269 0 : vNewIndex[oldInd] = i;
270 : }
271 : }
272 : else
273 : {
274 0 : for (size_t i = 0; i < newOrdSz; ++i)
275 : {
276 0 : size_t oldInd = vNewOrder[i];
277 0 : vNewIndex[oldInd] = i;
278 : }
279 : }
280 :
281 : // move unconnected indices to the bottom of the ordering
282 0 : for (size_t i = 0; i < nDoF; ++i)
283 : {
284 0 : if (vNewIndex[i] == (size_t)-1)
285 0 : vNewIndex[i] = newOrdSz++;
286 : }
287 : }
288 :
289 : #ifdef UG_DEBUG
290 : CheckPermutationBijective(vNewIndex);
291 : #endif
292 0 : }
293 :
294 : }
|
