Line data Source code
1 : /* Copyright (c) 2022: G-CSC, Goethe University Frankfurt
2 : * Felix Salfelder 2022
3 : * This file is part of UG4.
4 : *
5 : * UG4 is free software: you can redistribute it and/or modify it under the
6 : * terms of the GNU Lesser General Public License version 3 (as published by the
7 : * Free Software Foundation) with the following additional attribution
8 : * requirements (according to LGPL/GPL v3 ยง7):
9 : */
10 :
11 : // BGL interface for cpu sparse martrix (dynamic CRS).
12 :
13 : #ifndef UG_SPARSEMATRIX_BOOST_H
14 : #define UG_SPARSEMATRIX_BOOST_H
15 :
16 : #include "common/util/trace.h"
17 : #include "lib_algebra/cpu_algebra/sparsematrix.h"
18 :
19 : #include <boost/graph/properties.hpp> // put_get_helper
20 : #include <boost/iterator/counting_iterator.hpp>
21 :
22 : #include <cstdint> // for std::intmax_t
23 :
24 : namespace boost{
25 :
26 : using ug::iszero; // for now
27 :
28 : // a bit of a pointer dance,
29 : // because ug::SparseMatrix<T>::const_row_iterator is not default constructible.
30 : template<class T>
31 : class SM_adjacency_iterator : public iterator_facade<
32 : SM_adjacency_iterator<T>,
33 : typename ug::SparseMatrix<T>::const_row_iterator,
34 : std::input_iterator_tag,
35 : size_t, // <= reference
36 : std::intmax_t // difference_type
37 : >{ //
38 : typedef ug::SparseMatrix<T> M;
39 : typedef typename M::const_row_iterator iter_t;
40 : typedef iterator_facade<
41 : SM_adjacency_iterator<T>,
42 : typename ug::SparseMatrix<T>::const_row_iterator,
43 : std::input_iterator_tag,
44 : size_t, // <= reference
45 : std::intmax_t // difference_type
46 : > base_class;
47 :
48 : public:
49 : typedef iter_t* value_type;
50 : typedef intmax_t difference_type;
51 : typedef size_t reference;
52 :
53 : public:
54 0 : SM_adjacency_iterator() : base_class(), _base(nullptr), _end(nullptr){
55 : }
56 : SM_adjacency_iterator(SM_adjacency_iterator&& p) = delete;
57 0 : SM_adjacency_iterator(SM_adjacency_iterator const& p)
58 0 : : _base(p._base?(new iter_t(*p._base)) : nullptr)
59 0 : , _end(p._end?(new iter_t(*p._end)) : nullptr){
60 0 : }
61 0 : SM_adjacency_iterator(value_type p, value_type e) : base_class(){
62 : assert(p);
63 0 : _base = new iter_t(*p);
64 0 : _end = new iter_t(*e);
65 0 : skip_zeroes();
66 0 : }
67 0 : ~SM_adjacency_iterator() {
68 0 : delete _base;
69 0 : delete _end;
70 : _base = nullptr;
71 : _end = nullptr;
72 0 : }
73 : SM_adjacency_iterator& operator=(SM_adjacency_iterator&& other) = delete;
74 0 : SM_adjacency_iterator& operator=(const SM_adjacency_iterator& other) {
75 0 : if(other._base){
76 0 : delete _base;
77 0 : _base = new iter_t(*other._base);
78 : }else{ untested();
79 0 : _base = nullptr;
80 : }
81 0 : if(other._end){
82 0 : delete _end;
83 0 : _end = new iter_t(*other._end);
84 : }else{ untested();
85 0 : _end = nullptr;
86 : }
87 0 : return *this;
88 : }
89 :
90 : T const& value() const {
91 : assert(_base);
92 : return _base->value();
93 : }
94 : int row() const {
95 : assert(_base);
96 : return _base->row();
97 : }
98 : size_t idx() const { untested();
99 : assert(_base);
100 : return _base->idx();
101 : }
102 : reference dereference() const {
103 : assert(_base);
104 : return _base->index(); // row?
105 : }
106 : bool operator==(const SM_adjacency_iterator& other) const {
107 : assert(_base);
108 : assert(other._base);
109 : return *_base == *other._base;
110 : }
111 : bool operator!=(const SM_adjacency_iterator& other) const {
112 : assert(_base);
113 : assert(other._base);
114 0 : return *_base != *other._base;
115 : }
116 :
117 : private:
118 : // could use boost::filter_iterator, but does not work yet.
119 0 : void skip_zeroes(){
120 0 : while((*_base) != (*_end)){
121 0 : if(iszero(_base->value())){
122 : ++(*_base);
123 : }else{
124 : break;
125 : }
126 : }
127 0 : }
128 : bool equal(SM_adjacency_iterator const& other) const { untested();
129 : assert(_base);
130 : assert(other._base);
131 : return *_base == *other._base;
132 : }
133 0 : void increment() {
134 : assert(_base);
135 : assert(_end);
136 0 : ++(*_base);
137 0 : skip_zeroes();
138 0 : }
139 : void decrement() { untested();
140 : // incomplete(); // don't use. too complicated.
141 : assert(_base);
142 : --(*_base);
143 : }
144 :
145 : private:
146 : value_type _base;
147 : value_type _end;
148 : friend class iterator_core_access;
149 : }; // SM_adjacency_iterator
150 :
151 : template<class T>
152 : class SM_edge{
153 : public:
154 : // SM_edge(size_t v, size_t w) : _row(v), _idx(w) { untested();
155 : // }
156 : SM_edge(SM_edge const&) = default;
157 : SM_edge(SM_edge&&) = default;
158 : explicit SM_edge() { untested();
159 : }
160 : SM_edge(int s, int t, T const& v) : _source(s), _target(t), _value(v) {
161 : }
162 : SM_edge& operator=(SM_edge const&) = default;
163 : SM_edge& operator=(SM_edge&&) = default;
164 :
165 : public:
166 : int source() const{
167 : return _source;
168 : }
169 : int target() const{
170 : return _target;
171 : }
172 : template<class X>
173 : T const& value(X const&) const{
174 : return _value;
175 : }
176 :
177 : public:
178 : bool operator==(const SM_edge& other) const { untested();
179 : return _source == other._source && _target == other._target;
180 : }
181 : bool operator!=(const SM_edge& other) const { untested();
182 : return !operator==(other);
183 : }
184 :
185 : private:
186 : int _source;
187 : int _target;
188 : T _value;
189 : }; // SM_edge
190 :
191 : template<class T, class M>
192 : int source(SM_edge<T> const& e, M const&)
193 : { untested();
194 : return e.source();
195 : }
196 :
197 : template<class T, class M>
198 : int target(SM_edge<T> const& e, M const&)
199 : { untested();
200 : return e.target();
201 : }
202 :
203 : template<class T, bool out=true>
204 : class SM_out_edge_iterator : public iterator_facade<
205 : SM_out_edge_iterator<T, out>,
206 : SM_adjacency_iterator<T>,
207 : // bidirectional_traversal_tag, // breaks InputIterator (why?)
208 : std::input_iterator_tag,
209 : SM_edge<T>, // <= reference
210 : std::intmax_t // difference_type
211 : >{ //
212 : public: // types
213 : typedef iterator_facade<
214 : SM_out_edge_iterator<T, out>,
215 : SM_adjacency_iterator<T>,
216 : std::input_iterator_tag,
217 : SM_edge<T>, // <= reference
218 : std::intmax_t // difference_type
219 : > base_class;
220 : typedef ug::SparseMatrix<T> M;
221 : typedef SM_adjacency_iterator<T> value_type;
222 : typedef intmax_t difference_type;
223 : typedef SM_edge<T> reference;
224 : typedef SM_edge<T> edge_type;
225 :
226 : public: // construct
227 0 : explicit SM_out_edge_iterator() : base_class() {
228 : }
229 0 : explicit SM_out_edge_iterator(SM_adjacency_iterator<T> w, int src)
230 0 : : base_class(), _base(w), _src(src) {
231 : if(out){
232 : }else{
233 : }
234 : }
235 0 : /* explicit */ SM_out_edge_iterator(SM_out_edge_iterator const& p)
236 0 : : base_class(p), _base(p._base), _src(p._src){
237 : if(out){
238 : }else{
239 : }
240 : }
241 : SM_out_edge_iterator(SM_out_edge_iterator&& p) = delete; // why?
242 : ~SM_out_edge_iterator(){
243 0 : }
244 : SM_out_edge_iterator& operator=(SM_out_edge_iterator&& other) = default;
245 0 : SM_out_edge_iterator& operator=(SM_out_edge_iterator const& other) = default;
246 : #if 0
247 : public: // op
248 : SM_out_edge_iterator operator+(int a) const{ untested();
249 : SM_out_edge_iterator ret(*this);
250 : ret.base.first += a;
251 : return ret;
252 : }
253 : SM_out_edge_iterator operator-(int a) const{ untested();
254 : SM_out_edge_iterator ret(*this);
255 : ret.base.first -= a;
256 : return ret;
257 : }
258 : difference_type operator-(SM_out_edge_iterator const& other) const{ untested();
259 : SM_out_edge_iterator ret(*this);
260 : return ret.base.first - other.base.first;
261 : }
262 : #endif
263 : private:
264 : reference dereference() const {
265 : if(out){
266 : return edge_type(_src, *_base, _base.value());
267 : }else{
268 : return edge_type(*_base, _src, _base.value());
269 : }
270 : }
271 : bool equal(SM_out_edge_iterator const& other) const { untested();
272 : assert(_base.first == other._base.first);
273 : return _base.second == other._base.second;
274 : }
275 : void increment() { untested();
276 : ++_base;
277 : }
278 : void decrement() { untested();
279 : --_base;
280 : }
281 : // bool operator==(const SM_out_edge_iterator& other) const
282 : // { incomplete();
283 : // return false;
284 : // }
285 : public:
286 : SM_out_edge_iterator& operator++() {
287 : ++_base;
288 0 : return *this;
289 : }
290 : SM_out_edge_iterator operator++(int) { untested();
291 : SM_out_edge_iterator copy(*this);
292 : ++*this;
293 : return copy;
294 : }
295 : bool operator==(const SM_out_edge_iterator& other) const {
296 : return _base == other._base;
297 : }
298 : bool operator!=(const SM_out_edge_iterator& other) const {
299 : return _base != other._base;
300 : }
301 :
302 : private:
303 : value_type _base;
304 : int _src;
305 : friend class iterator_core_access;
306 : }; // SM_out_edge_iterator
307 :
308 : struct SM_traversal_tag
309 : : adjacency_graph_tag, bidirectional_graph_tag, vertex_list_graph_tag {};
310 :
311 : template <class T> struct graph_traits<ug::SparseMatrix<T>>{
312 : typedef int vertex_descriptor;
313 : typedef SM_edge<T> edge_descriptor;
314 : typedef directed_tag directed_category;
315 : typedef disallow_parallel_edge_tag edge_parallel_category;
316 : typedef SM_traversal_tag traversal_category;
317 : typedef counting_iterator<size_t> vertex_iterator;
318 : typedef SM_out_edge_iterator<T> out_edge_iterator;
319 : typedef SM_adjacency_iterator<T> adjacency_iterator;
320 : typedef int degree_size_type;
321 : typedef int vertices_size_type;
322 : };
323 :
324 : template<class T>
325 : std::pair<typename graph_traits<ug::SparseMatrix<T>>::adjacency_iterator,
326 : typename graph_traits<ug::SparseMatrix<T>>::adjacency_iterator>
327 0 : inline adjacent_vertices(size_t v, ug::SparseMatrix<T> const& M)
328 : {
329 : assert(v<M.num_rows());
330 : typedef typename graph_traits<ug::SparseMatrix<T>>::adjacency_iterator a;
331 :
332 : typename ug::SparseMatrix<T>::const_row_iterator b(M.begin_row(v));
333 : typename ug::SparseMatrix<T>::const_row_iterator e(M.end_row(v));
334 :
335 0 : return std::make_pair(a(&b, &e), a(&e, &e));
336 : }
337 :
338 : template<class T>
339 : inline std::pair<SM_out_edge_iterator<T>, SM_out_edge_iterator<T>>
340 0 : out_edges(int v, ug::SparseMatrix<T> const& g)
341 : {
342 : assert(size_t(v)<g.num_rows());
343 : typedef SM_out_edge_iterator<T> Iter;
344 0 : auto a = adjacent_vertices(v, g);
345 0 : return std::make_pair(Iter(a.first, v), Iter(a.second, v));
346 : }
347 :
348 : template<class T>
349 : class sparse_matrix_index_map
350 : : public put_get_helper<size_t, sparse_matrix_index_map<T> > { //
351 : public:
352 : typedef size_t vertex_index_type;
353 : typedef size_t vertex_descriptor;
354 : typedef readable_property_map_tag category;
355 : typedef vertex_index_type value_type;
356 : typedef vertex_index_type reference;
357 : typedef vertex_descriptor key_type;
358 :
359 : sparse_matrix_index_map(sparse_matrix_index_map const& p) {
360 : }
361 : sparse_matrix_index_map(ug::SparseMatrix<T>const&, boost::vertex_index_t) { untested();
362 : }
363 : template<class X>
364 : sparse_matrix_index_map(X const&) {
365 : }
366 : template <class T_>
367 : value_type operator[](T_ x) const {
368 : return x;
369 : }
370 : sparse_matrix_index_map& operator=(const sparse_matrix_index_map& s) { untested();
371 : return *this;
372 : }
373 : };
374 :
375 : template<class T>
376 : struct property_map<ug::SparseMatrix<T>, vertex_index_t>{
377 : typedef sparse_matrix_index_map<T> type;
378 : typedef type const_type;
379 : };
380 :
381 : template<class T>
382 : typename property_map<ug::SparseMatrix<T>, vertex_index_t>::const_type
383 : get(vertex_index_t, ug::SparseMatrix<T> const& m)
384 : {
385 : return sparse_matrix_index_map<T>(m);
386 : }
387 :
388 : template<class T>
389 : std::pair<counting_iterator<size_t>, counting_iterator<size_t> > vertices(
390 : ug::SparseMatrix<T> const& M)
391 : {
392 : counting_iterator<size_t> b(0);
393 : counting_iterator<size_t> e(M.num_rows());
394 :
395 : return std::make_pair(b,e);
396 : }
397 :
398 : template<class T>
399 : int num_vertices(ug::SparseMatrix<T> const& M)
400 : {
401 : assert(M.num_rows() == M.num_cols());
402 : return M.num_rows();
403 : }
404 :
405 : template<class T>
406 0 : int out_degree(int v, ug::SparseMatrix<T> const& M)
407 : {
408 : int c = 0;
409 0 : auto i = out_edges(v, M);
410 0 : for(; i.first != i.second; ++i.first) {
411 0 : ++c;
412 : }
413 0 : return c;
414 : }
415 :
416 : // template<class T>
417 : // T get(edge_weight_t, ug::SparseMatrix<T> const& M, SM_edge<T> e)
418 : // { untested();
419 : // return e.value();
420 : // }
421 :
422 : template<class T, class M=ug::SparseMatrix<T>>
423 : class SM_edge_weight_map :
424 : public put_get_helper< T /*algebraic connection?*/, SM_edge_weight_map<T> > { //
425 : public:
426 : typedef T edge_weight_type;
427 : typedef int vertex_descriptor;
428 : typedef readable_property_map_tag category;
429 : typedef edge_weight_type value_type;
430 : typedef T& reference;
431 : typedef vertex_descriptor key_type;
432 :
433 : SM_edge_weight_map(SM_edge_weight_map const& p) : _g(p._g) { untested();
434 : }
435 : SM_edge_weight_map(ug::SparseMatrix<T>const & g, boost::edge_weight_t) : _g(g) { untested();
436 : }
437 : // bug?
438 : SM_edge_weight_map(M const& g) : _g(g) {
439 : }
440 : template <class X>
441 : value_type operator[](X const& x) const {
442 : // assert(x == _g.position(x));
443 : return x.value(_g);
444 : }
445 : SM_edge_weight_map& operator=(const SM_edge_weight_map& s) { untested();
446 : assert(&s._g==&_g); (void)s;
447 : return *this;
448 : }
449 : private:
450 : M const& _g;
451 : }; // SM_edge_weight_map
452 :
453 : // g++ 'enumeral_type' in template unification not implemented workaround
454 : template<class T, class Tag>
455 : struct property_map<ug::SparseMatrix<T>, Tag> {
456 : // typedef typename gala_graph_property_map<Tag>::template bind_<T> map_gen;
457 : // typedef typename map_gen::type type;
458 : // typedef typename map_gen::const_type const_type;
459 : };
460 :
461 : template<class T>
462 : inline SM_edge_weight_map<T>
463 : //inline typename property_map<ug::SparseMatrix<T>, edge_weight_t>::type
464 : get(edge_weight_t, ug::SparseMatrix<T> const & g) {
465 : return SM_edge_weight_map<T>(g);
466 : }
467 :
468 : #if 0
469 : template<class T>
470 : inline typename property_map<ug::SparseMatrix<T>, edge_all_t>::type
471 : get(edge_all_t, ug::SparseMatrix<T> & g) { incomplete();
472 : typedef typename property_map<ug::SparseMatrix<T>, edge_all_t>::type
473 : pmap_type;
474 : return pmap_type(&g);
475 : }
476 : #endif
477 :
478 : } // boost
479 :
480 : namespace ug{
481 :
482 : // used from boost::print_graph, graph_utility.hpp.
483 : // must be in ug, because of ADL. why don't they call boost::{out_edges,vertices}?
484 : using boost::vertices;
485 : using boost::out_edges;
486 :
487 : }// ug
488 :
489 : #endif // guard
|
