ug4
Loading...
Searching...
No Matches
block_dense.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2013: 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
34#ifndef __H__UG__SMALL_ALGEBRA__BLOCK_DENSE__
35#define __H__UG__SMALL_ALGEBRA__BLOCK_DENSE__
36
37#include "densematrix.h"
38#include "densevector.h"
39#include <algorithm>
40
41namespace ug{
42
43template<typename A>
44inline double BlockNorm2(const DenseMatrix<A> &mat)
45{
46 double sum=0;
47 for(size_t r=0; r < mat.num_rows(); r++)
48 for(size_t c=0; c < mat.num_cols(); c++)
49 sum += BlockNorm2(mat(r, c));
50
51 return sum;
52}
53
54
55template<typename A>
56inline double BlockNorm2(const DenseVector<A> &v)
57{
58 double sum=0;
59 for(size_t i=0; i < v.size(); i++)
60 sum += BlockNorm2(v[i]);
61 return sum;
62}
63
64
65
66template<typename A>
67inline double BlockMaxNorm(const DenseVector<A> &v)
68{
69 double max=0;
70 for(size_t i=0; i < v.size(); i++)
71 max = std::max(max, BlockMaxNorm(v[i]));
72 return max;
73}
74
75
76
77
79// algebra stuff to avoid temporary variables
80
81// vec = mat * vec
82// todo: replace add_mult etc. with template expressions
83// dest = b*vec
84template<typename A, typename B, typename C>
85inline void AssignMult(DenseVector<A> &dest, const DenseMatrix<B> &mat, const DenseVector<C> &vec)
86{
87 UG_ASSERT(mat.num_cols() == vec.size(), "");
88 dest.resize(mat.num_rows());
89 for(size_t r=0; r < mat.num_rows(); r++)
90 {
91 AssignMult(dest(r), mat(r, 0), vec(0));
92 for(size_t c=1; c < mat.num_cols(); c++)
93 AddMult(dest(r), mat(r, c), vec(c));
94 }
95}
96
97template<typename A, typename B, typename C>
98inline void AssignMult(DenseMatrix<A> &dest, const DenseMatrix<B> &mA, const DenseMatrix<C> &mB)
99{
100 UG_ASSERT(mA.num_cols() == mB.num_rows(), "");
101 dest.resize(mA.num_rows(), mB.num_cols());
102 for(size_t r=0; r < mA.num_rows(); r++)
103 for(size_t c=0; c < mB.num_cols(); c++)
104 {
105 AssignMult(dest(r, c), mA(r, 0), mB(0, c));
106 for(size_t k=1; k < mB.num_rows(); k++)
107 AddMult(dest(r, c), mA(r, k), mB(k, c));
108 }
109}
110
111// dest += b*vec
112template<typename A, typename B, typename C>
113inline void AddMult(DenseVector<A> &dest, const DenseMatrix<B> &mat, const DenseVector<C> &vec)
114{
115 UG_ASSERT(mat.num_cols() == vec.size(), "");
116 dest.resize(mat.num_rows());
117 for(size_t r=0; r < mat.num_rows(); r++)
118 {
119 for(size_t c=0; c < mat.num_cols(); c++)
120 AddMult(dest(r), mat(r, c), vec(c));
121 }
122}
123
124
125// dest += b*vec
126template<typename A, typename B, typename C>
127inline void AddMult(DenseMatrix<A> &dest, const DenseMatrix<B> &mA, const DenseMatrix<C> &mB)
128{
129 UG_ASSERT(mA.num_cols() == mB.num_rows(), "");
130 UG_ASSERT(dest.num_rows()==mA.num_rows() && dest.num_cols()==mB.num_cols(), "");
131 for(size_t r=0; r < mA.num_rows(); r++)
132 for(size_t c=0; c < mB.num_cols(); c++)
133 {
134 for(size_t k=0; k < mB.num_rows(); k++)
135 AddMult(dest(r, c), mA(r, k), mB(k, c));
136 }
137}
138
139
140
141// dest -= b*vec
142template<typename A, typename B, typename C>
143inline void SubMult(DenseVector<A> &dest, const DenseMatrix<B> &mat, const DenseVector<C> &vec)
144{
145 UG_ASSERT(mat.num_cols() == vec.size(), "");
146 dest.resize(mat.num_rows());
147 for(size_t r=0; r < mat.num_rows(); r++)
148 {
149 for(size_t c=0; c < mat.num_cols(); c++)
150 SubMult(dest(r), mat(r, c), vec(c));
151 }
152}
153
154
155// mat = alpha * mat
156
157// dest = b*vec
158template<typename A, typename B>
159inline void AssignMult(DenseMatrix<A> &dest, const double &alpha, const DenseMatrix<B> &mat)
160{
161 dest.resize(mat.num_rows(), mat.num_cols());
162 for(size_t r=0; r < mat.num_rows(); r++)
163 {
164 for(size_t c=0; c < mat.num_cols(); c++)
165 AssignMult(dest(r, c), alpha, mat(r, c));
166 }
167}
168
169template<typename A, typename B>
170inline void AddMult(DenseMatrix<A> &dest, const double &alpha, const DenseMatrix<B> &mat)
171{
172 dest.resize(mat.num_rows(), mat.num_cols());
173 for(size_t r=0; r < mat.num_rows(); r++)
174 {
175 for(size_t c=0; c < mat.num_cols(); c++)
176 AddMult(dest(r, c), alpha, mat(r, c));
177 }
178}
179
180template<typename A, typename B>
181inline void SubMult(DenseMatrix<A> &dest, const double &alpha, const DenseMatrix<B> &mat)
182{
183 dest.resize(mat.num_rows(), mat.num_cols());
184 for(size_t r=0; r < mat.num_rows(); r++)
185 {
186 for(size_t c=0; c < mat.num_cols(); c++)
187 SubMult(dest(r, c), alpha, mat(r, c));
188 }
189}
190
191
192// VECTORs
193
194// dest = b*vec
195template<typename A, typename B>
196inline void AssignMult(DenseVector<A> &dest, const double &b, const DenseVector<B> &vec)
197{
198 dest.resize(vec.size());
199 for(size_t i=0; i<vec.size(); i++)
200 AssignMult(dest[i], b, vec[i]);
201}
202
203// dest += b*vec
204template<typename A, typename B>
205inline void AddMult(DenseVector<A> &dest, const double &b, const A &vec)
206{
207 dest.resize(vec.size());
208 for(size_t i=0; i<vec.size(); i++)
209 AddMult(dest[i], b, vec[i]);
210}
211
212// dest -= b*vec
213template<typename A, typename B>
214inline void SubMult(DenseVector<A> &dest, const double &b, const DenseVector<B> &vec)
215{
216 dest.resize(vec.size());
217 for(size_t i=0; i<vec.size(); i++)
218 SubMult(dest[i], b, vec[i]);
219}
220
221// dest = vec*b
222template<typename A> inline void AssignMult(A &dest, const A &vec, const double &b)
223{
224 AssignMult(dest, b, vec);
225}
226// dest += vec*b
227template<typename A> inline void AddMult(A &dest, const A &vec, const double &b)
228{
229 AddMult(dest, b, vec);
230}
231// dest -= vec*b
232template<typename A> inline void SubMult(A &dest, const A &vec, const double &b)
233{
234 SubMult(dest, b, vec);
235}
236
237
238template<typename T>
239inline void SetSize(DenseMatrix<T> &t, size_t a, size_t b)
240{
241 t.resize(a, b);
242}
243
244//setSize(t, a) for vectors
245template<typename T>
246inline void SetSize(DenseVector<T> &t, size_t a)
247{
248 t.resize(a);
249}
250
251// getSize
252template<typename T>
253inline size_t GetSize(const DenseVector<T> &t)
254{
255 return t.size();
256}
257
258//getRows
259template<typename T>
260inline size_t GetRows(const DenseMatrix<T> &t)
261{
262 return t.num_rows();
263}
264
265template<typename T>
266inline size_t GetCols(const DenseMatrix<T> &t)
267{
268 return t.num_cols();
269}
270
271template<typename T>
272struct block_traits;
273
274template<typename T>
276{
277 enum { ordering = DenseMatrix<T>::ordering };
279 enum { static_num_rows = DenseMatrix<T>::static_num_rows};
280 enum { static_num_cols = DenseMatrix<T>::static_num_cols};
281
282 // todo: to be implemented
283 //typedef DenseMatrixInverse inverse_type;
284};
285
286template<typename T>
292
293template<typename T1, typename T2>
294struct block_multiply_traits;
295
296template<typename T1, typename T2>
301
303// block_traits
304
305template<typename T> class DenseMatrixInverse;
306
308// variable matrix
309template<eMatrixOrdering TOrdering>
311{
312 enum { ordering = TOrdering };
313 enum { is_static = false};
314 enum { static_num_rows = 0};
315 enum { static_num_cols = 0};
316 enum { depth = 1 };
317
319};
321// fixed matrix
322template<size_t TBlockSize, eMatrixOrdering TOrdering>
323struct block_traits< DenseMatrix< FixedArray2<number, TBlockSize, TBlockSize, TOrdering> > >
324{
325 enum { ordering = TOrdering };
326 enum { is_static = true};
327 enum { static_num_rows = TBlockSize};
328 enum { static_num_cols = TBlockSize};
329 enum { depth = 1 };
330
332};
333
335// variable block matrix
336template<typename TValue, eMatrixOrdering TOrdering>
337struct block_traits< DenseMatrix< VariableArray2<TValue, TOrdering> > >
338{
339 enum { ordering = TOrdering };
340 enum { is_static = false};
341 enum { static_num_rows = 0};
342 enum { static_num_cols = 0};
343 enum { depth = block_traits<TValue>::depth+1 };
344
346};
347
348template<typename TValue, size_t TBlockSize, eMatrixOrdering TOrdering>
349struct block_traits< DenseMatrix< FixedArray2<TValue, TBlockSize, TBlockSize, TOrdering> > >
350{
351 enum { ordering = TOrdering };
352 enum { is_static = false};
353 enum { static_num_rows = 0};
354 enum { static_num_cols = 0};
355 enum { depth = block_traits<TValue>::depth+1 };
356
358};
359
361// fixed 1x1 to 3x3 : inverse is matrix
362template<eMatrixOrdering TOrdering>
363struct block_traits< DenseMatrix< FixedArray2<number, 1, 1, TOrdering> > >
364{
365 enum { ordering = DenseMatrix< FixedArray2<number, 1, 1> >::ordering };
366 enum { is_static = true};
367 enum { static_num_rows = 1};
368 enum { static_num_cols = 1};
369
371};
372
373template<eMatrixOrdering TOrdering>
374struct block_traits< DenseMatrix< FixedArray2<number, 2, 2, TOrdering> > >
375{
376 enum { ordering = DenseMatrix< FixedArray2<number, 2, 2> >::ordering };
377 enum { is_static = true};
378 enum { static_num_rows = 2};
379 enum { static_num_cols = 2};
380
382};
383
384template<eMatrixOrdering TOrdering>
385struct block_traits< DenseMatrix< FixedArray2<number, 3, 3, TOrdering> > >
386{
387 enum { ordering = DenseMatrix< FixedArray2<number, 3, 3> >::ordering };
388 enum { is_static = true};
389 enum { static_num_rows = 3};
390 enum { static_num_cols = 3};
391
393};
394
395
396template<typename T> struct block_multiply_traits<DenseMatrix<T>, DenseMatrix<T> >
397{
399};
400
401
403
404
405
406
407} // namespace ug
408
409#endif // __H__UG__SMALL_ALGEBRA__BLOCK_DENSE__
Definition densematrix.h:57
Definition lapack_densematrix_inverse.h:51
Definition densevector.h:101
Definition fixed_array.h:135
Definition variable_array.h:139
#define UG_ASSERT(expr, msg)
Definition assert.h:70
double number
Definition types.h:124
the ug namespace
size_t GetCols(const T &t)
void AssignMult(A &dest, const B &b, const C &vec)
void SubMult(A &dest, const B &b, const C &vec)
void SetSize(T &t, size_t a, size_t b)
size_t GetRows(const T &t)
number BlockMaxNorm(const T &t)
double BlockNorm2(const TYPE &v)
Definition blocks.h:51
void AddMult(A &dest, const B &b, const C &vec)
size_t GetSize(const T &t)
DenseVector< T2 > ReturnType
Definition block_dense.h:299
DenseMatrix< T > ReturnType
Definition block_dense.h:398
Definition smallalgebra_interface.h:42
DenseMatrixInverse< VariableArray2< number, TOrdering > > inverse_type
Definition block_dense.h:357
DenseMatrix< FixedArray2< number, 1, 1, TOrdering > > inverse_type
Definition block_dense.h:370
DenseMatrix< FixedArray2< number, 2, 2, TOrdering > > inverse_type
Definition block_dense.h:381
DenseMatrix< FixedArray2< number, 3, 3, TOrdering > > inverse_type
Definition block_dense.h:392
DenseMatrixInverse< FixedArray2< number, TBlockSize, TBlockSize, TOrdering > > inverse_type
Definition block_dense.h:331
DenseMatrixInverse< VariableArray2< number, TOrdering > > inverse_type
Definition block_dense.h:345
DenseMatrixInverse< VariableArray2< number, TOrdering > > inverse_type
Definition block_dense.h:318
Definition communication_policies.h:58
@ is_static
Definition communication_policies.h:60