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 : #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 : //
44 : namespace ug {
45 :
46 : /// \addtogroup small_algebra
47 : /// \{
48 :
49 : //////////////////////////////////////////////////////
50 : // 1x1
51 :
52 : template<typename T>
53 : inline 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 :
61 : template<typename T>
62 : bool Invert1(DenseMatrix<T> &mat)
63 : {
64 : //UG_ASSERT(mat(0,0)!=0.0, "Determinant zero, cannot invert matrix.");
65 0 : if(mat(0,0) == 0.0) return false;
66 0 : mat(0,0) = 1/mat(0,0);
67 0 : return true;
68 : };
69 :
70 : inline bool GetInverse(DenseMatrix<FixedArray2<double, 1, 1> > &inv, const DenseMatrix<FixedArray2<double, 1, 1> > &mat)
71 : {
72 : return GetInverse1(inv, mat);
73 : }
74 :
75 : inline bool Invert(DenseMatrix< FixedArray2<double, 1, 1> > &mat)
76 : {
77 : return Invert1(mat);
78 : }
79 :
80 : template<typename vector_t, typename matrix_t>
81 : inline bool InverseMatMult1(DenseVector<vector_t> &dest, double beta1,
82 : const DenseMatrix<matrix_t> &A1, const DenseVector<vector_t> &w1)
83 : {
84 0 : if(A1(0,0) == 0.0) return false;
85 : UG_ASSERT(&dest != &w1, "");
86 0 : dest[0] = beta1*w1[0]/A1(0,0);
87 0 : return true;
88 : }
89 :
90 : inline bool InverseMatMult(DenseVector< FixedArray1<double, 1> > &dest, double beta1,
91 : const DenseMatrix< FixedArray2<double, 1, 1> > &A1, const DenseVector< FixedArray1<double, 1> > &w1)
92 : {
93 : return InverseMatMult1(dest, beta1, A1, w1);
94 : }
95 :
96 : //////////////////////
97 : // 2x2
98 :
99 : template<typename T>
100 : inline double GetDet2(const DenseMatrix<T> &mat)
101 : {
102 : UG_ASSERT(mat.num_rows() == 2 && mat.num_cols() == 2, "only for 2x2-matrices");
103 0 : return mat(0,0)*mat(1,1) - mat(1,0)*mat(0,1);
104 : }
105 :
106 : template<typename T>
107 : inline 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 0 : if(invdet == 0.0) return false;
113 0 : invdet = 1.0/invdet;
114 0 : inv(0,0) = mat(1,1) * invdet;
115 0 : inv(1,1) = mat(0,0) * invdet;
116 0 : inv(0,1) = mat(0,1) * -invdet;
117 0 : inv(1,0) = mat(1,0) * -invdet;
118 0 : return true;
119 : }
120 :
121 : template<typename T>
122 0 : bool Invert2(DenseMatrix<T> &mat)
123 : {
124 : double invdet = GetDet2(mat);
125 : //UG_ASSERT(invdet != 0, "Determinant zero, cannot invert matrix.");
126 0 : if(invdet == 0.0) return false;
127 0 : invdet = 1.0/invdet;
128 :
129 : std::swap(mat(0,0), mat(1,1));
130 :
131 0 : mat(0,0) *= invdet;
132 0 : mat(0,1) *= -invdet;
133 0 : mat(1,0) *= -invdet;
134 0 : mat(1,1) *= invdet;
135 0 : return true;
136 : };
137 :
138 :
139 : inline bool GetInverse(DenseMatrix<FixedArray2<double, 2, 2> > &inv, const DenseMatrix<FixedArray2<double, 2, 2> > &mat)
140 : {
141 : return GetInverse2(inv, mat);
142 : }
143 :
144 : inline bool Invert(DenseMatrix< FixedArray2<double, 2, 2> > &mat)
145 : {
146 : return Invert2(mat);
147 : }
148 :
149 : template<typename vector_t, typename matrix_t>
150 0 : inline 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 0 : if(det == 0.0) return false;
157 0 : dest[0] = beta * (mat(1,1)*vec[0] - mat(0,1)*vec[1]) / det;
158 0 : dest[1] = beta * (-mat(1,0)*vec[0] + mat(0,0)*vec[1]) / det;
159 0 : return true;
160 : }
161 :
162 : template<typename T>
163 : inline bool InverseMatMult(DenseVector< FixedArray1<double, 2> > &dest, double beta,
164 : const DenseMatrix< FixedArray2<double, 2, 2> > &mat, const DenseVector< FixedArray1<double, 2> > &vec)
165 : {
166 : return InverseMatMult2(dest, beta, mat, vec);
167 : }
168 :
169 : //////////////////////
170 : // 3x3
171 :
172 :
173 : template<typename T>
174 0 : inline double GetDet3(const DenseMatrix<T> &mat)
175 : {
176 : UG_ASSERT(mat.num_rows() == 3 && mat.num_cols() == 3, "only for 3x3-matrices");
177 0 : 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 0 : - 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 :
181 : template<typename T>
182 0 : inline 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 0 : double invdet = GetDet3(mat);
186 : //UG_ASSERT(invdet != 0, "Determinant zero, cannot invert matrix.");
187 0 : if(invdet == 0.0) return false;
188 0 : invdet = 1.0/invdet;
189 :
190 0 : inv(0,0) = ( mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)) * invdet;
191 0 : inv(0,1) = (-mat(0,1)*mat(2,2) + mat(0,2)*mat(2,1)) * invdet;
192 0 : inv(0,2) = ( mat(0,1)*mat(1,2) - mat(0,2)*mat(1,1)) * invdet;
193 0 : inv(1,0) = (-mat(1,0)*mat(2,2) + mat(1,2)*mat(2,0)) * invdet;
194 0 : inv(1,1) = ( mat(0,0)*mat(2,2) - mat(0,2)*mat(2,0)) * invdet;
195 0 : inv(1,2) = (-mat(0,0)*mat(1,2) + mat(0,2)*mat(1,0)) * invdet;
196 0 : inv(2,0) = ( mat(1,0)*mat(2,1) - mat(1,1)*mat(2,0)) * invdet;
197 0 : inv(2,1) = (-mat(0,0)*mat(2,1) + mat(0,1)*mat(2,0)) * invdet;
198 0 : inv(2,2) = ( mat(0,0)*mat(1,1) - mat(0,1)*mat(1,0)) * invdet;
199 0 : return true;
200 : }
201 :
202 0 : inline bool Invert3(DenseMatrix<FixedArray2<double, 3, 3> > & mat)
203 : {
204 : DenseMatrix<FixedArray2<double, 3, 3> > inv;
205 0 : if(GetInverse3(inv, mat) == false) return false;
206 : mat = inv;
207 : return true;
208 : }
209 :
210 0 : inline bool Invert3(DenseMatrix<VariableArray2<double> > & mat)
211 : {
212 : DenseMatrix<VariableArray2<double> > inv;
213 0 : inv.resize(3,3);
214 0 : if(GetInverse3(inv, mat) == false) return false;
215 0 : mat = inv;
216 : return true;
217 : }
218 :
219 : inline bool GetInverse(DenseMatrix<FixedArray2<double, 3, 3> > &inv, const DenseMatrix<FixedArray2<double, 3, 3> > &mat)
220 : {
221 0 : return GetInverse3(inv, mat);
222 : }
223 :
224 : inline bool Invert(DenseMatrix< FixedArray2<double, 3, 3> > &mat)
225 : {
226 0 : return Invert3(mat);
227 : }
228 :
229 : template<typename vector_t, typename matrix_t>
230 0 : inline bool InverseMatMult3(DenseVector<vector_t> &dest, double beta,
231 : const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
232 : {
233 0 : number det = GetDet3(mat);
234 : //UG_ASSERT(det != 0, "Determinant zero, cannot invert matrix.");
235 : UG_ASSERT(&dest != &vec, "");
236 0 : if(det == 0.0) return false;
237 0 : dest[0] = ( ( mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)) *vec[0] +
238 0 : (-mat(0,1)*mat(2,2) + mat(0,2)*mat(2,1)) *vec[1] +
239 0 : ( mat(0,1)*mat(1,2) - mat(0,2)*mat(1,1)) *vec[2] ) * beta / det;
240 0 : dest[1] = ( (-mat(1,0)*mat(2,2) + mat(1,2)*mat(2,0)) * vec[0] +
241 0 : ( mat(0,0)*mat(2,2) - mat(0,2)*mat(2,0)) * vec[1] +
242 0 : (-mat(0,0)*mat(1,2) + mat(0,2)*mat(1,0)) * vec[2] ) * beta / det;
243 0 : dest[2] = ( ( mat(1,0)*mat(2,1) - mat(1,1)*mat(2,0)) * vec[0] +
244 0 : (-mat(0,0)*mat(2,1) + mat(0,1)*mat(2,0)) * vec[1] +
245 0 : ( mat(0,0)*mat(1,1) - mat(0,1)*mat(1,0)) * vec[2] ) * beta / det;
246 0 : return true;
247 : }
248 :
249 : template<typename T>
250 : inline bool InverseMatMult(DenseVector< FixedArray1<double, 3> > &dest, double beta,
251 : const DenseMatrix< FixedArray2<double, 3, 3> > &mat, const DenseVector< FixedArray1<double, 3> > &vec)
252 : {
253 : return InverseMatMult3(dest, beta, mat, vec);
254 : }
255 : //////////////////////
256 :
257 :
258 :
259 : //! calculates dest = beta1 * A1;
260 : template<typename vector_t, typename matrix_t>
261 0 : inline void MatMult(DenseVector<vector_t> &dest,
262 : const number &beta1, const DenseMatrixInverse<matrix_t> &A1, const DenseVector<vector_t> &w1)
263 : {
264 0 : if(beta1 == 1.0)
265 : {
266 0 : dest = w1;
267 : A1.apply(dest);
268 : }
269 : else
270 : {
271 : DenseVector<vector_t> tmp;
272 0 : tmp = w1;
273 : A1.apply(tmp);
274 0 : VecScaleAssign(dest, beta1, dest);
275 : }
276 0 : }
277 :
278 :
279 : //! calculates dest = alpha1*v1 + beta1 * A1 *w1;
280 : template<typename vector_t, typename matrix_t>
281 : inline void MatMultAdd(DenseVector<vector_t> &dest,
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
286 : DenseVector<vector_t> tmp;
287 : tmp = w1;
288 : A1.apply(tmp);
289 : VecScaleAdd(dest, alpha1, v1, beta1, tmp);
290 : }
291 :
292 :
293 :
294 : template<typename T, eMatrixOrdering TOrdering>
295 : inline bool GetInverse(DenseMatrixInverse<VariableArray2<T, TOrdering> > &inv, const DenseMatrix<VariableArray2<T, TOrdering> > &mat)
296 : {
297 : return inv.set_as_inverse_of(mat);
298 : }
299 :
300 : template<typename T, size_t TBlockSize, eMatrixOrdering TOrdering>
301 : inline bool GetInverse(DenseMatrixInverse<FixedArray2<T, TBlockSize, TBlockSize, TOrdering> > &inv, const DenseMatrix<FixedArray2<T, TBlockSize, TBlockSize, TOrdering> > &mat)
302 : {
303 : return inv.set_as_inverse_of(mat);
304 : }
305 :
306 : //////////////////////
307 : template<typename T>
308 0 : inline bool Invert(DenseMatrix<T> &mat)
309 : {
310 0 : switch(mat.num_rows())
311 : {
312 : case 1: return Invert1(mat);
313 0 : case 2: return Invert2(mat);
314 0 : case 3: return Invert3(mat);
315 0 : default: return InvertNdyn(mat);
316 : }
317 : }
318 :
319 : template<typename vector_t, typename matrix_t>
320 0 : inline 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 0 : if(!GetInverse(inv, mat)) return false;
325 0 : MatMult(dest, beta, inv, vec);
326 : return true;
327 0 : }
328 :
329 :
330 : template<typename vector_t, typename matrix_t>
331 0 : inline bool InverseMatMult(DenseVector<vector_t> &dest, double beta,
332 : const DenseMatrix<matrix_t> &mat, const DenseVector<vector_t> &vec)
333 : {
334 0 : switch(mat.num_rows())
335 : {
336 : case 1: return InverseMatMult1(dest, beta, mat, vec);
337 0 : case 2: return InverseMatMult2(dest, beta, mat, vec);
338 0 : case 3: return InverseMatMult3(dest, beta, mat, vec);
339 0 : default: return InverseMatMultN(dest, beta, mat, vec);
340 : }
341 : }
342 :
343 : ///////////////////////////////////////////
344 :
345 : // end group small_algebra
346 : /// \}
347 :
348 : }
349 :
350 : #endif // __H__UG__CPU_ALGEBRA__DENSEMATRIX_SMALL_INVERSE_H__
|
