MFEM  v4.6.0
Finite element discretization library
gecko.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "gecko.hpp"
13 
14 // This file collects the sources of the Gecko library as a single module.
15 // The original library can be found at https://github.com/LLNL/gecko
16 // Used here with permission.
17 
18 // ------------------------------------------------------------------------------
19 
20 // BSD 3-Clause License
21 //
22 // Copyright (c) 2019-2022, Lawrence Livermore National Security, LLC
23 // All rights reserved.
24 //
25 // Redistribution and use in source and binary forms, with or without
26 // modification, are permitted provided that the following conditions are met:
27 //
28 // * Redistributions of source code must retain the above copyright notice, this
29 // list of conditions and the following disclaimer.
30 //
31 // * Redistributions in binary form must reproduce the above copyright notice,
32 // this list of conditions and the following disclaimer in the documentation
33 // and/or other materials provided with the distribution.
34 //
35 // * Neither the name of the copyright holder nor the names of its
36 // contributors may be used to endorse or promote products derived from
37 // this software without specific prior written permission.
38 //
39 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
40 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
41 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
42 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
43 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
44 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
45 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
46 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
47 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
48 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
49 
50 // ------------------------------------------------------------------------------
51 
52 // Copyright (c) 2019-2022, Lawrence Livermore National Security, LLC and other
53 // gecko project contributors. See the above license for details.
54 // SPDX-License-Identifier: BSD-3-Clause
55 // LLNL-CODE-800597
56 
57 #include <algorithm>
58 #include <functional>
59 #include <cstdlib>
60 #include <cstddef>
61 #include <iomanip>
62 #include <iostream>
63 #include <sstream>
64 #include <stdexcept>
65 #include <map>
66 
67 // ----- options.h --------------------------------------------------------------
68 
69 // ratio of max to min weight for aggregation
70 #ifndef GECKO_PART_FRAC
71 #define GECKO_PART_FRAC 4
72 #endif
73 
74 // number of compatible relaxation sweeps
75 #ifndef GECKO_CR_SWEEPS
76 #define GECKO_CR_SWEEPS 1
77 #endif
78 
79 // number of Gauss-Seidel relaxation sweeps
80 #ifndef GECKO_GS_SWEEPS
81 #define GECKO_GS_SWEEPS 1
82 #endif
83 
84 // max number of nodes in subgraph
85 #ifndef GECKO_WINDOW_MAX
86 #define GECKO_WINDOW_MAX 16
87 #endif
88 
89 // use adjacency list (1) or adjacency matrix (0)
90 #ifndef GECKO_WITH_ADJLIST
91 #define GECKO_WITH_ADJLIST 0
92 #endif
93 
94 // use nonrecursive permutation algorithm
95 #ifndef GECKO_WITH_NONRECURSIVE
96 #define GECKO_WITH_NONRECURSIVE 0
97 #endif
98 
99 // use double-precision computations
100 #ifndef GECKO_WITH_DOUBLE_PRECISION
101 #define GECKO_WITH_DOUBLE_PRECISION 0
102 #endif
103 
104 // ----- heap.h -----------------------------------------------------------------
105 
106 template <
107  typename T, // data type
108  typename P, // priority type
109  class C = std::less<P>, // comparator for priorities
110  class M = std::map<T, unsigned int> // maps type T to unsigned integer
111  >
112 class DynamicHeap
113 {
114 public:
115  DynamicHeap(size_t count = 0);
116  ~DynamicHeap() {}
117  void insert(T data, P priority);
118  void update(T data, P priority);
119  bool top(T& data);
120  bool top(T& data, P& priority);
121  bool pop();
122  bool extract(T& data);
123  bool extract(T& data, P& priority);
124  bool erase(T data);
125  bool find(T data) const;
126  bool find(T data, P& priority) const;
127  bool empty() const { return heap.empty(); }
128  size_t size() const { return heap.size(); }
129 private:
130  struct HeapEntry
131  {
132  HeapEntry(P p, T d) : priority(p), data(d) {}
133  P priority;
134  T data;
135  };
136  std::vector<HeapEntry> heap;
137  M index;
138  C lower;
139  void ascend(unsigned int i);
140  void descend(unsigned int i);
141  void swap(unsigned int i, unsigned int j);
142  bool ordered(unsigned int i, unsigned int j) const
143  {
144  return !lower(heap[i].priority, heap[j].priority);
145  }
146  unsigned int parent(unsigned int i) const { return (i - 1) / 2; }
147  unsigned int left(unsigned int i) const { return 2 * i + 1; }
148  unsigned int right(unsigned int i) const { return 2 * i + 2; }
149 };
150 
151 template < typename T, typename P, class C, class M >
152 DynamicHeap<T, P, C, M>::DynamicHeap(size_t count)
153 {
154  heap.reserve(count);
155 }
156 
157 template < typename T, typename P, class C, class M >
158 void
159 DynamicHeap<T, P, C, M>::insert(T data, P priority)
160 {
161  if (index.find(data) != index.end())
162  {
163  update(data, priority);
164  }
165  else
166  {
167  unsigned int i = (unsigned int)heap.size();
168  heap.push_back(HeapEntry(priority, data));
169  ascend(i);
170  }
171 }
172 
173 template < typename T, typename P, class C, class M >
174 void
175 DynamicHeap<T, P, C, M>::update(T data, P priority)
176 {
177  unsigned int i = index[data];
178  heap[i].priority = priority;
179  ascend(i);
180  descend(i);
181 }
182 
183 template < typename T, typename P, class C, class M >
184 bool
185 DynamicHeap<T, P, C, M>::top(T& data)
186 {
187  if (!heap.empty())
188  {
189  data = heap[0].data;
190  return true;
191  }
192  else
193  {
194  return false;
195  }
196 }
197 
198 template < typename T, typename P, class C, class M >
199 bool
200 DynamicHeap<T, P, C, M>::top(T& data, P& priority)
201 {
202  if (!heap.empty())
203  {
204  data = heap[0].data;
205  priority = heap[0].priority;
206  return true;
207  }
208  else
209  {
210  return false;
211  }
212 }
213 
214 template < typename T, typename P, class C, class M >
215 bool
216 DynamicHeap<T, P, C, M>::pop()
217 {
218  if (!heap.empty())
219  {
220  T data = heap[0].data;
221  swap(0, (unsigned int)heap.size() - 1);
222  index.erase(data);
223  heap.pop_back();
224  if (!heap.empty())
225  {
226  descend(0);
227  }
228  return true;
229  }
230  else
231  {
232  return false;
233  }
234 }
235 
236 template < typename T, typename P, class C, class M >
237 bool
238 DynamicHeap<T, P, C, M>::extract(T& data)
239 {
240  if (!heap.empty())
241  {
242  data = heap[0].data;
243  return pop();
244  }
245  else
246  {
247  return false;
248  }
249 }
250 
251 template < typename T, typename P, class C, class M >
252 bool
253 DynamicHeap<T, P, C, M>::extract(T& data, P& priority)
254 {
255  if (!heap.empty())
256  {
257  data = heap[0].data;
258  priority = heap[0].priority;
259  return pop();
260  }
261  else
262  {
263  return false;
264  }
265 }
266 
267 template < typename T, typename P, class C, class M >
268 bool
269 DynamicHeap<T, P, C, M>::erase(T data)
270 {
271  if (index.find(data) == index.end())
272  {
273  return false;
274  }
275  unsigned int i = index[data];
276  swap(i, heap.size() - 1);
277  index.erase(data);
278  heap.pop_back();
279  if (i < heap.size())
280  {
281  ascend(i);
282  descend(i);
283  }
284  return true;
285 }
286 
287 template < typename T, typename P, class C, class M >
288 bool
289 DynamicHeap<T, P, C, M>::find(T data) const
290 {
291  return index.find(data) != index.end();
292 }
293 
294 template < typename T, typename P, class C, class M >
295 bool
296 DynamicHeap<T, P, C, M>::find(T data, P& priority) const
297 {
298  typename M::const_iterator p;
299  if ((p = index.find(data)) == index.end())
300  {
301  return false;
302  }
303  unsigned int i = p->second;
304  priority = heap[i].priority;
305  return true;
306 }
307 
308 template < typename T, typename P, class C, class M >
309 void
310 DynamicHeap<T, P, C, M>::ascend(unsigned int i)
311 {
312  for (unsigned int j; i && !ordered(j = parent(i), i); i = j)
313  {
314  swap(i, j);
315  }
316  index[heap[i].data] = i;
317 }
318 
319 template < typename T, typename P, class C, class M >
320 void
321 DynamicHeap<T, P, C, M>::descend(unsigned int i)
322 {
323  for (unsigned int j, k;
324  (j = ((k = left(i)) < heap.size() && !ordered(i, k) ? k : i),
325  j = ((k = right(i)) < heap.size() && !ordered(j, k) ? k : j)) != i;
326  i = j)
327  {
328  swap(i, j);
329  }
330  index[heap[i].data] = i;
331 }
332 
333 template < typename T, typename P, class C, class M >
334 void
335 DynamicHeap<T, P, C, M>::swap(unsigned int i, unsigned int j)
336 {
337  std::swap(heap[i], heap[j]);
338  index[heap[i].data] = i;
339 }
340 
341 // ----- subgraph.h -------------------------------------------------------------
342 
343 namespace Gecko
344 {
345 
346 // Node in a subgraph.
347 class Subnode
348 {
349 public:
350  typedef unsigned char Index;
351  Float pos; // node position
352  WeightedSum cost; // external cost at this position
353 };
354 
355 class Subgraph
356 {
357 public:
358  Subgraph(Graph* g, uint n);
359  ~Subgraph() { delete[] cache; }
360  void optimize(uint k);
361 
362 private:
363  Graph* const g; // full graph
364  const uint n; // number of subgraph nodes
365  Functional* const f; // ordering functional
366  WeightedSum min; // minimum cost so far
367  Subnode::Index best[GECKO_WINDOW_MAX]; // best permutation so far
368  Subnode::Index perm[GECKO_WINDOW_MAX]; // current permutation
369  const Subnode* node[GECKO_WINDOW_MAX]; // pointers to precomputed nodes
370  Subnode* cache; // precomputed node positions and costs
371 #if GECKO_WITH_ADJLIST
372  Subnode::Index
373  adj[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal adjacency list
374 #else
375  uint adj[GECKO_WINDOW_MAX]; // internal adjacency matrix
376 #endif
377  Float weight[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal arc weights
378  WeightedSum cost(uint k) const;
379  void swap(uint k);
380  void swap(uint k, uint l);
381  void optimize(WeightedSum c, uint i);
382 };
383 
384 }
385 
386 // ----- subgraph.cpp -----------------------------------------------------------
387 
388 using namespace Gecko;
389 
390 // Constructor.
391 Subgraph::Subgraph(Graph* g, uint n) : g(g), n(n), f(g->functional)
392 {
393  if (n > GECKO_WINDOW_MAX)
394  {
395  throw std::out_of_range("optimization window too large");
396  }
397  cache = new Subnode[n << n];
398 }
399 
400 // Cost of k'th node's edges to external nodes and nodes at {k+1, ..., n-1}.
402 Subgraph::cost(uint k) const
403 {
404  Subnode::Index i = perm[k];
405  WeightedSum c = node[i]->cost;
406  Float p = node[i]->pos;
407 #if GECKO_WITH_ADJLIST
408  for (k = 0; adj[i][k] != i; k++)
409  {
410  Subnode::Index j = adj[i][k];
411  Float l = node[j]->pos - p;
412  if (l > 0)
413  {
414  Float w = weight[i][k];
415  f->accumulate(c, WeightedValue(l, w));
416  }
417  }
418 #else
419  uint m = adj[i];
420  while (++k < n)
421  {
422  Subnode::Index j = perm[k];
423  if (m & (1u << j))
424  {
425  Float l = node[j]->pos - p;
426  Float w = weight[i][j];
427  f->accumulate(c, WeightedValue(l, w));
428  }
429  }
430 #endif
431  return c;
432 }
433 
434 // Swap the two nodes in positions k and k + 1.
435 void
436 Subgraph::swap(uint k)
437 {
438  uint l = k + 1;
439  Subnode::Index i = perm[k];
440  Subnode::Index j = perm[l];
441  perm[k] = j;
442  perm[l] = i;
443  node[i] -= ptrdiff_t(1) << j;
444  node[j] += ptrdiff_t(1) << i;
445 }
446 
447 // Swap the two nodes in positions k and l, k <= l.
448 void
449 Subgraph::swap(uint k, uint l)
450 {
451  Subnode::Index i = perm[k];
452  Subnode::Index j = perm[l];
453  perm[k] = j;
454  perm[l] = i;
455  // Update node positions.
456  uint m = 0;
457  while (++k < l)
458  {
459  Subnode::Index h = perm[k];
460  node[h] += ptrdiff_t(1) << i;
461  node[h] -= ptrdiff_t(1) << j;
462  m += 1u << h;
463  }
464  node[i] -= (1u << j) + m;
465  node[j] += (1u << i) + m;
466 }
467 
468 #if GECKO_WITH_NONRECURSIVE
469 // Evaluate all permutations generated by Heap's nonrecursive algorithm.
470 void
471 Subgraph::optimize(WeightedSum, uint)
472 {
473  WeightedSum c[GECKO_WINDOW_MAX + 1];
474  uint j[GECKO_WINDOW_MAX + 1];
475  j[n] = 1;
476  c[n] = 0;
477  uint i = n;
478  do
479  {
480  i--;
481  j[i] = i;
482  loop:
483  c[i] = f->sum(c[i + 1], cost(i));
484  }
485  while (i);
486  if (f->less(c[0], min))
487  {
488  min = c[0];
489  for (uint k = 0; k < n; k++)
490  {
491  best[k] = perm[k];
492  }
493  }
494  do
495  {
496  if (++i == n)
497  {
498  return;
499  }
500  swap(i & 1 ? i - j[i] : 0, i);
501  }
502  while (!j[i]--);
503  goto loop;
504 }
505 #else
506 // Apply branch-and-bound to permutations generated by Heap's algorithm.
507 void
508 Subgraph::optimize(WeightedSum c, uint i)
509 {
510  i--;
511  if (f->less(c, min))
512  {
513  if (i)
514  {
515  uint j = i;
516  do
517  {
518  optimize(f->sum(c, cost(i)), i);
519  swap(i & 1 ? i - j : 0, i);
520  }
521  while (j--);
522  }
523  else
524  {
525  f->accumulate(c, cost(0));
526  if (f->less(c, min))
527  {
528  min = c;
529  for (uint j = 0; j < n; j++)
530  {
531  best[j] = perm[j];
532  }
533  }
534  }
535  }
536  else if (i & 1)
537  do { swap(--i); }
538  while (i);
539 }
540 #endif
541 
542 // Optimize layout of nodes {p, ..., p + n - 1}.
543 void
544 Subgraph::optimize(uint p)
545 {
546  // Initialize subgraph.
547  const Float q = g->node[g->perm[p]].pos - g->node[g->perm[p]].hlen;
548  min = WeightedSum(GECKO_FLOAT_MAX, 1);
549  // Subnode::Index's "k" is unsigned char, "n" is unsigned int
550  for (Subnode::Index k = 0; (uint) k < n; k++)
551  {
552  best[k] = perm[k] = k;
553  Node::Index i = g->perm[p + k];
554  // Copy i's outgoing arcs. We distinguish between internal
555  // and external arcs to nodes within and outside the subgraph,
556  // respectively.
557 #if GECKO_WITH_ADJLIST
558  uint m = 0;
559 #else
560  adj[k] = 0;
561 #endif
562  std::vector<Arc::Index> external;
563  for (Arc::Index a = g->node_begin(i); a < g->node_end(i); a++)
564  {
565  Node::Index j = g->adj[a];
566  Subnode::Index l;
567  for (l = 0; (uint) l < n && g->perm[p + l] != j; l++);
568  if (l == n)
569  {
570  external.push_back(a);
571  }
572  else
573  {
574  // Copy internal arc to subgraph.
575 #if GECKO_WITH_ADJLIST
576  adj[k][m] = l;
577  weight[k][m] = g->weight[a];
578  m++;
579 #else
580  adj[k] += 1u << l;
581  weight[k][l] = g->weight[a];
582 #endif
583  }
584  }
585 #if GECKO_WITH_ADJLIST
586  adj[k][m] = k;
587 #endif
588  // Precompute external costs associated with all possible positions
589  // of this node. Since node lengths can be arbitrary, there are as
590  // many as 2^(n-1) possible positions, each corresponding to an
591  // (n-1)-bit string that specifies whether the remaining n-1 nodes
592  // succeed this node or not. Caching the
593  // n
594  // 2^(n-1) n = sum k C(n, k) = A001787
595  // k=1
596  // external costs is exponentially cheaper than recomputing the
597  // n-1 n
598  // n! sum 1/k! = sum k! C(n, k) = A007526
599  // k=0 k=1
600  // costs associated with all permutations.
601  node[k] = cache + (k << n);
602  for (uint m = 0; m < (1u << n); m++)
603  if (!(m & (1u << k)))
604  {
605  Subnode* s = cache + (k << n) + m;
606  s->pos = q + g->node[i].hlen;
607  for (Subnode::Index l = 0; (uint) l < n; l++)
608  if (l != k && !(m & (1u << l)))
609  {
610  s->pos += 2 * g->node[g->perm[p + l]].hlen;
611  }
612  s->cost = g->cost(external, s->pos);
613  }
614  else
615  {
616  m += (1u << k) - 1;
617  }
618  node[k] += (1u << n) - (2u << k);
619  }
620 
621  // Find optimal permutation of the n nodes.
622  optimize(0, n);
623 
624  // Apply permutation to original graph.
625  for (uint i = 0; i < n; i++)
626  {
627  g->swap(p + i, p + best[i]);
628  for (uint j = i + 1; j < n; j++)
629  if (best[j] == i)
630  {
631  best[j] = best[i];
632  }
633  }
634 }
635 
636 // ----- graph.cpp --------------------------------------------------------------
637 
638 using namespace std;
639 using namespace Gecko;
640 
641 // Constructor.
642 void
643 Graph::init(uint nodes)
644 {
645  node.push_back(Node(-1, 0, 1, Node::null));
646  adj.push_back(Node::null);
647  weight.push_back(0);
648  bond.push_back(0);
649  while (nodes--)
650  {
651  insert_node();
652  }
653 }
654 
655 // Insert node.
657 Graph::insert_node(Float length)
658 {
659  Node::Index p = Node::Index(node.size());
660  perm.push_back(p);
661  node.push_back(Node(-1, length));
662  return p;
663 }
664 
665 // Return nodes adjacent to i.
666 std::vector<Node::Index>
667 Graph::node_neighbors(Node::Index i) const
668 {
669  std::vector<Node::Index> neighbor;
670  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
671  {
672  neighbor.push_back(adj[a]);
673  }
674  return neighbor;
675 }
676 
677 // Insert directed edge (i, j).
679 Graph::insert_arc(Node::Index i, Node::Index j, Float w, Float b)
680 {
681  if (!i || !j || i == j || !(last_node <= i && i <= nodes()))
682  {
683  return Arc::null;
684  }
685  last_node = i;
686  for (Node::Index k = i - 1; node[k].arc == Arc::null; k--)
687  {
688  node[k].arc = Arc::Index(adj.size());
689  }
690  adj.push_back(j);
691  weight.push_back(w);
692  bond.push_back(b);
693  node[i].arc = Arc::Index(adj.size());
694  return Arc::Index(adj.size() - 1);
695 }
696 
697 // Remove arc a.
698 bool
699 Graph::remove_arc(Arc::Index a)
700 {
701  if (a == Arc::null)
702  {
703  return false;
704  }
705  Node::Index i = arc_source(a);
706  adj.erase(adj.begin() + a);
707  weight.erase(weight.begin() + a);
708  bond.erase(bond.begin() + a);
709  for (Node::Index k = i; k < node.size(); k++)
710  {
711  node[k].arc--;
712  }
713  return true;
714 }
715 
716 // Remove directed edge (i, j).
717 bool
718 Graph::remove_arc(Node::Index i, Node::Index j)
719 {
720  return remove_arc(arc_index(i, j));
721 }
722 
723 // Remove edge {i, j}.
724 bool
725 Graph::remove_edge(Node::Index i, Node::Index j)
726 {
727  bool success = remove_arc(i, j);
728  if (success)
729  {
730  success = remove_arc(j, i);
731  }
732  return success;
733 }
734 
735 // Index of arc (i, j) or null if not a valid arc.
737 Graph::arc_index(Node::Index i, Node::Index j) const
738 {
739  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
740  if (adj[a] == j)
741  {
742  return a;
743  }
744  return Arc::null;
745 }
746 
747 // Return source node i in arc a = (i, j).
749 Graph::arc_source(Arc::Index a) const
750 {
751  Node::Index j = adj[a];
752  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
753  {
754  Node::Index i = adj[b];
755  if (node_begin(i) <= a && a < node_end(i))
756  {
757  return i;
758  }
759  }
760  // should never get here
761  throw std::runtime_error("internal data structure corrupted");
762 }
763 
764 // Return reverse arc (j, i) of arc a = (i, j).
766 Graph::reverse_arc(Arc::Index a) const
767 {
768  Node::Index j = adj[a];
769  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
770  {
771  Node::Index i = adj[b];
772  if (node_begin(i) <= a && a < node_end(i))
773  {
774  return b;
775  }
776  }
777  return Arc::null;
778 }
779 
780 // Return first directed arc if one exists or null otherwise.
782 Graph::directed() const
783 {
784  for (Node::Index i = 1; i < node.size(); i++)
785  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
786  {
787  Node::Index j = adj[a];
788  if (!arc_index(j, i))
789  {
790  return a;
791  }
792  }
793  return Arc::null;
794 }
795 
796 // Add contribution of fine arc to coarse graph.
797 void
798 Graph::update(Node::Index i, Node::Index j, Float w, Float b)
799 {
800  Arc::Index a = arc_index(i, j);
801  if (a == Arc::null)
802  {
803  insert_arc(i, j, w, b);
804  }
805  else
806  {
807  weight[a] += w;
808  bond[a] += b;
809  }
810 }
811 
812 // Transfer contribution of fine arc a to coarse node p.
813 void
814 Graph::transfer(Graph* g, const vector<Float>& part, Node::Index p,
815  Arc::Index a, Float f) const
816 {
817  Float w = f * weight[a];
818  Float m = f * bond[a];
819  Node::Index j = arc_target(a);
820  Node::Index q = node[j].parent;
821  if (q == Node::null)
822  {
823  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
824  if (part[b] > 0)
825  {
826  q = node[adj[b]].parent;
827  if (q != p)
828  {
829  g->update(p, q, w * part[b], m * part[b]);
830  }
831  }
832  }
833  else
834  {
835  g->update(p, q, w, m);
836  }
837 }
838 
839 // Compute cost of a subset of arcs incident on node placed at pos.
841 Graph::cost(const vector<Arc::Index>& subset, Float pos) const
842 {
843  WeightedSum c;
844  for (Arc::ConstPtr ap = subset.begin(); ap != subset.end(); ap++)
845  {
846  Arc::Index a = *ap;
847  Node::Index j = arc_target(a);
848  Float l = fabs(node[j].pos - pos);
849  Float w = weight[a];
850  functional->accumulate(c, WeightedValue(l, w));
851  }
852  return c;
853 }
854 
855 // Compute cost of graph layout.
856 Float
857 Graph::cost() const
858 {
859  if (edges())
860  {
861  WeightedSum c;
862  Node::Index i = 1;
863  for (Arc::Index a = 1; a < adj.size(); a++)
864  {
865  while (node_end(i) <= a)
866  {
867  i++;
868  }
869  Node::Index j = arc_target(a);
870  Float l = length(i, j);
871  Float w = weight[a];
872  functional->accumulate(c, WeightedValue(l, w));
873  }
874  return functional->mean(c);
875  }
876  else
877  {
878  return Float(0);
879  }
880 }
881 
882 // Swap the two nodes in positions k and l, k <= l.
883 void
884 Graph::swap(uint k, uint l)
885 {
886  Node::Index i = perm[k];
887  perm[k] = perm[l];
888  perm[l] = i;
889  Float p = node[i].pos - node[i].hlen;
890  do
891  {
892  i = perm[k];
893  p += node[i].hlen;
894  node[i].pos = p;
895  p += node[i].hlen;
896  }
897  while (k++ != l);
898 }
899 
900 // Optimize continuous position of a single node.
901 Float
902 Graph::optimal(Node::Index i) const
903 {
904  vector<WeightedValue> v;
905  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
906  {
907  Node::Index j = adj[a];
908  if (placed(j))
909  {
910  v.push_back(WeightedValue(node[j].pos, weight[a]));
911  }
912  }
913  return v.empty() ? -1 : functional->optimum(v);
914 }
915 
916 // Compute coarse graph with roughly half the number of nodes.
917 Graph*
918 Graph::coarsen()
919 {
920  progress->beginphase(this, string("coarse"));
921  Graph* g = new Graph(0, level - 1);
922  g->functional = functional;
923  g->progress = progress;
924 
925  // Compute importance of nodes in fine graph.
926  DynamicHeap<Node::Index, Float> heap;
927  for (Node::Index i = 1; i < node.size(); i++)
928  {
929  node[i].parent = Node::null;
930  Float w = 0;
931  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
932  {
933  w += bond[a];
934  }
935  heap.insert(i, w);
936  }
937 
938  // Select set of important nodes from fine graph that will remain in
939  // coarse graph.
940  vector<Node::Index> child(1, Node::null);
941  while (!heap.empty())
942  {
943  Node::Index i;
944  Float w = 0;
945  heap.extract(i, w);
946  if (w < 0)
947  {
948  break;
949  }
950  child.push_back(i);
951  node[i].parent = g->insert_node(2 * node[i].hlen);
952 
953  // Reduce importance of neighbors.
954  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
955  {
956  Node::Index j = adj[a];
957  if (heap.find(j, w))
958  {
959  heap.update(j, w - 2 * bond[a]);
960  }
961  }
962  }
963 
964  // Assign parts of remaining nodes to aggregates.
965  vector<Float> part = bond;
966  for (Node::Index i = 1; i < node.size(); i++)
967  if (!persistent(i))
968  {
969  // Find all connections to coarse nodes.
970  Float w = 0;
971  Float max = 0;
972  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
973  {
974  Node::Index j = adj[a];
975  if (persistent(j))
976  {
977  w += part[a];
978  if (max < part[a])
979  {
980  max = part[a];
981  }
982  }
983  else
984  {
985  part[a] = -1;
986  }
987  }
988  max /= GECKO_PART_FRAC;
989 
990  // Weed out insignificant connections.
991  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
992  if (0 < part[a] && part[a] < max)
993  {
994  w -= part[a];
995  part[a] = -1;
996  }
997 
998  // Compute node fractions (interpolation matrix) and assign
999  // partial nodes to aggregates.
1000  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1001  if (part[a] > 0)
1002  {
1003  part[a] /= w;
1004  Node::Index p = node[adj[a]].parent;
1005  g->node[p].hlen += part[a] * node[i].hlen;
1006  }
1007  }
1008 
1009  // Transfer arcs to coarse graph.
1010  for (Node::Index p = 1; p < g->node.size(); p++)
1011  {
1012  Node::Index i = child[p];
1013  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1014  {
1015  transfer(g, part, p, a);
1016  Node::Index j = adj[a];
1017  if (!persistent(j))
1018  {
1019  Arc::Index b = arc_index(j, i);
1020  if (part[b] > 0)
1021  for (Arc::Index c = node_begin(j); c < node_end(j); c++)
1022  {
1023  Node::Index k = adj[c];
1024  if (k != i)
1025  {
1026  transfer(g, part, p, c, part[b]);
1027  }
1028  }
1029  }
1030  }
1031  }
1032 
1033 #if DEBUG
1034  if (g->directed())
1035  {
1036  throw runtime_error("directed edge found");
1037  }
1038 #endif
1039 
1040  // Free memory.
1041  vector<Float> t = bond;
1042  bond.swap(t);
1043 
1044  progress->endphase(this, false);
1045 
1046  return g;
1047 }
1048 
1049 // Order nodes according to coarsened graph layout.
1050 void
1051 Graph::refine(const Graph* graph)
1052 {
1053  progress->beginphase(this, string("refine"));
1054 
1055  // Place persistent nodes.
1056  DynamicHeap<Node::Index, Float> heap;
1057  for (Node::Index i = 1; i < node.size(); i++)
1058  if (persistent(i))
1059  {
1060  Node::Index p = node[i].parent;
1061  node[i].pos = graph->node[p].pos;
1062  }
1063  else
1064  {
1065  node[i].pos = -1;
1066  Float w = 0;
1067  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1068  {
1069  Node::Index j = adj[a];
1070  if (persistent(j))
1071  {
1072  w += weight[a];
1073  }
1074  }
1075  heap.insert(i, w);
1076  }
1077 
1078  // Place remaining nodes in order of decreasing connectivity with
1079  // already placed nodes.
1080  while (!heap.empty())
1081  {
1082  Node::Index i = 0;
1083  heap.extract(i);
1084  node[i].pos = optimal(i);
1085  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1086  {
1087  Node::Index j = adj[a];
1088  Float w;
1089  if (heap.find(j, w))
1090  {
1091  heap.update(j, w + weight[a]);
1092  }
1093  }
1094  }
1095 
1096  place(true);
1097  progress->endphase(this, true);
1098 }
1099 
1100 // Perform m sweeps of compatible or Gauss-Seidel relaxation.
1101 void
1102 Graph::relax(bool compatible, uint m)
1103 {
1104  progress->beginphase(this, compatible ? string("crelax") : string("frelax"));
1105  while (m--)
1106  for (uint k = 0; k < perm.size() && !progress->quit(); k++)
1107  {
1108  Node::Index i = perm[k];
1109  if (!compatible || !persistent(i))
1110  {
1111  node[i].pos = optimal(i);
1112  }
1113  }
1114  place(true);
1115  progress->endphase(this, true);
1116 }
1117 
1118 // Optimize successive n-node subgraphs.
1119 void
1120 Graph::optimize(uint n)
1121 {
1122  if (n > perm.size())
1123  {
1124  n = uint(perm.size());
1125  }
1126  ostringstream count;
1127  count << setw(2) << n;
1128  progress->beginphase(this, string("perm") + count.str());
1129  Subgraph* subgraph = new Subgraph(this, n);
1130  for (uint k = 0; k <= perm.size() - n && !progress->quit(); k++)
1131  {
1132  subgraph->optimize(k);
1133  }
1134  delete subgraph;
1135  progress->endphase(this, true);
1136 }
1137 
1138 // Place all nodes according to their positions.
1139 void
1140 Graph::place(bool sort)
1141 {
1142  place(sort, 0, uint(perm.size()));
1143 }
1144 
1145 // Place nodes {k, ..., k + n - 1} according to their positions.
1146 void
1147 Graph::place(bool sort, uint k, uint n)
1148 {
1149  // Place nodes.
1150  if (sort)
1151  {
1152  stable_sort(perm.begin() + k, perm.begin() + k + n,
1153  Node::Comparator(node.begin()));
1154  }
1155 
1156  // Assign node positions according to permutation.
1157  for (Float p = k ? node[perm[k - 1]].pos + node[perm[k - 1]].hlen : 0; n--;
1158  k++)
1159  {
1160  Node::Index i = perm[k];
1161  p += node[i].hlen;
1162  node[i].pos = p;
1163  p += node[i].hlen;
1164  }
1165 }
1166 
1167 // Perform one V-cycle.
1168 void
1169 Graph::vcycle(uint n, uint work)
1170 {
1171  if (n < nodes() && nodes() < edges() && level && !progress->quit())
1172  {
1173  Graph* graph = coarsen();
1174  graph->vcycle(n, work + edges());
1175  refine(graph);
1176  delete graph;
1177  }
1178  else
1179  {
1180  place();
1181  }
1182  if (edges())
1183  {
1184  relax(true, GECKO_CR_SWEEPS);
1185  relax(false, GECKO_GS_SWEEPS);
1186  for (uint w = edges(); w * (n + 1) < work; w *= ++n);
1187  n = std::min(n, uint(GECKO_WINDOW_MAX));
1188  if (n)
1189  {
1190  optimize(n);
1191  }
1192  }
1193 }
1194 
1195 // Custom random-number generator for reproducibility.
1196 // LCG from doi:10.1090/S0025-5718-99-00996-5.
1197 uint
1198 Graph::random(uint seed)
1199 {
1200  static uint state = 1;
1201  state = (seed ? seed : 0x1ed0675 * state + 0xa14f);
1202  return state;
1203 }
1204 
1205 // Generate a random permutation of the nodes.
1206 void
1207 Graph::shuffle(uint seed)
1208 {
1209  random(seed);
1210  for (uint k = 0; k < perm.size(); k++)
1211  {
1212  uint r = random() >> 8;
1213  uint l = k + r % (uint(perm.size()) - k);
1214  std::swap(perm[k], perm[l]);
1215  }
1216  place();
1217 }
1218 
1219 // Recompute bonds for k'th V-cycle.
1220 void
1221 Graph::reweight(uint k)
1222 {
1223  bond.resize(weight.size());
1224  for (Arc::Index a = 1; a < adj.size(); a++)
1225  {
1226  bond[a] = functional->bond(weight[a], length(a), k);
1227  }
1228 }
1229 
1230 // Linearly order graph.
1231 void
1232 Graph::order(Functional* functional_, uint iterations, uint window, uint period,
1233  uint seed, Progress* progress_)
1234 {
1235  // Initialize graph.
1236  this->functional = functional_;
1237  progress_ = this->progress = progress_ ? progress_ : new Progress;
1238  for (level = 0; (1u << level) < nodes(); level++);
1239  place();
1240  Float mincost = cost();
1241  vector<Node::Index> minperm = perm;
1242  if (seed)
1243  {
1244  shuffle(seed);
1245  }
1246 
1247  progress->beginorder(this, mincost);
1248  if (edges())
1249  {
1250  // Perform specified number of V-cycles.
1251  for (uint k = 1; k <= iterations && !progress->quit(); k++)
1252  {
1253  progress->beginiter(this, k, iterations, window);
1254  reweight(k);
1255  vcycle(window);
1256  Float c = cost();
1257  if (c < mincost)
1258  {
1259  mincost = c;
1260  minperm = perm;
1261  }
1262  progress->enditer(this, mincost, c);
1263  if (period && !(k % period))
1264  {
1265  window++;
1266  }
1267  }
1268  perm = minperm;
1269  place();
1270  }
1271  progress->endorder(this, mincost);
1272 
1273  if (!progress)
1274  {
1275  delete this->progress;
1276  this->progress = 0;
1277  }
1278 }
std::vector< Index >::const_iterator ConstPtr
Definition: gecko.hpp:587
Arc::Index node_end(Node::Index i) const
Definition: gecko.hpp:636
std::vector< Float > weight
Definition: gecko.hpp:735
unsigned int uint
Definition: gecko.hpp:204
std::vector< Node::Index > adj
Definition: gecko.hpp:734
STL namespace.
Arc::Index directed() const
Definition: gecko.cpp:782
std::vector< Node > node
Definition: gecko.hpp:733
uint Index
Definition: gecko.hpp:595
std::vector< Node::Index > perm
Definition: gecko.hpp:732
double b
Definition: lissajous.cpp:42
double Float
Definition: gecko.hpp:208
Arc::Index node_begin(Node::Index i) const
Definition: gecko.hpp:635
Functional * functional
Definition: gecko.hpp:730
Definition: gecko.cpp:343
uint Index
Definition: gecko.hpp:586
void vcycle(uint n, uint work=0)
Definition: gecko.cpp:1169
double a
Definition: lissajous.cpp:41
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
Float cost() const
Definition: gecko.cpp:857
RefCoord t[3]
RefCoord s[3]
Node::Index insert_node(Float length=1)
Definition: gecko.cpp:657
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Progress * progress
Definition: gecko.hpp:731
double f(const Vector &p)