ug4
Loading...
Searching...
No Matches
densematrix_inverse.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#ifndef __H__UG__CPU_ALGEBRA__DENSEMATRIX_SMALL_INVERSE_H__
34#define __H__UG__CPU_ALGEBRA__DENSEMATRIX_SMALL_INVERSE_H__
35
36#include "densematrix.h"
37#include "densevector.h"
38#include "block_dense.h"
39#include "common/common.h"
40#include "lib_algebra/small_algebra/small_algebra.h" // for InvertNdyn
41#include <algorithm>
42
43//
44namespace ug {
45
48
50// 1x1
51
52template<typename T>
53inline bool GetInverse1(DenseMatrix<T> &inv, const DenseMatrix<T> &mat)
54{
55 //UG_ASSERT(mat(0,0)!=0.0, "Determinant zero, cannot invert matrix.");
56 if(mat(0,0) == 0.0) return false;
57 inv(0,0) = 1/mat(0,0);
58 return true;
59}
60
61template<typename T>
63{
64 //UG_ASSERT(mat(0,0)!=0.0, "Determinant zero, cannot invert matrix.");
65 if(mat(0,0) == 0.0) return false;
66 mat(0,0) = 1/mat(0,0);
67 return true;
68};
69
71{
72 return GetInverse1(inv, mat);
73}
74
76{
77 return Invert1(mat);
78}
79
80template<typename vector_t, typename matrix_t>
81inline bool InverseMatMult1(DenseVector<vector_t> &dest, double beta1,
83{
84 if(A1(0,0) == 0.0) return false;
85 UG_ASSERT(&dest != &w1, "");
86 dest[0] = beta1*w1[0]/A1(0,0);
87 return true;
88}
89
90inline bool InverseMatMult(DenseVector< FixedArray1<double, 1> > &dest, double beta1,
92{
93 return InverseMatMult1(dest, beta1, A1, w1);
94}
95
97// 2x2
98
99template<typename T>
100inline double GetDet2(const DenseMatrix<T> &mat)
101{
102 UG_ASSERT(mat.num_rows() == 2 && mat.num_cols() == 2, "only for 2x2-matrices");
103 return mat(0,0)*mat(1,1) - mat(1,0)*mat(0,1);
104}
105
106template<typename T>
107inline bool GetInverse2(DenseMatrix<T> &inv, const DenseMatrix<T> &mat)
108{
109 UG_ASSERT(&inv != &mat, "inv and mat have to be different. Otherwise use Invert/Invert2");
110 double invdet = GetDet2(mat);
111 //UG_ASSERT(invdet != 0, "Determinant zero, cannot invert matrix.");
112 if(invdet == 0.0) return false;
113 invdet = 1.0/invdet;
114 inv(0,0) = mat(1,1) * invdet;
115 inv(1,1) = mat(0,0) * invdet;
116 inv(0,1) = mat(0,1) * -invdet;
117 inv(1,0) = mat(1,0) * -invdet;
118 return true;
119}
120
121template<typename T>
123{
124 double invdet = GetDet2(mat);
125 //UG_ASSERT(invdet != 0, "Determinant zero, cannot invert matrix.");
126 if(invdet == 0.0) return false;
127 invdet = 1.0/invdet;
128
129 std::swap(mat(0,0), mat(1,1));
130
131 mat(0,0) *= invdet;
132 mat(0,1) *= -invdet;
133 mat(1,0) *= -invdet;
134 mat(1,1) *= invdet;
135 return true;
136};
137
138
140{
141 return GetInverse2(inv, mat);
142}
143
145{
146 return Invert2(mat);
147}
148
149template<typename vector_t, typename matrix_t>
150inline bool InverseMatMult2(DenseVector<vector_t> &dest, double beta,
151 const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
152{
153 number det = GetDet2(mat);
154 //UG_ASSERT(det != 0, "Determinant zero, cannot invert matrix.");
155 UG_ASSERT(&dest != &vec, "");
156 if(det == 0.0) return false;
157 dest[0] = beta * (mat(1,1)*vec[0] - mat(0,1)*vec[1]) / det;
158 dest[1] = beta * (-mat(1,0)*vec[0] + mat(0,0)*vec[1]) / det;
159 return true;
160}
161
162template<typename T>
163inline bool InverseMatMult(DenseVector< FixedArray1<double, 2> > &dest, double beta,
165{
166 return InverseMatMult2(dest, beta, mat, vec);
167}
168
170// 3x3
171
172
173template<typename T>
174inline double GetDet3(const DenseMatrix<T> &mat)
175{
176 UG_ASSERT(mat.num_rows() == 3 && mat.num_cols() == 3, "only for 3x3-matrices");
177 return mat(0,0)*mat(1,1)*mat(2,2) + mat(0,1)*mat(1,2)*mat(2,0) + mat(0,2)*mat(1,0)*mat(2,1)
178 - mat(0,0)*mat(1,2)*mat(2,1) - mat(0,1)*mat(1,0)*mat(2,2) - mat(0,2)*mat(1,1)*mat(2,0);
179}
180
181template<typename T>
182inline bool GetInverse3(DenseMatrix<T> &inv, const DenseMatrix<T> &mat)
183{
184 UG_ASSERT(&inv != &mat, "inv and mat have to be different. Otherwise use Invert/Invert3");
185 double invdet = GetDet3(mat);
186 //UG_ASSERT(invdet != 0, "Determinant zero, cannot invert matrix.");
187 if(invdet == 0.0) return false;
188 invdet = 1.0/invdet;
189
190 inv(0,0) = ( mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)) * invdet;
191 inv(0,1) = (-mat(0,1)*mat(2,2) + mat(0,2)*mat(2,1)) * invdet;
192 inv(0,2) = ( mat(0,1)*mat(1,2) - mat(0,2)*mat(1,1)) * invdet;
193 inv(1,0) = (-mat(1,0)*mat(2,2) + mat(1,2)*mat(2,0)) * invdet;
194 inv(1,1) = ( mat(0,0)*mat(2,2) - mat(0,2)*mat(2,0)) * invdet;
195 inv(1,2) = (-mat(0,0)*mat(1,2) + mat(0,2)*mat(1,0)) * invdet;
196 inv(2,0) = ( mat(1,0)*mat(2,1) - mat(1,1)*mat(2,0)) * invdet;
197 inv(2,1) = (-mat(0,0)*mat(2,1) + mat(0,1)*mat(2,0)) * invdet;
198 inv(2,2) = ( mat(0,0)*mat(1,1) - mat(0,1)*mat(1,0)) * invdet;
199 return true;
200}
201
203{
205 if(GetInverse3(inv, mat) == false) return false;
206 mat = inv;
207 return true;
208}
209
211{
213 inv.resize(3,3);
214 if(GetInverse3(inv, mat) == false) return false;
215 mat = inv;
216 return true;
217}
218
220{
221 return GetInverse3(inv, mat);
222}
223
225{
226 return Invert3(mat);
227}
228
229template<typename vector_t, typename matrix_t>
230inline bool InverseMatMult3(DenseVector<vector_t> &dest, double beta,
231 const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
232{
233 number det = GetDet3(mat);
234 //UG_ASSERT(det != 0, "Determinant zero, cannot invert matrix.");
235 UG_ASSERT(&dest != &vec, "");
236 if(det == 0.0) return false;
237 dest[0] = ( ( mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)) *vec[0] +
238 (-mat(0,1)*mat(2,2) + mat(0,2)*mat(2,1)) *vec[1] +
239 ( mat(0,1)*mat(1,2) - mat(0,2)*mat(1,1)) *vec[2] ) * beta / det;
240 dest[1] = ( (-mat(1,0)*mat(2,2) + mat(1,2)*mat(2,0)) * vec[0] +
241 ( mat(0,0)*mat(2,2) - mat(0,2)*mat(2,0)) * vec[1] +
242 (-mat(0,0)*mat(1,2) + mat(0,2)*mat(1,0)) * vec[2] ) * beta / det;
243 dest[2] = ( ( mat(1,0)*mat(2,1) - mat(1,1)*mat(2,0)) * vec[0] +
244 (-mat(0,0)*mat(2,1) + mat(0,1)*mat(2,0)) * vec[1] +
245 ( mat(0,0)*mat(1,1) - mat(0,1)*mat(1,0)) * vec[2] ) * beta / det;
246 return true;
247}
248
249template<typename T>
250inline bool InverseMatMult(DenseVector< FixedArray1<double, 3> > &dest, double beta,
252{
253 return InverseMatMult3(dest, beta, mat, vec);
254}
256
257
258
260template<typename vector_t, typename matrix_t>
262 const number &beta1, const DenseMatrixInverse<matrix_t> &A1, const DenseVector<vector_t> &w1)
263{
264 if(beta1 == 1.0)
265 {
266 dest = w1;
267 A1.apply(dest);
268 }
269 else
270 {
272 tmp = w1;
273 A1.apply(tmp);
274 VecScaleAssign(dest, beta1, dest);
275 }
276}
277
278
280template<typename vector_t, typename matrix_t>
282 const number &alpha1, const DenseVector<vector_t> &v1,
283 const number &beta1, const DenseMatrixInverse<matrix_t> &A1, const DenseVector<vector_t> &w1)
284{
285 // todo: use dynamic stack here
287 tmp = w1;
288 A1.apply(tmp);
289 VecScaleAdd(dest, alpha1, v1, beta1, tmp);
290}
291
292
293
294template<typename T, eMatrixOrdering TOrdering>
296{
297 return inv.set_as_inverse_of(mat);
298}
299
300template<typename T, size_t TBlockSize, eMatrixOrdering TOrdering>
305
307template<typename T>
308inline bool Invert(DenseMatrix<T> &mat)
309{
310 switch(mat.num_rows())
311 {
312 case 1: return Invert1(mat);
313 case 2: return Invert2(mat);
314 case 3: return Invert3(mat);
315 default: return InvertNdyn(mat);
316 }
317}
318
319template<typename vector_t, typename matrix_t>
320inline bool InverseMatMultN(DenseVector<vector_t> &dest, double beta,
321 const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
322{
323 typename block_traits<DenseMatrix<matrix_t> >::inverse_type inv;
324 if(!GetInverse(inv, mat)) return false;
325 MatMult(dest, beta, inv, vec);
326 return true;
327}
328
329
330template<typename vector_t, typename matrix_t>
331inline bool InverseMatMult(DenseVector<vector_t> &dest, double beta,
332 const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
333{
334 switch(mat.num_rows())
335 {
336 case 1: return InverseMatMult1(dest, beta, mat, vec);
337 case 2: return InverseMatMult2(dest, beta, mat, vec);
338 case 3: return InverseMatMult3(dest, beta, mat, vec);
339 default: return InverseMatMultN(dest, beta, mat, vec);
340 }
341}
342
344
345// end group small_algebra
347
348}
349
350#endif // __H__UG__CPU_ALGEBRA__DENSEMATRIX_SMALL_INVERSE_H__
Definition densematrix.h:57
Definition lapack_densematrix_inverse.h:51
void apply(DenseVector< vector_t > &vec) const
Definition lapack_densematrix_inverse.h:115
Definition densevector.h:101
Definition fixed_array.h:56
Definition fixed_array.h:135
Definition variable_array.h:139
bool GetInverse3(DenseMatrix< T > &inv, const DenseMatrix< T > &mat)
Definition densematrix_inverse.h:182
bool InverseMatMultN(DenseVector< vector_t > &dest, double beta, const DenseMatrix< matrix_t > &mat, const DenseVector< vector_t > &vec)
Definition densematrix_inverse.h:320
bool InverseMatMult3(DenseVector< vector_t > &dest, double beta, const DenseMatrix< matrix_t > &mat, const DenseVector< vector_t > &vec)
Definition densematrix_inverse.h:230
bool GetInverse2(DenseMatrix< T > &inv, const DenseMatrix< T > &mat)
Definition densematrix_inverse.h:107
bool InverseMatMult1(DenseVector< vector_t > &dest, double beta1, const DenseMatrix< matrix_t > &A1, const DenseVector< vector_t > &w1)
Definition densematrix_inverse.h:81
double GetDet2(const DenseMatrix< T > &mat)
Definition densematrix_inverse.h:100
bool InverseMatMult2(DenseVector< vector_t > &dest, double beta, const DenseMatrix< matrix_t > &mat, const DenseVector< vector_t > &vec)
Definition densematrix_inverse.h:150
double GetDet3(const DenseMatrix< T > &mat)
Definition densematrix_inverse.h:174
bool GetInverse1(DenseMatrix< T > &inv, const DenseMatrix< T > &mat)
Definition densematrix_inverse.h:53
bool Invert1(DenseMatrix< T > &mat)
Definition densematrix_inverse.h:62
bool Invert3(DenseMatrix< FixedArray2< double, 3, 3 > > &mat)
Definition densematrix_inverse.h:202
bool Invert2(DenseMatrix< T > &mat)
Definition densematrix_inverse.h:122
#define UG_ASSERT(expr, msg)
Definition assert.h:70
double number
Definition types.h:124
void VecScaleAdd(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1, typename vector_t::value_type s2, const vector_t &v2)
Scales two Vectors, adds them and returns the sum in a third vector.
Definition math_vector_functions_common_impl.hpp:265
the ug namespace
bool GetInverse(block_traits< T >::inverse_type &inv, const T &m)
bool Invert(T &m)
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:68
bool InvertNdyn(DenseMatrix< T > &mat)
Definition lapack_invert.h:47
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
void VecScaleAssign(double &dest, double alpha1, const double &v1)
calculates dest = alpha1*v1. for doubles
Definition operations_vec.h:49
Definition communication_policies.h:58