Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
parallel_matrix_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 * Author: Andreas Vogel
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__LIB_ALGEBRA__PARALLELIZATION__PARALLEL_MATRIX_IMPL__
34#define __H__LIB_ALGEBRA__PARALLELIZATION__PARALLEL_MATRIX_IMPL__
35
36#include "parallel_matrix.h"
37
38namespace ug
39{
40
41template <typename TMatrix>
44{
45// forward to sequential matrices
46 TMatrix::operator= (*dynamic_cast<const TMatrix*>(&M));
47
48// copy storage type and layouts
49 this->set_storage_type(M.get_storage_mask());
50 this->set_layouts(M.layouts());
51
52// we're done
53 return *this;
54}
55
56template <typename TMatrix>
57bool
60{
61 // can only change if current state is defined
62 if(has_storage_type(PST_UNDEFINED))
63 {
64 UG_LOG("Current Storage Type is undefined. "
65 "Cannot change storage type.\n");
66 return false;
67 }
68
69 // if already in that type
70 if(has_storage_type(type))
71 return true;
72
73// todo: Implement more types.
74 UG_THROW("ParallelMatrix::change_storage_type:"
75 " Currently no storage conversion supported.");
76}
77
78// calculate res = A x
79template <typename TMatrix>
80template<typename TPVector>
81bool
83apply(TPVector &res, const TPVector &x) const
84{
85 PROFILE_FUNC_GROUP("algebra");
86// check types combinations
87 int type = -1;
88 if(has_storage_type(PST_ADDITIVE)
89 && x.has_storage_type(PST_CONSISTENT)) type = 0;
90 if(has_storage_type(PST_CONSISTENT)
91 && x.has_storage_type(PST_ADDITIVE)) type = 1;
92 if(has_storage_type(PST_CONSISTENT)
93 && x.has_storage_type(PST_CONSISTENT)) type = 2;
94
95// if no admissible type is found, return error
96 if(type == -1)
97 {
98 UG_THROW("ParallelMatrix::apply (b = A*x): "
99 "Wrong storage type of Matrix/Vector: Possibilities are:\n"
100 " - A is PST_ADDITIVE and x is PST_CONSISTENT\n"
101 " - A is PST_CONSISTENT and x is PST_ADDITIVE\n"
102 " (storage type of A = " << get_storage_type() << ", x = " << x.get_storage_type() << ")");
103 }
104
105// apply on single process vector
106 TMatrix::axpy(res, 0.0, res, 1.0, x);
107
108// set outgoing vector to additive storage
109 switch(type)
110 {
111 case 0: res.set_storage_type(PST_ADDITIVE); break;
112 case 1: res.set_storage_type(PST_ADDITIVE); break;
113 case 2: res.set_storage_type(PST_CONSISTENT); break;
114 }
115
116// we're done.
117 return true;
118}
119
120// calculate res = A.T x
121template <typename TMatrix>
122template<typename TPVector>
123bool
125apply_transposed(TPVector &res, const TPVector &x) const
126{
127 PROFILE_FUNC_GROUP("algebra");
128// check types combinations
129 int type = -1;
130 if(has_storage_type(PST_ADDITIVE)
131 && x.has_storage_type(PST_CONSISTENT)) type = 0;
132 if(has_storage_type(PST_CONSISTENT)
133 && x.has_storage_type(PST_ADDITIVE)) type = 1;
134 if(has_storage_type(PST_CONSISTENT)
135 && x.has_storage_type(PST_CONSISTENT)) type = 2;
136
137// if no admissible type is found, return error
138 if(type == -1)
139 {
140 UG_THROW("ParallelMatrix::apply_transposed (b = A^T*x): "
141 "Wrong storage type of Matrix/Vector: Possibilities are:\n"
142 " - A is PST_ADDITIVE and x is PST_CONSISTENT\n"
143 " - A is PST_CONSISTENT and x is PST_ADDITIVE\n"
144 " (storage type of A = " << get_storage_type() << ", x = " << x.get_storage_type() << ")");
145 }
146
147// apply on single process vector
148 TMatrix::axpy_transposed(res, 0.0, res, 1.0, x);
149
150// set outgoing vector to additive storage
151 switch(type)
152 {
153 case 0: res.set_storage_type(PST_ADDITIVE); break;
154 case 1: res.set_storage_type(PST_ADDITIVE); break;
155 case 2: res.set_storage_type(PST_CONSISTENT); break;
156 }
157
158// we're done.
159 return true;
160}
161
162// calculate res -= A x
163template <typename TMatrix>
164template<typename TPVector>
165bool
167matmul_minus(TPVector &res, const TPVector &x) const
168{
169 PROFILE_FUNC_GROUP("algebra");
170// check types combinations
171 int type = -1;
172 if(this->has_storage_type(PST_ADDITIVE)
173 && x.has_storage_type(PST_CONSISTENT)
174 && res.has_storage_type(PST_ADDITIVE)) type = 0;
175
176// if no admissible type is found, return error
177 if(type == -1)
178 {
179 UG_THROW("ParallelMatrix::matmul_minus (b -= A*x):"
180 " Wrong storage type of Matrix/Vector: Possibilities are:\n"
181 " - A is PST_ADDITIVE and x is PST_CONSISTENT and b is PST_ADDITIVE\n"
182 " (storage type of A = " << this->get_storage_type() << ", x = " << x.get_storage_type() << ", b = " << res.get_storage_type() << ")");
183 }
184
185// apply on single process vector
186 TMatrix::axpy(res, 1.0, res, -1.0, x);
187
188// set outgoing vector to additive storage
189// (it could have been PST_UNIQUE before)
190 switch(type)
191 {
192 case 0: res.set_storage_type(PST_ADDITIVE); break;
193 }
194
195// we're done.
196 return true;
197}
198
199
200template<typename matrix_type, typename vector_type>
202{
207 && x.has_storage_type(PST_ADDITIVE)) type = PST_ADDITIVE; // 1
210
211 // if no admissible type is found, return error
212 if(type == PST_UNDEFINED)
213 {
214 UG_THROW("ParallelMatrix::apply (b = A*x): "
215 "Wrong storage type of Matrix/Vector: Possibilities are:\n"
216 " - A is PST_ADDITIVE and x is PST_CONSISTENT\n"
217 " - A is PST_CONSISTENT and x is PST_ADDITIVE\n"
218 " (storage type of A = " << A1.get_storage_type() << ", x = " << x.get_storage_type() << ")");
219 }
220 return type;
221}
222
223template<typename matrix_type, typename vector_type>
225 const number &beta1, const ParallelMatrix<matrix_type> &A1, const ParallelVector<vector_type> &w1)
226{
227 PROFILE_FUNC_GROUP("algebra");
228 // check types combinations
230 if(type == PST_UNDEFINED) return false;
231
232 MatMult(dynamic_cast<vector_type&>(dest), beta1, dynamic_cast<const matrix_type&>(A1), dynamic_cast<const vector_type&>(w1));
233
234 // set outgoing vector storage type
235 dest.set_storage_type(type);
236
237 // we're done.
238 return true;
239}
240
241template<typename matrix_type, typename vector_type>
243 const number &alpha1, const ParallelVector<vector_type> &v1,
244 const number &beta1, const ParallelMatrix<matrix_type> &A1, const ParallelVector<vector_type> &w1)
245{
246 PROFILE_FUNC_GROUP("algebra");
247 // check types combinations
249 if(type == PST_UNDEFINED) return false;
250
251 if(!v1.has_storage_type(type))
252 {
253 UG_THROW("MatMultAdd(dest, alpha1, v1, beta1, A1, w1): Storage type of A1*w1 = " << type << ", storage type of v1 = " << v1.get_storage_type() << ".");
254 }
255
256 MatMultAdd(dynamic_cast<vector_type&>(dest), alpha1, dynamic_cast<const vector_type&>(v1),
257 beta1, dynamic_cast<const matrix_type&>(A1), dynamic_cast<const vector_type&>(w1));
258
259 // set outgoing vector storage type
260 dest.set_storage_type(type);
261
262 // we're done.
263 return true;
264}
265
266template<typename matrix_type, typename vector_type>
268 const number &alpha1, const ParallelVector<vector_type> &v1,
269 const number &alpha2, const ParallelVector<vector_type> &v2,
270 const number &beta1, const ParallelMatrix<matrix_type> &A1, const ParallelVector<vector_type> &w1)
271{
272 PROFILE_FUNC_GROUP("algebra");
273 // check types combinations
275 if(type == PST_UNDEFINED) return false;
276
277 if(!v1.has_storage_type(type) || !v2.has_storage_type(type))
278 {
279 UG_THROW("MatMultAdd(dest, alpha1, v1, alpha2, v2, beta1, A1, w1):"
280 " Storage type of A1*w1 doesnt match storage type of v1 or v2.\n"
281 "Storage type of A1*w1 = " << type << " v1 = " << v1.get_storage_type() << ", v2 = " << v2.get_storage_type() << ".");
282 }
283
284 MatMultAdd(dynamic_cast<vector_type&>(dest), alpha1, dynamic_cast<const vector_type&>(v1), alpha2, dynamic_cast<const vector_type&>(v2),
285 beta1, dynamic_cast<const matrix_type&>(A1), dynamic_cast<const vector_type&>(w1));
286
287 // set outgoing vector storage type
288 dest.set_storage_type(type);
289
290 // we're done.
291 return true;
292}
293
294template<typename matrix_type, typename vector_type>
296 const number &beta1, const ParallelMatrix<matrix_type> &A1, const ParallelVector<vector_type> &w1)
297{
298 PROFILE_FUNC_GROUP("algebra");
299 // check types combinations
301 if(type == PST_UNDEFINED) UG_THROW("Wrong storage type");
302
303 MatMultTransposed(dynamic_cast<vector_type&>(dest), beta1, dynamic_cast<const matrix_type&>(A1), dynamic_cast<const vector_type&>(w1));
304
305 // set outgoing vector storage type
306 dest.set_storage_type(type);
307}
308
309
310} // end namespace ug
311
312#endif /* __H__LIB_ALGEBRA__PARALLELIZATION__PARALLEL_MATRIX_IMPL__ */
Wrapper for sequential matrices to handle them in parallel.
Definition parallel_matrix.h:65
bool has_storage_type(uint type) const
returns if the current storage type has a given representation
Definition parallel_matrix.h:116
bool apply_transposed(TPVector &res, const TPVector &x) const
calculate res = A.T x
Definition parallel_matrix_impl.h:125
bool matmul_minus(TPVector &res, const TPVector &x) const
calculate res -= A x
Definition parallel_matrix_impl.h:167
ParallelStorageType get_storage_type() const
Definition parallel_matrix.h:121
ConstSmartPtr< AlgebraLayouts > layouts() const
returns the algebra layouts
Definition parallel_matrix.h:94
bool apply(TPVector &res, const TPVector &x) const
calculate res = A x
Definition parallel_matrix_impl.h:83
uint get_storage_mask() const
returns storage type mask
Definition parallel_matrix.h:120
bool change_storage_type(ParallelStorageType type)
changes to the requested storage type if possible
Definition parallel_matrix_impl.h:59
ParallelMatrix< TMatrix > this_type
own type
Definition parallel_matrix.h:76
this_type & operator=(const this_type &M)
assignment
Definition parallel_matrix_impl.h:43
Definition parallel_vector.h:60
bool has_storage_type(uint type) const
returns if the current storage type has a given representation
Definition parallel_vector.h:119
ParallelStorageType get_storage_type() const
Definition parallel_vector.h:124
void set_storage_type(uint type)
sets the storage type
Definition parallel_vector.h:104
ParallelStorageType
Definition parallel_storage_type.h:66
@ PST_UNDEFINED
Definition parallel_storage_type.h:67
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
void MatMultTransposedDirect(ParallelVector< vector_type > &dest, const number &beta1, const ParallelMatrix< matrix_type > &A1, const ParallelVector< vector_type > &w1)
Definition parallel_matrix_impl.h:295
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:68
bool MatMultDirect(ParallelVector< vector_type > &dest, const number &beta1, const ParallelMatrix< matrix_type > &A1, const ParallelVector< vector_type > &w1)
Definition parallel_matrix_impl.h:224
bool MatMultAddDirect(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1, const number &alpha1, const vector_t &v1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition matrix_use_global_functions.h:43
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
ug::ParallelStorageType GetMultType(const ParallelMatrix< matrix_type > &A1, const ParallelVector< vector_type > &x)
Definition parallel_matrix_impl.h:201
bool MatMultTransposed(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:111
#define PROFILE_FUNC_GROUP(groups)
Definition profiler.h:258