Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__CPU_ALGEBRA__PERMUTATION_UTIL__
34 : #define __H__UG__CPU_ALGEBRA__PERMUTATION_UTIL__
35 :
36 : #include "common/common.h"
37 : #include "common/profiler/profiler.h"
38 : #include "common/error.h"
39 : #include "lib_algebra/ordering_strategies/algorithms/native_cuthill_mckee.h"
40 : #include <vector>
41 :
42 : namespace ug{
43 : /**
44 : * Function to return a permutation of a matrix
45 : * @param[out] PA the permuted matrix PA(perm[r], perm[c]) = A(r, c)
46 : * @param[in] A the input matrix
47 : * @param[in] perm array mapping i -> perm[i]
48 : */
49 : template<typename TMatrix>
50 0 : static void SetMatrixAsPermutation(TMatrix &PA, const TMatrix &A, const std::vector<size_t> &perm)
51 : {
52 : PROFILE_FUNC_GROUP("algebra");
53 0 : PA.resize_and_clear(A.num_rows(), A.num_cols());
54 :
55 0 : for(size_t r=0; r<A.num_rows(); r++)
56 : {
57 0 : size_t Pr=perm[r];
58 0 : for(typename TMatrix::const_row_iterator it = A.begin_row(r); it != A.end_row(r); ++it)
59 : {
60 0 : size_t Pc = perm[it.index()];
61 0 : PA(Pr, Pc) = it.value();
62 : }
63 : }
64 0 : }
65 :
66 :
67 : /**
68 : * Function to compute a permutation of a vector
69 : * @param[out] Pv the permuted vector: Pv[perm[i]] = v[i]
70 : * @param[in] v the input vector
71 : * @param[in] perm array mapping oldindices i -> new index perm[i]
72 : */
73 : template<typename TVector>
74 0 : static void SetVectorAsPermutation(TVector &Pv, const TVector &v, const std::vector<size_t> &perm)
75 : {
76 0 : if(Pv.size() != v.size()) Pv.resize(v.size());
77 0 : for(size_t i=0; i<v.size(); i++)
78 0 : Pv[ perm[i] ] = v[i];
79 0 : }
80 :
81 :
82 : /**
83 : * Function to calculate the inverse of a permutation
84 : * @param[in] perm array mapping i -> perm[i]
85 : * @param[out] invPerm array mapping i -> invPerm[i] so that perm[ invPerm[i] ] = i
86 : * @return true if perm[i] == i forall i.
87 : */
88 : bool GetInversePermutation(const std::vector<size_t> &perm, std::vector<size_t> &invPerm);
89 :
90 : /**
91 : * @param mat A sparse matrix
92 : * @param newIndex the cuthill-mckee ordered new indices
93 : */
94 : // /todo move this to a proper place and remove ordering_strategies dependency
95 :
96 : template<typename TSparseMatrix>
97 : void GetCuthillMcKeeOrder(const TSparseMatrix &mat, std::vector<size_t> &newIndex, bool reverse=true, bool bPreserveConsec=false)
98 : {
99 : std::vector<std::vector<size_t> > neighbors;
100 : neighbors.resize(mat.num_rows());
101 :
102 : for(size_t i=0; i<mat.num_rows(); i++)
103 : {
104 : for(typename TSparseMatrix::const_row_iterator i_it = mat.begin_row(i); i_it != mat.end_row(i); ++i_it)
105 : neighbors[i].push_back(i_it.index());
106 :
107 : // make sure there are no disconnected DoFs
108 : //UG_ASSERT(neighbors[i].size(), "Index "<< i << " does not have any connections. This will most probably "
109 : // "lead to problems and is therefore disallowed.");
110 : }
111 :
112 : ComputeCuthillMcKeeOrder(newIndex, neighbors, reverse, bPreserveConsec);
113 : }
114 : /// @}
115 :
116 :
117 : #ifndef HAVE_BOOL
118 : #define HAVE_BOOL
119 :
120 : class BOOL{
121 : public:
122 : BOOL() : value_(bool()){}
123 : BOOL(bool const& t): value_(t) {}
124 : operator bool() const { return value_; }
125 : private:
126 : char value_;
127 : };
128 :
129 : #endif
130 :
131 : #ifndef HAVE_IS_PERMUTATION
132 : #define HAVE_IS_PERMUTATION
133 :
134 : template <typename O_t>
135 : bool is_permutation(O_t &o){
136 : std::vector<BOOL> container(o.size(), false);
137 : for(unsigned i = 0; i < o.size(); ++i){
138 : if(!container[o[i]]){
139 : container[o[i]] = true;
140 : }
141 : else{
142 : return false; //no doubles allowed
143 : }
144 : }
145 :
146 : for(unsigned i = 0; i < o.size(); ++i){
147 : if(!container[i]){
148 : return false;
149 : }
150 : }
151 :
152 : return true;
153 : }
154 :
155 : #endif
156 :
157 : } // end namespace ug
158 :
159 :
160 : #endif
|
