Line data Source code
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 :
41 : namespace ug{
42 :
43 : template<typename A>
44 : inline double BlockNorm2(const DenseMatrix<A> &mat)
45 : {
46 : double sum=0;
47 0 : for(size_t r=0; r < mat.num_rows(); r++)
48 0 : for(size_t c=0; c < mat.num_cols(); c++)
49 0 : sum += BlockNorm2(mat(r, c));
50 :
51 : return sum;
52 : }
53 :
54 :
55 : template<typename A>
56 : inline double BlockNorm2(const DenseVector<A> &v)
57 : {
58 : double sum=0;
59 0 : for(size_t i=0; i < v.size(); i++)
60 0 : sum += BlockNorm2(v[i]);
61 : return sum;
62 : }
63 :
64 :
65 :
66 : template<typename A>
67 : inline double BlockMaxNorm(const DenseVector<A> &v)
68 : {
69 0 : double max=0;
70 0 : for(size_t i=0; i < v.size(); i++)
71 0 : max = std::max(max, BlockMaxNorm(v[i]));
72 0 : return max;
73 : }
74 :
75 :
76 :
77 :
78 : //////////////////////////////////////////////////////
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
84 : template<typename A, typename B, typename C>
85 : inline 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 :
97 : template<typename A, typename B, typename C>
98 0 : inline 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 0 : for(size_t r=0; r < mA.num_rows(); r++)
103 0 : for(size_t c=0; c < mB.num_cols(); c++)
104 : {
105 : AssignMult(dest(r, c), mA(r, 0), mB(0, c));
106 0 : for(size_t k=1; k < mB.num_rows(); k++)
107 : AddMult(dest(r, c), mA(r, k), mB(k, c));
108 : }
109 0 : }
110 :
111 : // dest += b*vec
112 : template<typename A, typename B, typename C>
113 : inline 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
126 : template<typename A, typename B, typename C>
127 0 : inline 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 0 : for(size_t r=0; r < mA.num_rows(); r++)
132 0 : for(size_t c=0; c < mB.num_cols(); c++)
133 : {
134 0 : for(size_t k=0; k < mB.num_rows(); k++)
135 : AddMult(dest(r, c), mA(r, k), mB(k, c));
136 : }
137 0 : }
138 :
139 :
140 :
141 : // dest -= b*vec
142 : template<typename A, typename B, typename C>
143 : inline 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
158 : template<typename A, typename B>
159 : inline 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 :
169 : template<typename A, typename B>
170 : inline void AddMult(DenseMatrix<A> &dest, const double &alpha, const DenseMatrix<B> &mat)
171 : {
172 : dest.resize(mat.num_rows(), mat.num_cols());
173 0 : for(size_t r=0; r < mat.num_rows(); r++)
174 : {
175 0 : for(size_t c=0; c < mat.num_cols(); c++)
176 : AddMult(dest(r, c), alpha, mat(r, c));
177 : }
178 : }
179 :
180 : template<typename A, typename B>
181 : inline 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
195 : template<typename A, typename B>
196 : inline 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
204 : template<typename A, typename B>
205 : inline 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
213 : template<typename A, typename B>
214 : inline 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
222 : template<typename A> inline void AssignMult(A &dest, const A &vec, const double &b)
223 : {
224 : AssignMult(dest, b, vec);
225 : }
226 : // dest += vec*b
227 : template<typename A> inline void AddMult(A &dest, const A &vec, const double &b)
228 : {
229 : AddMult(dest, b, vec);
230 : }
231 : // dest -= vec*b
232 : template<typename A> inline void SubMult(A &dest, const A &vec, const double &b)
233 : {
234 : SubMult(dest, b, vec);
235 : }
236 :
237 :
238 : template<typename T>
239 : inline 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
245 : template<typename T>
246 : inline void SetSize(DenseVector<T> &t, size_t a)
247 : {
248 : t.resize(a);
249 : }
250 :
251 : // getSize
252 : template<typename T>
253 : inline size_t GetSize(const DenseVector<T> &t)
254 : {
255 : return t.size();
256 : }
257 :
258 : //getRows
259 : template<typename T>
260 : inline size_t GetRows(const DenseMatrix<T> &t)
261 : {
262 : return t.num_rows();
263 : }
264 :
265 : template<typename T>
266 : inline size_t GetCols(const DenseMatrix<T> &t)
267 : {
268 : return t.num_cols();
269 : }
270 :
271 : template<typename T>
272 : struct block_traits;
273 :
274 : template<typename T>
275 : struct block_traits< DenseMatrix<T> >
276 : {
277 : enum { ordering = DenseMatrix<T>::ordering };
278 : enum { is_static = DenseMatrix<T>::is_static};
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 :
286 : template<typename T>
287 : struct block_traits< DenseVector<T> >
288 : {
289 : enum { is_static = DenseVector<T>::is_static};
290 : enum { static_size = DenseVector<T>::static_size};
291 : };
292 :
293 : template<typename T1, typename T2>
294 : struct block_multiply_traits;
295 :
296 : template<typename T1, typename T2>
297 : struct block_multiply_traits<DenseMatrix<T1>, DenseVector<T2> >
298 : {
299 : typedef DenseVector<T2> ReturnType;
300 : };
301 :
302 : //////////////////////////////////////////////////////////////////////////////////////////////
303 : // block_traits
304 :
305 : template<typename T> class DenseMatrixInverse;
306 :
307 : //////////////////////////////////////////////////////////////////////////////////////////////
308 : // variable matrix
309 : template<eMatrixOrdering TOrdering>
310 : struct block_traits< DenseMatrix< VariableArray2<number, 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 :
318 : typedef DenseMatrixInverse< VariableArray2<number, TOrdering> > inverse_type;
319 : };
320 : //////////////////////////////////////////////////////////////////////////////////////////////
321 : // fixed matrix
322 : template<size_t TBlockSize, eMatrixOrdering TOrdering>
323 : struct 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 :
331 : typedef DenseMatrixInverse< FixedArray2<number, TBlockSize, TBlockSize, TOrdering> > inverse_type;
332 : };
333 :
334 : //////////////////////////////////////////////////////////////////////////////////////////////
335 : // variable block matrix
336 : template<typename TValue, eMatrixOrdering TOrdering>
337 : struct 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 :
345 : typedef DenseMatrixInverse< VariableArray2<number, TOrdering> > inverse_type;
346 : };
347 :
348 : template<typename TValue, size_t TBlockSize, eMatrixOrdering TOrdering>
349 : struct 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 :
357 : typedef DenseMatrixInverse< VariableArray2<number, TOrdering> > inverse_type;
358 : };
359 :
360 : //////////////////////////////////////////////////////////////////////////////////////////////
361 : // fixed 1x1 to 3x3 : inverse is matrix
362 : template<eMatrixOrdering TOrdering>
363 : struct 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 :
370 : typedef DenseMatrix< FixedArray2<number, 1, 1, TOrdering> > inverse_type;
371 : };
372 :
373 : template<eMatrixOrdering TOrdering>
374 : struct 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 :
381 : typedef DenseMatrix< FixedArray2<number, 2, 2, TOrdering> > inverse_type;
382 : };
383 :
384 : template<eMatrixOrdering TOrdering>
385 : struct 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 :
392 : typedef DenseMatrix< FixedArray2<number, 3, 3, TOrdering> > inverse_type;
393 : };
394 :
395 :
396 : template<typename T> struct block_multiply_traits<DenseMatrix<T>, DenseMatrix<T> >
397 : {
398 : typedef DenseMatrix<T> ReturnType;
399 : };
400 :
401 :
402 : //////////////////////////////////////////////////////////////////////////////////////////////
403 :
404 :
405 :
406 :
407 : } // namespace ug
408 :
409 : #endif // __H__UG__SMALL_ALGEBRA__BLOCK_DENSE__
|
