ug4
Loading...
Searching...
No Matches
additional_math.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: 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#ifndef FAMG_ADDITIONAL_MATH_H_
34#define FAMG_ADDITIONAL_MATH_H_
36namespace ug{
37
38inline void vecSum(double &erg, double alpha, double vec)
39{
40 erg = alpha * vec;
41}
42
43template<typename T>
44inline void vecSum(typename T::value_type &erg, double alpha, const T &vec)
45{
46 erg = alpha * vec[0];
47 for (size_t i = 1; i < vec.size(); i++)
48 erg += vec[i];
49 erg *= alpha;
50}
51
52inline double vecSum(double alpha, double vec)
53{
54 return alpha * vec;
55}
56template<typename T>
57inline typename T::value_type vecSum(double alpha, const T &vec)
58{
59 typename T::value_type erg;
60 vecSum(erg, alpha, vec);
61 return erg;
62}
63
64inline void matSum(double &erg, double alpha, double vec)
65{
66 erg = alpha * vec;
67}
68
69template<typename T1, typename T2>
70inline void matSum(T1 &erg, double alpha, T2 &mat)
71{
72 for (size_t r = 0; mat.num_rows(); r++)
73 {
74 erg[r] = mat(r, 0);
75 for (size_t c = 1; c < mat.num_cols(); c++)
76 erg[r] += mat(r, c);
77 erg[r] *= alpha;
78 }
79}
80
81template<typename T1, typename T2>
82inline T1 matSum(double alpha, T2 &mat)
83{
84 T1 erg;
85 matSum(erg, alpha, mat);
86 return erg;
87}
88
89template<typename T1>
91{
92 typename DenseMatrix<T1>::value_type ret = 0.0;
93 for (size_t r = 0; mat.num_rows(); r++)
94 for (size_t c = 1; c < mat.num_cols(); c++)
95 ret += mat(r, c);
96 return ret;
97}
98
99inline double matTrace(const double d)
100{
101 return d;
102}
103
104template<typename T1>
105inline double matTrace(const DenseMatrix<T1> &mat)
106{
107 double tr=0.0;
108 const size_t rk = (mat.num_rows() < mat.num_cols()) ? mat.num_rows() : mat.num_cols();
109 for (size_t k = 0; k<rk; k++)
110 {
111 tr += mat(k, k);
112 }
113 return tr;
114}
115
116inline double matDiagMax(const double d)
117{
118 return d;
119}
120
121template<typename T1>
122inline double matDiagMax(const DenseMatrix<T1> &mat)
123{
124 double val=0.0;
125 double max=0.0;
126 const size_t rk = (mat.num_rows() < mat.num_cols()) ? mat.num_rows() : mat.num_cols();
127 for (size_t k = 0; k<rk; k++)
128 {
129 double abs=fabs(mat(k, k));
130 if (abs>max) {
131 val = mat(k, k);
132 max = abs;
133 }
134 }
135 return val;
136}
137
138
139inline double Sum1Mat1(double d)
140{
141 return d;
142}
143
144inline void GetDiag(double &a, double b)
145{
146 a = b;
147}
148
149template<typename T1, typename T2>
150inline void GetDiag(T1 &m1, const T2 &m)
151{
152 UG_ASSERT(m.num_rows()==m.num_cols(), "");
153 m1.resize(m.num_rows(), m.num_rows());
154 m1=0.0;
155 for (size_t i = 0; i < m.num_rows(); i++)
156 m1(i, i) = m(i, i);
157}
158
159template<typename T1, typename T2>
160inline void GetDiagSqrt(T1 &v, const T2 &m)
161{
162 UG_ASSERT(m.num_rows()==m.num_cols(), "");
163 v.resize(m.num_rows());
164 for (size_t i = 0; i < m.num_rows(); i++)
165 v[i] = sqrt(m(i, i));
166}
167
168inline void GetDiagSqrt(double &a, double b)
169{
170 a = sqrt(b);
171}
172
173inline double EnergyProd(double v1, double M, double v2)
174{
175 return v1 * M * v2;
176}
177
178template<typename T1, typename T2>
179inline double EnergyProd(const T1 &v1, const DenseMatrix<T2> &M, const T1 &v2)
180{
181 double sum = 0;
182 for (size_t r = 0; r < M.num_rows(); r++)
183 {
184 double t = 0;
185 for (size_t c = 0; c < M.num_cols(); c++)
186 t += M(r, c) * v2[c];
187 sum += t * v1[r];
188 }
189 return sum;
190}
191
192template<typename TMatrix>
194{
196 Ad.resize(blockSize*Ab.num_rows(), blockSize*Ab.num_cols());
197
198 for(size_t i=0; i<Ab.num_rows(); i++)
199 {
200 for(size_t j=0; j<Ab.num_cols(); j++)
201 for(size_t r=0; r<blockSize; r++)
202 for(size_t c=0; c<blockSize; c++)
203 Ad(i*blockSize+r, j*blockSize+c) = BlockRef(Ab(i, j), r, c);
204 }
205}
206
207}
208#endif /* ADDITIONAL_MATH_H_ */
Definition densematrix.h:57
TStorage::value_type value_type
Definition densematrix.h:59
Definition variable_array.h:139
#define UG_ASSERT(expr, msg)
Definition assert.h:70
the ug namespace
void matSum(double &erg, double alpha, double vec)
Definition additional_math.h:64
void GetDiagSqrt(T1 &v, const T2 &m)
Definition additional_math.h:160
double matDiagMax(const double d)
Definition additional_math.h:116
double matTrace(const double d)
Definition additional_math.h:99
double EnergyProd(vector_type &v1, matrix_type &A, vector_type &v2, vector_type &tmp)
Definition additional_math.h:90
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
void vecSum(double &erg, double alpha, double vec)
Definition additional_math.h:38
void BlockMatrixToDoubleMatrix(DenseMatrix< VariableArray2< double > > &Ad, TMatrix &Ab)
Definition additional_math.h:193
void GetDiag(double &a, double b)
Definition additional_math.h:144
DenseMatrix< T1 >::value_type Sum1Mat1(const DenseMatrix< T1 > &mat)
Definition additional_math.h:90
Definition communication_policies.h:58