ug4
densematrix_impl.h
Go to the documentation of this file.
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 
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>
57 {
58 }
59 
60 template<typename TStorage>
62 {
63 }
64 
65 template<typename TStorage>
66 DenseMatrix<TStorage>::DenseMatrix(double value) : TStorage()
67 {
68  operator =(value);
69 }
70 
71 // matrix assignment operators
73 
74 template<typename TStorage>
77 {
79  for(size_t r1=0; r1<t.num_rows(); r1++)
80  for(size_t c1=0; c1<t.num_cols(); c1++)
81  entry(r1, c1) = t(r1, c1);
82  return *this;
83 }
84 
85 template<typename TStorage>
86 template<typename T2>
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>
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)
103 {
104  for(size_t r=0; r<num_rows(); ++r)
105  for(size_t c=0; c<num_cols(); ++c)
106  {
107  if(r==c) entry(r, c) = rhs;
108  else entry(r, c) = 0.0;
109  }
110  return *this;
111 }
112 
114 template<typename TStorage>
115 template<typename T2>
116 DenseMatrix<TStorage> &
118 {
119  UG_ASSERT(t.num_rows() == num_rows() && t.num_cols() == num_cols(), "");
120  for(size_t r1=0; r1<t.num_rows(); r1++)
121  for(size_t c1=0; c1<t.num_cols(); c1++)
122  entry(r1, c1) += t(r1, c1);
123  return *this;
124 }
125 
126 template<typename TStorage>
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 
137 
138 template<typename TStorage>
139 template<typename T2>
142 {
143  UG_ASSERT(t.num_rows() == num_rows() && t.num_cols() == num_cols(), "");
144  for(size_t r1=0; r1<t.num_rows(); r1++)
145  for(size_t c1=0; c1<t.num_cols(); c1++)
146  entry(r1, c1) -= t(r1, c1);
147  return *this;
148 }
149 
150 template<typename TStorage>
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 
161 
162 template<typename TStorage>
165 {
166  for(size_t r=0; r<num_rows(); ++r)
167  for(size_t c=0; c<num_cols(); ++c)
168  entry(r, c) *= alpha;
169  return *this;
170 }
171 
172 template<typename TStorage>
175 {
176  operator=(operator*(mat));
177  return *this;
178 }
179 
181 template<typename TStorage>
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>
194 {
195  this_type tmp = other;
196  bool success = Invert(tmp);
197  UG_COND_THROW(!success, "Failed to invert dense matrix.");
198  (*this) = (*this) * tmp;
199  return *this;
200 }
201 
202 
204 template<typename TStorage>
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  for(size_t r=0; r<num_rows(); r++)
212  for(size_t c=0; c<num_cols(); c++)
213  erg(r, c) = entry(r, c) + other(r,c);
214  return erg;
215 }
217 template<typename TStorage>
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 
232 template<typename TStorage>
235 {
236  this_type erg;
237  erg.resize(num_rows(), num_cols());
238  for(size_t r=0; r < num_rows(); r++)
239  for(size_t c=0; c < num_cols(); c++)
240  {
241  erg(r,c) = entry(r, c);
242  erg(r,c) *= -1.0;
243  }
244  return erg;
245 }
246 
247 // multiply
249 template<typename TStorage>
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  for(size_t r=0; r < num_rows(); r++)
260  for(size_t c=0; c < other.num_cols(); c++)
261  {
262  erg(r,c) = 0.0;
263  for(size_t i=0; i < num_cols(); i++)
264  AddMult(erg(r,c), at(r, i), other.at(i, c));
265  }
266  return erg;
267 }
268 
269 
270 template<typename TStorage>
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>
286 {
287  UG_ASSERT(num_cols() == vec.size(), "");
289  erg.resize(num_rows());
290 
291  for(size_t r=0; r < num_rows(); r++)
292  {
293  erg[r] = 0.0;
294  for(size_t c=0; c < num_cols(); c++)
295  erg[r] += at(r,c) * vec[c];
296  }
297  return erg;
298 }
299 
300 template<typename TStorage>
303 {
304  this_type erg;
305  erg.resize(num_rows(), num_cols());
306 
307  for(size_t r=0; r < num_rows(); r++)
308  for(size_t c=0; c < num_cols(); c++)
309  erg(r,c) = at(r,c)*alpha;
310  return erg;
311 }
312 
313 
314 
316 template<typename TStorage>
319 {
320  this_type tmp = other;
321  Invert(tmp);
322 
323  return (*this) * tmp;
324 }
325 
326 // compare operators
327 
328 template<typename TStorage>
329 bool
331 {
332  for(size_t r=0; r<num_rows(); ++r)
333  for(size_t c=0; c<num_cols(); ++c)
334  {
335  if(r==c)
336  {
337  if(entry(r,c) != t) return false;
338  }
339  else
340  if(entry(r,c) != 0.0) return false;
341  }
342  return true;
343 }
344 
345 template<typename TStorage>
346 template<typename TStorage2>
347 bool
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  return !(operator == (t));
362 }
363 
364 template<typename TStorage>
366 {
367  UG_LOG(MatlabString(*this, name));
368 }
369 
370 
371 template<typename TStorage>
372 std::ostream &operator << (std::ostream &out, const DenseMatrix<TStorage> &mat)
373 {
374  out << "[ ";
375  typedef size_t size_type;
376  for(size_type r=0; r<mat.num_rows(); ++r)
377  {
378  for(size_type c=0; c<mat.num_cols(); ++c)
379  out << mat(r, c) << " ";
380  if(r != mat.num_rows()-1) out << "| ";
381  }
382  out << "]";
383 // out << "(DenseMatrix " << mat.num_rows() << "x" << mat.num_cols() << ", " << ((DenseMatrix<TStorage>::ordering == ColMajor) ? "ColMajor)" : "RowMajor)");
384 
385  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
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 >
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__
location name
Definition: checkpoint_util.lua:128
Definition: densematrix.h:57
this_type & operator+=(const T &t)
bool operator!=(const TStorage2 &t) const
this_type & operator-=(const T &t)
this_type & operator*=(double alpha)
Definition: densematrix_impl.h:164
this_type & operator/=(double alpha)
Definition: densematrix_impl.h:183
void maple_print(const char *name)
Definition: densematrix_impl.h:365
bool operator==(double t) const
Definition: densematrix_impl.h:330
this_type operator/(this_type &other)
Definition: densematrix_impl.h:318
DenseMatrix()
Definition: densematrix_impl.h:56
DenseVector< TStorage2 > operator*(const DenseVector< TStorage2 > &vec) const
Definition: densematrix_impl.h:285
this_type T() const
Definition: densematrix_impl.h:272
this_type operator-() const
Definition: densematrix_impl.h:234
this_type & operator=(const this_type &t)
Definition: densematrix_impl.h:76
this_type operator+(const this_type &other) const
Definition: densematrix_impl.h:206
Definition: densevector.h:101
Definition: fixed_array.h:135
Definition: variable_array.h:139
size_type num_rows() const
Definition: variable_array_impl.h:222
size_type num_cols() const
Definition: variable_array_impl.h:229
std::ostream & operator<<(std::ostream &outStream, const ug::MathMatrix< 2, 2 > &m)
Definition: math_matrix.cpp:38
number alpha
#define UG_ASSERT(expr, msg)
Definition: assert.h:70
#define UG_LOG(msg)
Definition: log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition: error.h:61
bool operator==(const MathVector< N, T > &v, const MathVector< N, T > &w)
Definition: math_vector.h:523
the ug namespace
bool BlockSerialize(const DenseMatrix< FixedArray2< number, Tr, Tc > > &mat, std::ostream &buff)
Definition: densematrix_impl.h:392
std::string MatlabString(const DenseMatrix< TStorage > &A, const char *name)
Definition: print.h:65
bool Invert(T &m)
bool BlockDeserialize(std::istream &buff, const DenseMatrix< FixedArray2< number, Tr, Tc > > &mat)
Definition: densematrix_impl.h:401
void Deserialize(TIStream &buf, ParallelVector< T > &v)
Deerialize for ParallelVector<T>
Definition: restart_bridge.cpp:112
bool IsFiniteAndNotTooBig(double d)
Definition: number_util.h:39
void Serialize(TOStream &buf, const ParallelVector< T > &v)
Serialize for ParallelVector<T>
Definition: restart_bridge.cpp:103
void AddMult(A &dest, const B &b, const C &vec)
bool resize(size_t newRows, size_t newCols)
size_t num_rows() const
Definition: sparsematrix_interface.h:38
size_t num_cols() const
Definition: sparsematrix_interface.h:39