Line data Source code
1 : /*
2 : * Copyright (c) 2010-2014: 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 : * \file densematrix_impl.h
35 : *
36 : * \author Martin Rupp
37 : *
38 : * \date 21.07.2010
39 : *
40 : * Goethe-Center for Scientific Computing 2010.
41 : *
42 : * awesome!
43 : */
44 :
45 : #ifndef __H__UG__COMMON__DENSEMATRIX_IMPL_H__
46 : #define __H__UG__COMMON__DENSEMATRIX_IMPL_H__
47 :
48 : #include "print.h"
49 : #include "../blocks.h"
50 :
51 : // constructors
52 : namespace ug{
53 :
54 :
55 : template<typename TStorage>
56 : DenseMatrix<TStorage>::DenseMatrix() : TStorage()
57 : {
58 : }
59 :
60 : template<typename TStorage>
61 0 : DenseMatrix<TStorage>::DenseMatrix(const DenseMatrix<TStorage> &rhs) : TStorage(rhs)
62 : {
63 : }
64 :
65 : template<typename TStorage>
66 : DenseMatrix<TStorage>::DenseMatrix(double value) : TStorage()
67 : {
68 : operator =(value);
69 : }
70 :
71 : // matrix assignment operators
72 : ////// =
73 :
74 : template<typename TStorage>
75 : DenseMatrix<TStorage> &
76 0 : DenseMatrix<TStorage>::operator =(const this_type &t)
77 : {
78 0 : resize(t.num_rows(), t.num_cols());
79 0 : for(size_t r1=0; r1<t.num_rows(); r1++)
80 0 : for(size_t c1=0; c1<t.num_cols(); c1++)
81 0 : entry(r1, c1) = t(r1, c1);
82 0 : return *this;
83 : }
84 :
85 : template<typename TStorage>
86 : template<typename T2>
87 : DenseMatrix<TStorage> &
88 : DenseMatrix<TStorage>::operator =(const T2 &t)
89 : {
90 : resize(t.num_rows(), t.num_cols());
91 : for(size_t r1=0; r1<t.num_rows(); r1++)
92 : for(size_t c1=0; c1<t.num_cols(); c1++)
93 : entry(r1, c1) = t(r1, c1);
94 : return *this;
95 : }
96 :
97 : template<typename TStorage>
98 : DenseMatrix<TStorage> &
99 : // this is stupid since i like to make it that way that i have one ::value_type and one double,
100 : // but then there are two (double) functions...
101 : //DenseMatrix<TStorage>::operator = (const typename DenseMatrix<TStorage>::value_type &rhs)
102 0 : DenseMatrix<TStorage>::operator = (double rhs)
103 : {
104 0 : for(size_t r=0; r<num_rows(); ++r)
105 0 : for(size_t c=0; c<num_cols(); ++c)
106 : {
107 0 : if(r==c) entry(r, c) = rhs;
108 0 : else entry(r, c) = 0.0;
109 : }
110 0 : return *this;
111 : }
112 :
113 : ////// +=
114 : template<typename TStorage>
115 : template<typename T2>
116 : DenseMatrix<TStorage> &
117 : DenseMatrix<TStorage>::operator += (const T2 &t)
118 : {
119 : UG_ASSERT(t.num_rows() == num_rows() && t.num_cols() == num_cols(), "");
120 0 : for(size_t r1=0; r1<t.num_rows(); r1++)
121 0 : for(size_t c1=0; c1<t.num_cols(); c1++)
122 0 : entry(r1, c1) += t(r1, c1);
123 : return *this;
124 : }
125 :
126 : template<typename TStorage>
127 : DenseMatrix<TStorage> &
128 : DenseMatrix<TStorage>::operator += (double alpha)
129 : {
130 : size_t minimum=num_rows() > num_cols() ? num_cols() : num_rows();
131 : for(size_t i=0; i<minimum; i++)
132 : entry(i, i) += alpha;
133 : return *this;
134 : }
135 :
136 : ////// -=
137 :
138 : template<typename TStorage>
139 : template<typename T2>
140 : DenseMatrix<TStorage> &
141 : DenseMatrix<TStorage>::operator -= (const T2 &t)
142 : {
143 : UG_ASSERT(t.num_rows() == num_rows() && t.num_cols() == num_cols(), "");
144 0 : for(size_t r1=0; r1<t.num_rows(); r1++)
145 0 : for(size_t c1=0; c1<t.num_cols(); c1++)
146 0 : entry(r1, c1) -= t(r1, c1);
147 : return *this;
148 : }
149 :
150 : template<typename TStorage>
151 : DenseMatrix<TStorage> &
152 : DenseMatrix<TStorage>::operator -= (double alpha)
153 : {
154 : for(size_t r=0; r<num_rows(); ++r)
155 : for(size_t c=0; c<num_cols(); ++c)
156 : entry(r, c) -= alpha;
157 : return *this;
158 : }
159 :
160 : ////// *=
161 :
162 : template<typename TStorage>
163 : DenseMatrix<TStorage> &
164 : DenseMatrix<TStorage>::operator *= (double alpha)
165 : {
166 0 : for(size_t r=0; r<num_rows(); ++r)
167 0 : for(size_t c=0; c<num_cols(); ++c)
168 0 : entry(r, c) *= alpha;
169 : return *this;
170 : }
171 :
172 : template<typename TStorage>
173 : DenseMatrix<TStorage> &
174 : DenseMatrix<TStorage>::operator *= (const this_type &mat)
175 : {
176 : operator=(operator*(mat));
177 : return *this;
178 : }
179 :
180 : ////// /=
181 : template<typename TStorage>
182 : DenseMatrix<TStorage> &
183 : DenseMatrix<TStorage>::operator /= (double alpha)
184 : {
185 : for(size_t r=0; r<num_rows(); ++r)
186 : for(size_t c=0; c<num_cols(); ++c)
187 : entry(r, c) /= alpha;
188 : return *this;
189 : }
190 :
191 : template<typename TStorage>
192 : DenseMatrix<TStorage> &
193 0 : DenseMatrix<TStorage>::operator /= (this_type &other)
194 : {
195 : this_type tmp = other;
196 : bool success = Invert(tmp);
197 0 : UG_COND_THROW(!success, "Failed to invert dense matrix.");
198 0 : (*this) = (*this) * tmp;
199 0 : return *this;
200 : }
201 :
202 :
203 : ////// +
204 : template<typename TStorage>
205 : DenseMatrix<TStorage>
206 : DenseMatrix<TStorage>::operator + (const this_type &other ) const
207 : {
208 : UG_ASSERT(num_rows() == other.num_rows() && num_cols() == other.num_cols(), "");
209 : this_type erg;
210 : erg.resize(num_rows(), num_cols());
211 0 : for(size_t r=0; r<num_rows(); r++)
212 0 : for(size_t c=0; c<num_cols(); c++)
213 0 : erg(r, c) = entry(r, c) + other(r,c);
214 : return erg;
215 : }
216 : ////// -
217 : template<typename TStorage>
218 : DenseMatrix<TStorage>
219 : DenseMatrix<TStorage>::operator - (const this_type &other ) const
220 : {
221 : UG_ASSERT(num_rows() == other.num_rows() && num_cols() == other.num_cols(), "");
222 : this_type erg;
223 : erg.resize(num_rows(), num_cols());
224 :
225 : for(size_t r=0; r<num_rows(); r++)
226 : for(size_t c=0; c<num_cols(); c++)
227 : erg(r, c) = entry(r, c) - other(r,c);
228 : return erg;
229 : }
230 :
231 : ////// unary -
232 : template<typename TStorage>
233 : DenseMatrix<TStorage>
234 : DenseMatrix<TStorage>::operator - () const
235 : {
236 : this_type erg;
237 : erg.resize(num_rows(), num_cols());
238 0 : for(size_t r=0; r < num_rows(); r++)
239 0 : for(size_t c=0; c < num_cols(); c++)
240 : {
241 0 : erg(r,c) = entry(r, c);
242 0 : erg(r,c) *= -1.0;
243 : }
244 : return erg;
245 : }
246 :
247 : // multiply
248 : ////// *
249 : template<typename TStorage>
250 : DenseMatrix<TStorage>
251 0 : DenseMatrix<TStorage>::operator * (const this_type &other ) const
252 : {
253 : // that aint 100% correct
254 : UG_ASSERT(num_cols() == other.num_rows(), "");
255 :
256 : this_type erg;
257 : erg.resize(num_rows(), other.num_cols());
258 :
259 0 : for(size_t r=0; r < num_rows(); r++)
260 0 : for(size_t c=0; c < other.num_cols(); c++)
261 : {
262 0 : erg(r,c) = 0.0;
263 0 : for(size_t i=0; i < num_cols(); i++)
264 : AddMult(erg(r,c), at(r, i), other.at(i, c));
265 : }
266 0 : return erg;
267 : }
268 :
269 :
270 : template<typename TStorage>
271 : DenseMatrix<TStorage>
272 : DenseMatrix<TStorage>::T() const
273 : {
274 : this_type erg;
275 : erg.resize(num_rows(), num_cols());
276 : for(size_t r=0; r < num_rows(); r++)
277 : for(size_t c=0; c < num_cols(); c++)
278 : erg(r,c) = entry(c, r);
279 : return erg;
280 : }
281 :
282 : template<typename TStorage>
283 : template<typename TStorage2>
284 : DenseVector<TStorage2>
285 0 : DenseMatrix<TStorage>::operator * (const DenseVector<TStorage2> &vec) const
286 : {
287 : UG_ASSERT(num_cols() == vec.size(), "");
288 : DenseVector<TStorage2> erg;
289 0 : erg.resize(num_rows());
290 :
291 0 : for(size_t r=0; r < num_rows(); r++)
292 : {
293 0 : erg[r] = 0.0;
294 0 : for(size_t c=0; c < num_cols(); c++)
295 0 : erg[r] += at(r,c) * vec[c];
296 : }
297 0 : return erg;
298 : }
299 :
300 : template<typename TStorage>
301 : DenseMatrix<TStorage>
302 : DenseMatrix<TStorage>::operator * (double alpha ) const
303 : {
304 : this_type erg;
305 : erg.resize(num_rows(), num_cols());
306 :
307 0 : for(size_t r=0; r < num_rows(); r++)
308 0 : for(size_t c=0; c < num_cols(); c++)
309 0 : erg(r,c) = at(r,c)*alpha;
310 : return erg;
311 : }
312 :
313 :
314 :
315 : ///// /
316 : template<typename TStorage>
317 : DenseMatrix<TStorage>
318 0 : DenseMatrix<TStorage>::operator / (this_type &other)
319 : {
320 : this_type tmp = other;
321 : Invert(tmp);
322 :
323 0 : return (*this) * tmp;
324 : }
325 :
326 : // compare operators
327 :
328 : template<typename TStorage>
329 : bool
330 0 : DenseMatrix<TStorage>::operator == (double t) const
331 : {
332 0 : for(size_t r=0; r<num_rows(); ++r)
333 0 : for(size_t c=0; c<num_cols(); ++c)
334 : {
335 0 : if(r==c)
336 : {
337 0 : if(entry(r,c) != t) return false;
338 : }
339 : else
340 0 : if(entry(r,c) != 0.0) return false;
341 : }
342 : return true;
343 : }
344 :
345 : template<typename TStorage>
346 : template<typename TStorage2>
347 : bool
348 : DenseMatrix<TStorage>::operator == (const DenseMatrix<TStorage2> &t) const
349 : {
350 : for(size_t r=0; r<num_rows(); ++r)
351 : for(size_t c=0; c<num_cols(); ++c)
352 : if(entry(r,c) != t(r,c)) return false;
353 : return true;
354 : }
355 :
356 :
357 : template<typename TStorage>
358 : template<typename T2>
359 : bool DenseMatrix<TStorage>::operator != (const T2 &t) const
360 : {
361 0 : return !(operator == (t));
362 : }
363 :
364 : template<typename TStorage>
365 : void DenseMatrix<TStorage>::maple_print(const char *name)
366 : {
367 : UG_LOG(MatlabString(*this, name));
368 : }
369 :
370 :
371 : template<typename TStorage>
372 0 : std::ostream &operator << (std::ostream &out, const DenseMatrix<TStorage> &mat)
373 : {
374 0 : out << "[ ";
375 : typedef size_t size_type;
376 0 : for(size_type r=0; r<mat.num_rows(); ++r)
377 : {
378 0 : for(size_type c=0; c<mat.num_cols(); ++c)
379 0 : out << mat(r, c) << " ";
380 0 : if(r != mat.num_rows()-1) out << "| ";
381 : }
382 0 : out << "]";
383 : // out << "(DenseMatrix " << mat.num_rows() << "x" << mat.num_cols() << ", " << ((DenseMatrix<TStorage>::ordering == ColMajor) ? "ColMajor)" : "RowMajor)");
384 :
385 0 : return out;
386 : }
387 :
388 :
389 : template<size_t Tr, size_t Tc>
390 : inline
391 : bool
392 : BlockSerialize(const DenseMatrix<FixedArray2<number, Tr, Tc> > &mat, std::ostream &buff)
393 : {
394 : buff.write((char*)&mat, sizeof(mat));
395 : return true;
396 : }
397 :
398 : template<size_t Tr, size_t Tc>
399 : inline
400 : bool
401 : BlockDeserialize(std::istream &buff, const DenseMatrix<FixedArray2<number, Tr, Tc> > &mat)
402 : {
403 : buff.read((char*)&mat, sizeof(mat));
404 : return true;
405 : }
406 :
407 :
408 : template<typename T>
409 : inline
410 : void
411 : Serialize(std::ostream &buff, const DenseMatrix<VariableArray2<T> > &mat)
412 : {
413 : size_t rows = mat.num_rows();
414 : size_t cols = mat.num_cols();
415 : buff.write((char*)&rows, sizeof(rows));
416 : buff.write((char*)&cols, sizeof(cols));
417 : for(size_t r=0; r<rows; r++)
418 : for(size_t c=0; c<cols; c++)
419 : BlockSerialize(mat(r, c), buff);
420 : }
421 :
422 :
423 : template<typename T>
424 : inline
425 : void
426 : Deserialize(std::istream &buff, const DenseMatrix<VariableArray2<T> > &mat)
427 : {
428 : size_t rows, cols;
429 : buff.read((char*)&rows, sizeof(rows));
430 : buff.read((char*)&cols, sizeof(cols));
431 : mat.resize(rows, cols);
432 : for(size_t r=0; r<rows; r++)
433 : for(size_t c=0; c<cols; c++)
434 : BlockDeserialize(buff, mat(r, c));
435 : }
436 :
437 : template<typename T >
438 : inline bool IsFiniteAndNotTooBig(const DenseMatrix<T> &m)
439 : {
440 : for(size_t r=0; r<m.num_rows(); r++)
441 : for(size_t c=0; c<m.num_cols(); c++)
442 : if(IsFiniteAndNotTooBig(m(r, c)) == false) return false;
443 :
444 : return true;
445 : }
446 : } // namespace ug
447 :
448 : #endif // __H__UG__COMMON__DENSEMATRIX_IMPL_H__
|
