1 /** @file gsMatrixAddons.h
2
3 @brief Provides extra member declarations to the Eigen MatrixBase class
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s): A. Mantzaflaris
12 */
13
14 inline const internal::adjugate_impl<Derived> adjugate() const;
15
16 inline void adjugateInPlace();
17
18 typedef BlockDiag<Derived,Dynamic> BlockDiagReturnType;
19 inline const BlockDiagReturnType blockDiag(Index rowFactor) const;
20
21 typedef BlockTranspose<Derived,Dynamic> BlockTransposeReturnType;
22 inline const BlockTransposeReturnType blockTranspose(Index rowFactor) const;
23
24 template<typename IndicesType>
25 const RowSelection<Derived,IndicesType> selectRows(const IndicesType & ind) const;
26
27
28 /**
29 * \brief Simple (inplace) Gauss elimination without any pivoting
30 */
gaussElim()31 void gaussElim()
32 {
33 Derived & M = derived();
34 index_t piv = 0;
35 const index_t nr = M.rows();
36 const index_t nc = M.cols();
37 for (index_t r=0; r!=nr; ++r)
38 {
39 piv = 0;
40 while (piv!= nc && 0 == M(r, piv)) ++piv;
41 if (piv == nc ) continue;
42
43 const index_t br = nr-r-1;
44 const index_t bc = nc-piv-1;
45 M.block(r+1, piv+1, br, bc).noalias() -=
46 M.col(piv).tail(br) * M.row(r).tail(bc) / M(r, piv);
47 M.col(piv).tail(br).setZero();
48 }
49 }
50
51
52 /**
53 * \brief Inversion for small matrices using Cramer's Rule
54 */
cramerInverse()55 inline Matrix<Scalar, Dynamic, Dynamic> cramerInverse() const
56 {
57 const Derived & M = derived();
58 eigen_assert(M.rows() == M.cols() && "Matrix is not square.");
59
60 Matrix<Scalar, Dynamic, Dynamic> rvo1(M.rows(), M.rows());
61 switch (M.rows())
62 {
63 case 1:
64 rvo1 = M.template topLeftCorner<1, 1>().inverse();
65 break;
66 case 2:
67 rvo1 = M.template topLeftCorner<2, 2>().inverse();
68 break;
69 case 3:
70 rvo1 = M.template topLeftCorner<3, 3>().inverse();
71 break;
72 case 4:
73 rvo1 = M.template topLeftCorner<4, 4>().inverse();
74 break;
75 default:
76 eigen_assert(false && "Not implemented.");
77 };
78 return rvo1;
79 }
80
81 /**
82 * \brief Inplace inversion for small matrices using Cramer's Rule
83 */
cramerInverseInPlace()84 inline void cramerInverseInPlace()
85 {
86 // derived() = cramerInverse().eval();
87 derived().swap(cramerInverse().eval());
88 }