ug4
Loading...
Searching...
No Matches
math_matrix_functions_common_impl.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 * Author: Andreas Vogel
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__COMMON__MATH_MATRIX_FUNCTIONS_COMMON_IMPL__
34#define __H__UG__COMMON__MATH_MATRIX_FUNCTIONS_COMMON_IMPL__
35
36#include <cmath>
37#include <iostream>
38#include <iomanip>
39#include <cassert>
40#include "math_matrix.h"
41#include "common/assert.h"
43
44namespace ug{
45
47// Addition of Matrices
49
50template <typename matrix_t>
51inline void
52MatAdd(matrix_t& mOut, const matrix_t& m1, const matrix_t& m2)
53{
54 typedef typename matrix_t::size_type size_type;
55 for(size_type i = 0; i < mOut.num_rows(); ++i)
56 for(size_type j = 0; j < mOut.num_cols(); ++j)
57 {
58 mOut(i,j) = m1(i,j) + m2(i,j);
59 }
60}
61
63// Subtraction of Matrices
65
66template <typename matrix_t>
67inline void
68MatSubtract(matrix_t& mOut, const matrix_t& m1, const matrix_t& m2)
69{
70 typedef typename matrix_t::size_type size_type;
71 for(size_type i = 0; i < mOut.num_rows(); ++i)
72 for(size_type j = 0; j < mOut.num_cols(); ++j)
73 {
74 mOut(i,j) = m1(i,j) - m2(i,j);
75 }
76}
77
78
80// Multiplication of Matrices
82
83template <size_t N, size_t M, size_t L, typename T>
84inline void
86 const MathMatrix<L, M, T>& m2)
87{
88 for(size_t i = 0; i < N; ++i)
89 for(size_t j = 0; j < M; ++j)
90 {
91 mOut(i,j) = 0;
92 for(size_t k = 0; k < L; ++k)
93 {
94 mOut(i,j) += m1(i,k) * m2(k,j);
95 }
96 }
97}
98
99template <size_t N, size_t M, size_t L, size_t P, typename T>
100inline void
102 const MathMatrix<L, P, T>& m2, const MathMatrix<P, M, T>& m3)
103{
105
106 for(size_t i = 0; i < N; ++i)
107 for(size_t j = 0; j < M; ++j)
108 {
109 mOut(i,j) = 0;
110 for(size_t k = 0; k < L; ++k)
111 {
112 help(k,j) = 0;
113 for(size_t l = 0; l < P; ++l)
114 {
115 help(k,j) += m2(k,l) * m3(l,j);
116 }
117
118 mOut(i,j) += m1(i,k) * help(k,j);
119 }
120 }
121}
122
123template <size_t N, size_t M, size_t L, typename T>
124inline void
126 const MathMatrix<M, L, T>& m2)
127{
128 for(size_t i = 0; i < N; ++i)
129 for(size_t j = 0; j < M; ++j)
130 {
131 mOut(i,j) = 0;
132 for(size_t k = 0; k < L; ++k)
133 {
134 mOut(i,j) += m1(k,i) * m2(j,k);
135 }
136 }
137}
138
139template <size_t N, size_t M, typename T>
140inline void
142{
143 for(size_t i = 0; i < N; ++i)
144 {
145 for(size_t j = 0; j < i; ++j)
146 {
147 mOut(i,j) = 0;
148 for(size_t k = 0; k < M; ++k)
149 {
150 mOut(i,j) += m(k,i) * m(k,j);
151 }
152 mOut(j,i) = mOut(i,j);
153 }
154 mOut(i,i) = 0;
155 for(size_t k = 0; k < M; ++k)
156 {
157 mOut(i,i) += m(k,i) * m(k,i);
158 }
159 }
160}
161
162template <size_t N, size_t M, typename T>
163inline void
165{
166 for(size_t i = 0; i < M; ++i)
167 {
168 for(size_t j = 0; j < i; ++j)
169 {
170 mOut(i,j) = 0;
171 for(size_t k = 0; k < N; ++k)
172 {
173 mOut(i,j) += m(i,k) * m(j,k);
174 }
175 mOut(j,i) = mOut(i,j);
176 }
177 mOut(i,i) = 0;
178 for(size_t k = 0; k < N; ++k)
179 {
180 mOut(i,i) += m(i,k) * m(i,k);
181 }
182 }
183}
184
185template <size_t N, size_t M, size_t L, typename T>
186inline void
188 const MathMatrix<M, L, T>& m2)
189{
190 for(size_t i = 0; i < N; ++i)
191 {
192 for(size_t j = 0; j < M; ++j)
193 {
194 mOut(i,j) = 0;
195 for(size_t k = 0; k < L; ++k)
196 {
197 mOut(i,j) += m1(i,k) * m2(j,k);
198 }
199 }
200 }
201}
202
203template <size_t N, size_t M, size_t L, typename T>
204inline void
206 const MathMatrix<L, M, T>& m2)
207{
208 for(size_t i = 0; i < N; ++i)
209 {
210 for(size_t j = 0; j < M; ++j)
211 {
212 mOut(i,j) = 0;
213 for(size_t k = 0; k < L; ++k)
214 {
215 mOut(i,j) += m1(k,i) * m2(k,j);
216 }
217 }
218 }
219}
220
221template <size_t N, size_t M, typename T>
222inline void
224 const MathMatrix<M, M, T>& m2)
225{
227
228 for(size_t i = 0; i < N; ++i)
229 for(size_t j = 0; j < N; ++j)
230 {
231 mOut(i,j) = 0;
232 for(size_t k = 0; k < M; ++k)
233 {
234 help(k,j) = 0;
235 for(size_t l = 0; l < M; ++l)
236 {
237 help(k,j) += m2(k,l) * m1(j,l);
238 }
239
240 mOut(i,j) += m1(i,k) * help(k,j);
241 }
242 }
243}
244
245template <size_t N, size_t M, typename T>
246inline void
248 const MathMatrix<M, M, T>& m2)
249{
251
252 for(size_t i = 0; i < N; ++i)
253 for(size_t j = 0; j < N; ++j)
254 {
255 mOut(i,j) = 0;
256 for(size_t k = 0; k < M; ++k)
257 {
258 help(k,j) = 0;
259 for(size_t l = 0; l < M; ++l)
260 {
261 help(k,j) += m2(k,l) * m1(l,j);
262 }
263
264 mOut(i,j) += m1(k,i) * help(k,j);
265 }
266 }
267}
268
270// "Contraction" for Matrices (note: contraction is usually known regarding tensors!)
272
273template <typename matrix_t>
274inline typename matrix_t::value_type
275MatContraction(const matrix_t& m1, const matrix_t& m2)
276{
277 typename matrix_t::value_type norm = 0;
278 typedef typename matrix_t::size_type size_type;
279 for(size_type i = 0; i < m1.num_rows(); ++i)
280 for(size_type j = 0; j < m1.num_cols(); ++j)
281 {
282 norm += m1(i,j) * m2(i,j);
283 }
284
285 return norm;
286}
287
288
290// "Deviator" and trace for Matrices
292
293template <typename matrix_t>
294inline typename matrix_t::value_type
295MatDeviatorTrace(const matrix_t& m, matrix_t& dev)
296{
297 typename matrix_t::value_type trace = Trace(m);
298
299 typedef typename matrix_t::size_type size_type;
300 for(size_type i = 0; i < m.num_rows(); ++i)
301 {
302 for(size_type j = 0; j < m.num_cols(); ++j)
303 {
304 dev(i,j) = m(i,j);
305 }
306 dev(i,i) -= 1.0 / 3.0 * trace;
307 }
308 return trace;
309}
310
312// Scaling of Matrices
314
315template <typename matrix_t>
316inline void
317MatScale(matrix_t& mOut, typename matrix_t::value_type s, const matrix_t& m)
318{
319 typedef typename matrix_t::size_type size_type;
320 for(size_type i = 0; i < mOut.num_rows(); ++i)
321 for(size_type j = 0; j < mOut.num_cols(); ++j)
322 {
323 mOut(i, j) = m(i, j) * s;
324 }
325}
326
327template <typename matrix_t>
328inline void
329MatScaleAppend(matrix_t& mOut, typename matrix_t::value_type s, const matrix_t& m)
330{
331 typedef typename matrix_t::size_type size_type;
332 for(size_type i = 0; i < mOut.num_rows(); ++i)
333 for(size_type j = 0; j < mOut.num_cols(); ++j)
334 {
335 mOut(i, j) += m(i, j) * s;
336 }
337}
338
340// Transposed of Matrix
342
343template <size_t N, size_t M, typename T>
344inline void
346{
347 typedef typename MathMatrix<N,M,T>::size_type size_type;
348 for(size_type i = 0; i < mOut.num_rows(); ++i)
349 for(size_type j = 0; j < mOut.num_cols(); ++j)
350 {
351 mOut(i, j) = m(j, i);
352 }
353}
354
355template <typename matrix_t>
356inline void
357Transpose(matrix_t& m)
358{
359 UG_ASSERT(m.num_rows()==m.num_cols(), "Transpose: Square Matrix needed");
360
361 typedef typename matrix_t::size_type size_type;
362 matrix_t _temp;
363 for(size_type i = 1; i < m.num_rows(); ++i)
364 for(size_type j = 0; j < i; ++j)
365 _temp(i, j) = m(i, j);
366
367 for(size_type i = 1; i < m.num_rows(); ++i)
368 for(size_type j = 0; j < i; ++j)
369 m(i, j) = m(j, i);
370
371 for(size_type i = 0; i < m.num_rows()-1; ++i)
372 for(size_type j = i+1; j < m.num_cols(); ++j)
373 m(i, j) = _temp(j, i);
374}
375
377// Determinant of Matrix
379
380template <size_t N, typename T>
381inline typename MathMatrix<N,N,T>::value_type
383{
384 UG_THROW("Determinant for matrix of size "<<N<<"x"<<N<<" not implemented.");
385}
386
387template <typename T>
388inline typename MathMatrix<0,0,T>::value_type
390{
391 return 0;
392}
393
394template <typename T>
395inline typename MathMatrix<1,1,T>::value_type
397{
398 return m(0,0);
399}
400
401template <typename T>
402inline typename MathMatrix<2,2,T>::value_type
404{
405 return (m(0,0)*m(1,1) - m(1,0)*m(0,1));
406}
407
408template <typename T>
409inline typename MathMatrix<3,3,T>::value_type
411{
412 return m(0,0)*m(1,1)*m(2,2)
413 + m(0,1)*m(1,2)*m(2,0)
414 + m(0,2)*m(1,0)*m(2,1)
415 - m(0,0)*m(1,2)*m(2,1)
416 - m(0,1)*m(1,0)*m(2,2)
417 - m(0,2)*m(1,1)*m(2,0);
418}
419
420template <typename T>
421inline typename MathMatrix<4,4,T>::value_type
423{
424 return m(0,3)*m(1,2)*m(2,1)*m(3,0)-m(0,2)*m(1,3)*m(2,1)*m(3,0)
425 - m(0,3)*m(1,1)*m(2,2)*m(3,0)+m(0,1)*m(1,3)*m(2,2)*m(3,0)
426 + m(0,2)*m(1,1)*m(2,3)*m(3,0)-m(0,1)*m(1,2)*m(2,3)*m(3,0)
427 - m(0,3)*m(1,2)*m(2,0)*m(3,1)+m(0,2)*m(1,3)*m(2,0)*m(3,1)
428 + m(0,3)*m(1,0)*m(2,2)*m(3,1)-m(0,0)*m(1,3)*m(2,2)*m(3,1)
429 - m(0,2)*m(1,0)*m(2,3)*m(3,1)+m(0,0)*m(1,2)*m(2,3)*m(3,1)
430 + m(0,3)*m(1,1)*m(2,0)*m(3,2)-m(0,1)*m(1,3)*m(2,0)*m(3,2)
431 - m(0,3)*m(1,0)*m(2,1)*m(3,2)+m(0,0)*m(1,3)*m(2,1)*m(3,2)
432 + m(0,1)*m(1,0)*m(2,3)*m(3,2)-m(0,0)*m(1,1)*m(2,3)*m(3,2)
433 - m(0,2)*m(1,1)*m(2,0)*m(3,3)+m(0,1)*m(1,2)*m(2,0)*m(3,3)
434 + m(0,2)*m(1,0)*m(2,1)*m(3,3)-m(0,0)*m(1,2)*m(2,1)*m(3,3)
435 - m(0,1)*m(1,0)*m(2,2)*m(3,3)+m(0,0)*m(1,1)*m(2,2)*m(3,3);
436}
437
439// Gram Determinant of Matrix
441
442template <size_t N, size_t M, typename T>
443inline typename MathMatrix<N,M,T>::value_type
445{
446 if(N <= M)
447 {
449 MatMultiplyMMT(mmT, m);
450 return Determinant(mmT);
451 }
452 else
453 {
455 MatMultiplyMTM(mTm, m);
456 return Determinant(mTm);
457 }
458}
459
460template <typename T>
461inline typename MathMatrix<0,0,T>::value_type
463{
464 return 0;
465}
466
467template <size_t N, typename T>
468inline typename MathMatrix<N,0,T>::value_type
470{
471 return 0;
472}
473
474template <size_t M, typename T>
475inline typename MathMatrix<0,M,T>::value_type
477{
478 return 0;
479}
480
481template <typename T>
482inline typename MathMatrix<1,1,T>::value_type
484{
485 return pow(Determinant(m), 2);
486}
487
488template <typename T>
489inline typename MathMatrix<2,2,T>::value_type
491{
492 return pow(Determinant(m), 2);
493}
494
495template <typename T>
496inline typename MathMatrix<3,3,T>::value_type
498{
499 return pow(Determinant(m), 2);
500}
501
503// Sqrt Gram Determinant of Matrix
505
506template <size_t N, size_t M, typename T>
507inline typename MathMatrix<N,M,T>::value_type
509{
510 return sqrt(GramDeterminant(m));
511}
512
513template <typename T>
514inline typename MathMatrix<0,0,T>::value_type
516{
517 return 0;
518}
519
520template <size_t N, typename T>
521inline typename MathMatrix<N,0,T>::value_type
523{
524 return 0;
525}
526
527template <size_t M, typename T>
528inline typename MathMatrix<0,M,T>::value_type
530{
531 return 0;
532}
533
534template <typename T>
535inline typename MathMatrix<1,1,T>::value_type
537{
538 return fabs(Determinant(m));
539}
540
541template <typename T>
542inline typename MathMatrix<2,2,T>::value_type
544{
545 return fabs(Determinant(m));
546}
547
548template <typename T>
549inline typename MathMatrix<3,3,T>::value_type
551{
552 return fabs(Determinant(m));
553}
554
556// Inverse of Matrix
558template <size_t N, size_t M, typename T>
559inline typename MathMatrix<N,M,T>::value_type
561{
562 UG_THROW("Inverse for matrix of size "<<M<<"x"<<N<<" not implemented. You could use GeneralizedInverse for pseudo-Inverse.");
563}
564
565template <typename T>
566inline typename MathMatrix<1,1,T>::value_type
568{
569 const typename MathMatrix<1,1,T>::value_type det = m(0,0);
570 UG_ASSERT(&mOut != &m, "Inverse: mOut and m have to be different");
571 UG_ASSERT(det != 0, "Inverse: determinate is zero, can not Invert Matrix");
572 mOut(0,0) = 1./m(0,0);
573 return det;
574}
575
576
577template <typename T>
578inline typename MathMatrix<2,2,T>::value_type
580{
581 const typename MathMatrix<2,2,T>::value_type det = Determinant(m);
582 UG_ASSERT(&mOut != &m, "Inverse: mOut and m have to be different");
583 UG_ASSERT(det != 0, "Inverse: determinate is zero, can not Invert Matrix");
584 const typename MathMatrix<2,2,T>::value_type invdet = 1./det;
585
586 mOut(0,0) = m(1,1)*invdet;
587 mOut(1,0) = -m(1,0)*invdet;
588 mOut(0,1) = -m(0,1)*invdet;
589 mOut(1,1) = m(0,0)*invdet;
590
591 return det;
592}
593
594template <typename T>
595inline typename MathMatrix<3,3,T>::value_type
597{
598 const typename MathMatrix<3,3,T>::value_type det = Determinant(m);
599 UG_ASSERT(&mOut != &m, "Inverse: mOut and m have to be different");
600 UG_ASSERT(det != 0, "Inverse: determinate is zero, can not Invert Matrix");
601 const typename MathMatrix<3,3,T>::value_type invdet = 1./det;
602
603 mOut(0,0) = ( m(1,1)*m(2,2) - m(1,2)*m(2,1)) * invdet;
604 mOut(0,1) = (-m(0,1)*m(2,2) + m(0,2)*m(2,1)) * invdet;
605 mOut(0,2) = ( m(0,1)*m(1,2) - m(0,2)*m(1,1)) * invdet;
606 mOut(1,0) = (-m(1,0)*m(2,2) + m(1,2)*m(2,0)) * invdet;
607 mOut(1,1) = ( m(0,0)*m(2,2) - m(0,2)*m(2,0)) * invdet;
608 mOut(1,2) = (-m(0,0)*m(1,2) + m(0,2)*m(1,0)) * invdet;
609 mOut(2,0) = ( m(1,0)*m(2,1) - m(1,1)*m(2,0)) * invdet;
610 mOut(2,1) = (-m(0,0)*m(2,1) + m(0,1)*m(2,0)) * invdet;
611 mOut(2,2) = ( m(0,0)*m(1,1) - m(0,1)*m(1,0)) * invdet;
612
613 return det;
614}
615
617// Inverse Transposed of Matrix
619
620template <size_t N, size_t M, typename T>
621inline typename MathMatrix<N,M,T>::value_type
623{
624 UG_THROW("InverseTransposed for matrix of size "<<M<<"x"<<N<<" not implemented.");
625}
626
627template <typename T>
628inline typename MathMatrix<1,1,T>::value_type
630{
631 return Inverse(mOut, m);
632}
633
634template <typename T>
635inline typename MathMatrix<2,2,T>::value_type
637{
638 const typename MathMatrix<2,2,T>::value_type det = Determinant(m);
639 UG_ASSERT(&mOut != &m, "InverseTransposed: mOut and m have to be different");
640 UG_ASSERT(det != 0, "InverseTransposed: determinate is zero, can not Invert Matrix");
641 const typename MathMatrix<2,2,T>::value_type invdet = 1./det;
642
643 mOut(0,0) = m(1,1)*invdet;
644 mOut(1,0) = -m(0,1)*invdet;
645 mOut(0,1) = -m(1,0)*invdet;
646 mOut(1,1) = m(0,0)*invdet;
647
648 return det;
649}
650
651template <typename T>
652inline typename MathMatrix<3,3,T>::value_type
654{
655 const typename MathMatrix<3,3,T>::value_type det = Determinant(m);
656 UG_ASSERT(&mOut != &m, "InverseTransposed: mOut and m have to be different");
657 UG_ASSERT(det != 0, "InverseTransposed: determinate is zero, can not Invert Matrix");
658 const typename MathMatrix<3,3,T>::value_type invdet = 1./det;
659
660 mOut(0,0) = ( m(1,1)*m(2,2) - m(2,1)*m(1,2)) * invdet;
661 mOut(0,1) = (-m(1,0)*m(2,2) + m(2,0)*m(1,2)) * invdet;
662 mOut(0,2) = ( m(1,0)*m(2,1) - m(2,0)*m(1,1)) * invdet;
663 mOut(1,0) = (-m(0,1)*m(2,2) + m(2,1)*m(0,2)) * invdet;
664 mOut(1,1) = ( m(0,0)*m(2,2) - m(2,0)*m(0,2)) * invdet;
665 mOut(1,2) = (-m(0,0)*m(2,1) + m(2,0)*m(0,1)) * invdet;
666 mOut(2,0) = ( m(0,1)*m(1,2) - m(1,1)*m(0,2)) * invdet;
667 mOut(2,1) = (-m(0,0)*m(1,2) + m(1,0)*m(0,2)) * invdet;
668 mOut(2,2) = ( m(0,0)*m(1,1) - m(1,0)*m(0,1)) * invdet;
669
670 return det;
671}
672
674// Right-Inverse of Matrix
676
677template <size_t N, size_t M, typename T>
678inline typename MathMatrix<N,M,T>::value_type
680{
681 //UG_STATIC_ASSERT(M <= N, pseudo_inverse_does_not_exist);
682 if (M > N) // note: this is a 'static condition', it should be eliminated by the optimizer
683 UG_THROW ("RightInverse: Type mismatch, cannot right-invert a MxN-matrix with M > N!");
684
685 // H = m*mT (H is symmetric)
686 // TODO: Since H is symmetric, we could store only lower or upper elements
687 MathMatrix<M,M,T> H, HInv;
688 MatMultiplyMMT(H, m);
689 // Invert H
690 const number det = Inverse(HInv, H);
691
692 MatMultiplyTransposed(mOut, m, HInv);
693
694 return sqrt(det);
695}
696
697template <typename T>
698inline typename MathMatrix<1,1,T>::value_type
699RightInverse(MathMatrix<1,1>& mOut, const MathMatrix<1,1>& m){return fabs(Inverse(mOut, m));}
700
701template <typename T>
702inline typename MathMatrix<2,2,T>::value_type
703RightInverse(MathMatrix<2,2>& mOut, const MathMatrix<2,2>& m){return fabs(Inverse(mOut, m));}
704
705template <typename T>
706inline typename MathMatrix<3,3,T>::value_type
707RightInverse(MathMatrix<3,3>& mOut, const MathMatrix<3,3>& m){return fabs(Inverse(mOut, m));}
708
710// Left-Inverse of Matrix
712
713template <size_t N, size_t M, typename T>
714inline typename MathMatrix<N,M,T>::value_type
716{
717 //UG_STATIC_ASSERT(N <= M, pseudo_inverse_does_not_exist);
718 if (N > M) // note: this is a 'static condition', it should be eliminated by the optimizer
719 UG_THROW ("LeftInverse: Type mismatch, cannot right-invert a MxN-matrix with M < N!");
720
721 // H = mT*m (H is symmetric)
722 // TODO: Since H is symmetric, we could store only lower or upper elements
723 MathMatrix<N,N,T> H, HInv;
724 MatMultiplyMTM(H, m);
725
726 // Invert H
727 const number det = Inverse(HInv, H);
728
729 MatMultiplyTransposed(mOut, HInv, m);
730
731 return sqrt(det);
732}
733
734template <typename T>
735inline typename MathMatrix<1,1,T>::value_type
736LeftInverse(MathMatrix<1,1>& mOut, const MathMatrix<1,1>& m){return fabs(Inverse(mOut, m));}
737
738template <typename T>
739inline typename MathMatrix<2,2,T>::value_type
740LeftInverse(MathMatrix<2,2>& mOut, const MathMatrix<2,2>& m){return fabs(Inverse(mOut, m));}
741
742template <typename T>
743inline typename MathMatrix<3,3,T>::value_type
744LeftInverse(MathMatrix<3,3>& mOut, const MathMatrix<3,3>& m){return fabs(Inverse(mOut, m));}
745
747// Generalized-Inverse of Matrix
749template<size_t N, size_t M, typename T>
750inline typename MathMatrix<N,M,T>::value_type
752{
753 if(M<N){//UG_LOG("Right Inverse for matrix of size "<<M<<"x"<<N<<".");
754 return RightInverse(mOut,m);
755 }
756
757 if(M>N){//UG_LOG("Left Inverse for matrix of size "<<M<<"x"<<N<<".");
758 return LeftInverse(mOut,m);
759 }
760 return Inverse(mOut,m);
761}
762
764// Trace of Matrix
766
767template <typename T>
768inline typename MathMatrix<1,1,T>::value_type
770{
771 return m(0,0);
772}
773
774template <typename T>
775inline typename MathMatrix<2,2,T>::value_type
777{
778 return (m(0,0)+m(1,1));
779}
780
781template <typename T>
782inline typename MathMatrix<3,3,T>::value_type
784{
785 return (m(0,0)+m(1,1)+m(2,2));
786}
787
789// Scalar operations for Matrices
791
792template <typename matrix_t>
793inline void
794MatSet(matrix_t& mInOut, typename matrix_t::value_type s)
795{
796 typedef typename matrix_t::size_type size_type;
797 for(size_type i = 0; i < mInOut.num_rows(); ++i)
798 for(size_type j = 0; j < mInOut.num_cols(); ++j)
799 {
800 mInOut(i, j) = s;
801 }
802}
803
804template <typename matrix_t>
805inline void
806MatDiagSet(matrix_t& mInOut, typename matrix_t::value_type s)
807{
808 typedef typename matrix_t::size_type size_type;
809 for(size_type i = 0; i < mInOut.num_rows(); ++i)
810 {
811 mInOut(i, i) = s;
812 }
813}
814
815template <typename matrix_t>
816inline void
817MatAdd(matrix_t& mOut, const matrix_t& m, typename matrix_t::value_type s)
818{
819 typedef typename matrix_t::size_type size_type;
820 for(size_type i = 0; i < mOut.num_rows(); ++i)
821 for(size_type j = 0; j < mOut.num_cols(); ++j)
822 {
823 mOut(i, j) = m(i,j) + s;
824 }
825}
826
827template <typename matrix_t>
828inline void
829MatSubtract(matrix_t& mOut, const matrix_t& m, typename matrix_t::value_type s)
830{
831 typedef typename matrix_t::size_type size_type;
832 for(size_type i = 0; i < mOut.num_rows(); ++i)
833 for(size_type j = 0; j < mOut.num_cols(); ++j)
834 {
835 mOut(i, j) = m(i,j) - s;
836 }
837}
838
839template <typename matrix_t>
840inline void
841MatDivide(matrix_t& mOut, const matrix_t& m, typename matrix_t::value_type s)
842{
843 typedef typename matrix_t::size_type size_type;
844 for(size_type i = 0; i < mOut.num_rows(); ++i)
845 for(size_type j = 0; j < mOut.num_cols(); ++j)
846 {
847 mOut(i, j) = m(i,j) /s;
848 }
849}
850
851template <typename matrix_t>
852inline void
853MatMultiply(matrix_t& mOut, const matrix_t& m, typename matrix_t::value_type s)
854{
855 typedef typename matrix_t::size_type size_type;
856 for(size_type i = 0; i < mOut.num_rows(); ++i)
857 for(size_type j = 0; j < mOut.num_cols(); ++j)
858 {
859 mOut(i, j) = m(i,j) * s;
860 }
861}
862
863template <typename matrix_t>
864inline void
865MatIdentity(matrix_t& mOut)
866{
867 MatSet(mOut, 0);
868 MatDiagSet(mOut, 1);
869}
870
871template <typename matrix_t>
872inline void
873MatRotationX(matrix_t& mOut, typename matrix_t::value_type rads)
874{
875 //UG_STATIC_ASSERT(matrix_t::RowSize > 2, AT_LEAST_3_ROWS_REQUIRED);
876 //UG_STATIC_ASSERT(matrix_t::ColSize > 2, AT_LEAST_3_COLS_REQUIRED);
877
878 MatIdentity(mOut);
879 typename matrix_t::value_type s = sin(rads);
880 typename matrix_t::value_type c = cos(rads);
881
882 mOut(1, 1) = c; mOut(1, 2) = -s;
883 mOut(2, 1) = s; mOut(2, 2) = c;
884}
885
886template <typename matrix_t>
887inline void
888MatRotationY(matrix_t& mOut, typename matrix_t::value_type rads)
889{
890 //UG_STATIC_ASSERT(matrix_t::RowSize > 2, AT_LEAST_3_ROWS_REQUIRED);
891 //UG_STATIC_ASSERT(matrix_t::ColSize > 2, AT_LEAST_3_COLS_REQUIRED);
892
893 MatIdentity(mOut);
894 typename matrix_t::value_type s = sin(rads);
895 typename matrix_t::value_type c = cos(rads);
896
897 mOut(0, 0) = c; mOut(0, 2) = s;
898 mOut(2, 0) = -s; mOut(2, 2) = c;
899}
900
901template <typename matrix_t>
902inline void
903MatRotationZ(matrix_t& mOut, typename matrix_t::value_type rads)
904{
905 MatIdentity(mOut);
906 typename matrix_t::value_type s = sin(rads);
907 typename matrix_t::value_type c = cos(rads);
908
909 mOut(0, 0) = c; mOut(0, 1) = -s;
910 mOut(1, 0) = s; mOut(1, 1) = c;
911}
912
916
917template <typename matrix_t>
918inline void
920 typename matrix_t::value_type yaw,
921 typename matrix_t::value_type pitch,
922 typename matrix_t::value_type roll)
923{
924 //UG_STATIC_ASSERT(matrix_t::RowSize > 2, AT_LEAST_3_ROWS_REQUIRED);
925 //UG_STATIC_ASSERT(matrix_t::ColSize > 2, AT_LEAST_3_COLS_REQUIRED);
926
927 matrix_t tMat1, tMat2, tMat3;
928 MatRotationX(tMat1, yaw);
929 MatRotationY(tMat2, pitch);
930 MatMultiply(tMat3, tMat1, tMat2);
931 MatRotationZ(tMat1, roll);
932 MatMultiply(mOut, tMat3, tMat1);
933}
934
939
940template <typename matrix_t, typename vector_t>
941inline void
942MatHouseholder(matrix_t& mOut, const vector_t& orthoVec)
943{
944 assert(vector_t::Size == matrix_t::RowSize);
945 assert(vector_t::Size == matrix_t::ColSize);
946
947 typename vector_t::value_type scalarProd = VecDot(orthoVec, orthoVec);
948
949 typedef typename matrix_t::size_type size_type_mat;
950 for(size_type_mat i = 0; i < mOut.num_rows(); ++i)
951 {
952 for(size_type_mat j = 0; j < mOut.num_cols(); ++j){
953 mOut(i,j) = - 2.0/scalarProd * orthoVec[i] * orthoVec[j];
954 }
955 mOut(i,i) += 1.0;
956 }
957
958}
959
961// Norms for Matrices
963
964template <typename matrix_t>
965inline typename matrix_t::value_type
967{
968 typename matrix_t::value_type norm = 0;
969 typedef typename matrix_t::size_type size_type;
970 for(size_type i = 0; i < m.num_rows(); ++i)
971 for(size_type j = 0; j < m.num_cols(); ++j)
972 {
973 norm += m(i,j)*m(i,j);
974 }
975
976 return norm;
977}
978
979template <typename matrix_t>
980inline typename matrix_t::value_type
982{
983 return static_cast<typename matrix_t::value_type>(sqrt(MatFrobeniusNormSq(m)));
984}
985
986template <typename matrix_t>
987inline typename matrix_t::value_type
988MatOneNorm(matrix_t& m)
989{
990 typename matrix_t::value_type sum, max = 0;
991 typedef typename matrix_t::size_type size_type;
992 for(size_type j = 0; j < m.num_cols(); ++j)
993 {
994 sum = 0;
995 for(size_type i = 0; i < m.num_rows(); ++i)
996 {
997 sum += fabs(m(i,j));
998 }
999 max = (sum > max) ? sum : max;
1000 }
1001 return max;
1002}
1003
1004template <typename matrix_t>
1005inline typename matrix_t::value_type
1006MatInftyNorm(matrix_t& m)
1007{
1008 typename matrix_t::value_type sum, max = 0;
1009 typedef typename matrix_t::size_type size_type;
1010 for(size_type i = 0; i < m.num_rows(); ++i)
1011 {
1012 sum = 0;
1013 for(size_type j = 0; j < m.num_cols(); ++j)
1014 {
1015 sum += fabs(m(i,j));
1016 }
1017 max = (sum > max) ? sum : max;
1018 }
1019 return max;
1020}
1021
1022template <typename matrix_t>
1023inline typename matrix_t::value_type
1024MatMaxNorm(matrix_t& m)
1025{
1026 typename matrix_t::value_type max = 0;
1027 typedef typename matrix_t::size_type size_type;
1028 for(size_type i = 0; i < m.num_rows(); ++i)
1029 for(size_type j = 0; j < m.num_cols(); ++j)
1030 {
1031 max = (m(i,j) > max) ? m(i,j) : max;
1032 }
1033
1034 return max;
1035}
1036
1037
1038
1039template <size_t N, size_t M, typename T>
1040inline typename MathMatrix<N,M,T>::value_type
1042{
1043 UG_THROW("MaxAbsEigenvalue for matrix of size "<<N<<"x"<<M<<" not implemented.");
1044}
1045
1046template <typename T>
1047inline typename MathMatrix<1,1,T>::value_type
1049{
1050 const typename MathMatrix<1,1,T>::value_type val = m(0,0);
1051 return fabs(val);
1052}
1053
1054template <typename T>
1055inline typename MathMatrix<2,2,T>::value_type
1057{
1058 typename MathMatrix<2,2,T>::value_type minus_p_half, val;
1059 minus_p_half = m(0,0)+m(1,1);
1060 val = minus_p_half*minus_p_half - (m(0,0)*m(1,1) - m(0,1)*m(1,0));
1061 UG_ASSERT(val >= 0.0, "MaxAbsEigenvalues: Complex Eigenvalues???");
1062
1063 if (minus_p_half >=0.0) {return (minus_p_half + sqrt(val));}
1064 else {return fabs(minus_p_half-sqrt(val));}
1065}
1066
1067
1068template <typename matrix_t>
1069inline typename matrix_t::value_type
1071{
1072 matrix_t inv;
1073 Inverse(inv, m);
1074 typename matrix_t::value_type min=1.0/MaxAbsEigenvalue(inv);
1075 return min;
1076}
1077
1078} // end of namespace
1079
1080#endif /* __H__UG__COMMON__MATH_MATRIX_FUNCTIONS_COMMON_IMPL__ */
parameterString s
A class for fixed size, dense matrices.
Definition math_matrix.h:63
std::size_t size_type
Definition math_matrix.h:66
T value_type
Definition math_matrix.h:65
std::size_t num_cols() const
Definition math_matrix.h:216
std::size_t num_rows() const
Definition math_matrix.h:215
void MatMultiplyMMT(MathMatrix< M, M, T > &mOut, const MathMatrix< M, N, T > &m)
multiply a matrix with its transposed and stores the result in a second one
Definition math_matrix_functions_common_impl.hpp:164
MathMatrix< N, M, T >::value_type MinAbsEigenvalue(const MathMatrix< M, N, T > &m)
Computes minimum eigenvalue of a (symmetric) matrix.
matrix_t::value_type MatDeviatorTrace(const matrix_t &m, matrix_t &dev)
Definition math_matrix_functions_common_impl.hpp:295
MathMatrix< N, M, T >::value_type SqrtGramDeterminant(const MathMatrix< N, M, T > &m)
Square root of Gram Determinant of a matrix.
Definition math_matrix_functions_common_impl.hpp:508
matrix_t::value_type MatContraction(const matrix_t &m1, const matrix_t &m2)
Definition math_matrix_functions_common_impl.hpp:275
MathMatrix< N, M, T >::value_type RightInverse(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
Right-Inverse of a Matrix.
Definition math_matrix_functions_common_impl.hpp:679
void MatMultiplyMBMT(MathMatrix< N, N, T > &mOut, const MathMatrix< N, M, T > &m1, const MathMatrix< M, M, T > &m2)
Definition math_matrix_functions_common_impl.hpp:223
void MatRotationX(matrix_t &mOut, typename matrix_t::value_type rads)
Fills the matrix with a matrix that rotates around the x-axis in 3 dimensions.
Definition math_matrix_functions_common_impl.hpp:873
void Transpose(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
transpose a matrix
Definition math_matrix_functions_common_impl.hpp:345
MathMatrix< 1, 1, T >::value_type Trace(const MathMatrix< 1, 1, T > &m)
Trace of a Matrix.
Definition math_matrix_functions_common_impl.hpp:769
matrix_t::value_type MatInftyNorm(matrix_t &m)
Definition math_matrix_functions_common_impl.hpp:1006
void MatMultiplyMTB(MathMatrix< N, M, T > &mOut, const MathMatrix< L, N, T > &m1, const MathMatrix< L, M, T > &m2)
multiply the transposed of a matrix with a matrix and stores the result in mOut
Definition math_matrix_functions_common_impl.hpp:205
MathMatrix< N, M, T >::value_type GeneralizedInverse(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
Definition math_matrix_functions_common_impl.hpp:751
void MatMultiplyMBT(MathMatrix< N, M, T > &mOut, const MathMatrix< N, L, T > &m1, const MathMatrix< M, L, T > &m2)
multiply a matrix with the transposed of a second one and stores the result in mOut
Definition math_matrix_functions_common_impl.hpp:187
void MatMultiplyTransposed(MathMatrix< N, M, T > &mOut, const MathMatrix< L, N, T > &m1, const MathMatrix< M, L, T > &m2)
multiply two transposed matrices and stores the result in a third one
Definition math_matrix_functions_common_impl.hpp:125
void MatIdentity(matrix_t &mOut)
Fills the matrix with the identity matrix.
Definition math_matrix_functions_common_impl.hpp:865
void MatSet(matrix_t &mInOut, typename matrix_t::value_type s)
Set each matrix entry to a scalar (componentwise)
Definition math_matrix_functions_common_impl.hpp:794
MathMatrix< N, N, T >::value_type Determinant(const MathMatrix< N, N, T > &m)
Determinant of a matrix.
Definition math_matrix_functions_common_impl.hpp:382
void MatDiagSet(matrix_t &mInOut, typename matrix_t::value_type s)
Set diagonal entries of a matrix to a scalar (other entries are not changed)
Definition math_matrix_functions_common_impl.hpp:806
void MatRotationZ(matrix_t &mOut, typename matrix_t::value_type rads)
Fills the matrix with a matrix that rotates around the y-axis in 2 or 3 dimensions.
Definition math_matrix_functions_common_impl.hpp:903
matrix_t::value_type MatOneNorm(matrix_t &m)
Definition math_matrix_functions_common_impl.hpp:988
void MatRotationY(matrix_t &mOut, typename matrix_t::value_type rads)
Fills the matrix with a matrix that rotates around the y-axis in 3 dimensions.
Definition math_matrix_functions_common_impl.hpp:888
matrix_t::value_type MatFrobeniusNormSq(matrix_t &m)
Definition math_matrix_functions_common_impl.hpp:966
matrix_t::value_type MatMaxNorm(matrix_t &m)
Definition math_matrix_functions_common_impl.hpp:1024
void MatMultiplyMTM(MathMatrix< N, N, T > &mOut, const MathMatrix< M, N, T > &m)
multiply a transposed matrix with itself and stores the result in a second one
Definition math_matrix_functions_common_impl.hpp:141
MathMatrix< N, M, T >::value_type GramDeterminant(const MathMatrix< N, M, T > &m)
Gram Determinant of a matrix.
Definition math_matrix_functions_common_impl.hpp:444
void MatAdd(matrix_t &mOut, const matrix_t &m1, const matrix_t &m2)
adds two matrices and stores the result in a third one
Definition math_matrix_functions_common_impl.hpp:52
void MatScale(matrix_t &mOut, typename matrix_t::value_type s, const matrix_t &m)
scales a matrix_t
Definition math_matrix_functions_common_impl.hpp:317
void MatRotationYawPitchRoll(matrix_t &mOut, typename matrix_t::value_type yaw, typename matrix_t::value_type pitch, typename matrix_t::value_type roll)
Creates a rotation matrix given yaw, pitch and roll in radiants.
Definition math_matrix_functions_common_impl.hpp:919
MathMatrix< N, M, T >::value_type LeftInverse(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
Left-Inverse of a Matrix.
Definition math_matrix_functions_common_impl.hpp:715
matrix_t::value_type MatFrobeniusNorm(matrix_t &m)
Definition math_matrix_functions_common_impl.hpp:981
MathMatrix< N, M, T >::value_type InverseTransposed(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
Transposed-Inverse of a Matrix (= Inverse-Transposed of a Matrix)
Definition math_matrix_functions_common_impl.hpp:622
void MatMultiply(MathMatrix< N, M, T > &mOut, const MathMatrix< N, L, T > &m1, const MathMatrix< L, M, T > &m2)
multiply two matrices and stores the result in a third one
Definition math_matrix_functions_common_impl.hpp:85
void MatMultiplyMTBM(MathMatrix< N, N, T > &mOut, const MathMatrix< M, N, T > &m1, const MathMatrix< M, M, T > &m2)
Definition math_matrix_functions_common_impl.hpp:247
MathMatrix< N, M, T >::value_type Inverse(MathMatrix< N, M, T > &mOut, const MathMatrix< M, N, T > &m)
Inverse of a matrix.
Definition math_matrix_functions_common_impl.hpp:560
void MatSubtract(matrix_t &mOut, const matrix_t &m1, const matrix_t &m2)
subtracts m2 from m1 and stores the result in a mOut
Definition math_matrix_functions_common_impl.hpp:68
MathMatrix< N, M, T >::value_type MaxAbsEigenvalue(const MathMatrix< M, N, T > &m)
Computes maximum eigenvalue of a (symmetric) matrix.
Definition math_matrix_functions_common_impl.hpp:1041
void MatScaleAppend(matrix_t &mOut, typename matrix_t::value_type s, const matrix_t &m)
scales a matrix_t and adds to result to a second matrix
Definition math_matrix_functions_common_impl.hpp:329
void MatHouseholder(matrix_t &mOut, const vector_t &orthoVec)
Definition math_matrix_functions_common_impl.hpp:942
void MatDivide(matrix_t &mOut, const matrix_t &m, typename matrix_t::value_type s)
Devide a matrix by a scalar (componentwise)
Definition math_matrix_functions_common_impl.hpp:841
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_THROW(msg)
Definition error.h:57
double number
Definition types.h:124
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
returns the dot-product of two vector_ts
Definition math_vector_functions_common_impl.hpp:385
the ug namespace