1 #ifndef AMGCL_VALUE_TYPE_EIGEN_HPP
2 #define AMGCL_VALUE_TYPE_EIGEN_HPP
3
4 /*
5 The MIT License
6
7 Copyright (c) 2012-2021 Denis Demidov <dennis.demidov@gmail.com>
8
9 Permission is hereby granted, free of charge, to any person obtaining a copy
10 of this software and associated documentation files (the "Software"), to deal
11 in the Software without restriction, including without limitation the rights
12 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 copies of the Software, and to permit persons to whom the Software is
14 furnished to do so, subject to the following conditions:
15
16 The above copyright notice and this permission notice shall be included in
17 all copies or substantial portions of the Software.
18
19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
25 THE SOFTWARE.
26 */
27
28 /**
29 * \file amgcl/value_type/eigen.hpp
30 * \author Denis Demidov <dennis.demidov@gmail.com>
31 * \brief Enable statically sized eigen matrices as value types.
32 */
33
34 #include <Eigen/Dense>
35
36 #include <amgcl/backend/builtin.hpp>
37 #include <amgcl/value_type/interface.hpp>
38
39 namespace amgcl {
40 namespace backend {
41
42 /// Enable Eigen matrix as a value-type.
43 template <typename T, int N, int M>
44 struct is_builtin_vector< std::vector<Eigen::Matrix<T, N, M> > > : std::true_type {};
45
46 } // namespace backend
47
48 namespace math {
49
50 /// Scalar type of a non-scalar type.
51 template <class T, int N, int M>
52 struct scalar_of< Eigen::Matrix<T, N, M> > {
53 typedef typename math::scalar_of<T>::type type;
54 };
55
56 /// Replace scalar type in the static matrix
57 template <class T, int N, int M, class S>
58 struct replace_scalar< Eigen::Matrix<T, N, M>, S> {
59 typedef Eigen::Matrix<S, N, M> type;
60 };
61
62 /// RHS type corresponding to a non-scalar type.
63 template <class T, int N>
64 struct rhs_of< Eigen::Matrix<T, N, N> > {
65 typedef Eigen::Matrix<T, N, 1> type;
66 };
67
68 /// Whether the value type is a statically sized matrix.
69 template <class T, int N, int M>
70 struct is_static_matrix< Eigen::Matrix<T, N, M> > : std::true_type {};
71
72 /// Number of rows for statically sized matrix types.
73 template <class T, int N, int M>
74 struct static_rows< Eigen::Matrix<T, N, M> > : std::integral_constant<int, N> {};
75
76 /// Number of columns for statically sized matrix types.
77 template <class T, int N, int M>
78 struct static_cols< Eigen::Matrix<T, N, M> > : std::integral_constant<int, M> {};
79
80 /// Specialization of conjugate transpose for eigen matrices.
81 template <typename T, int N, int M>
82 struct adjoint_impl< Eigen::Matrix<T, N, M> >
83 {
84 typedef typename Eigen::Matrix<T, N, M>::AdjointReturnType return_type;
85
getamgcl::math::adjoint_impl86 static return_type get(const Eigen::Matrix<T, N, M> &x) {
87 return x.adjoint();
88 }
89 };
90
91 /// Inner-product result of two Eigen vectors.
92 template <class T, int N>
93 struct inner_product_impl< Eigen::Matrix<T, N, 1> >
94 {
95 typedef T return_type;
getamgcl::math::inner_product_impl96 static T get(const Eigen::Matrix<T, N, 1> &x, const Eigen::Matrix<T, N, 1> &y) {
97 return x.adjoint() * y;
98 }
99 };
100
101 /// Inner-product result of two Eigen matrices.
102 template <class T, int N, int M>
103 struct inner_product_impl< Eigen::Matrix<T, N, M> >
104 {
105 typedef Eigen::Matrix<T, M, M> return_type;
106
getamgcl::math::inner_product_impl107 static return_type get(const Eigen::Matrix<T, N, M> &x, const Eigen::Matrix<T, N, M> &y) {
108 return x.adjoint() * y;
109 }
110 };
111
112 /// Specialization of element norm for eigen matrices.
113 template <typename T, int N, int M>
114 struct norm_impl< Eigen::Matrix<T, N, M> >
115 {
getamgcl::math::norm_impl116 static typename math::scalar_of<T>::type get(const Eigen::Matrix<T, N, M> &x) {
117 return x.norm();
118 }
119 };
120
121 /// Specialization of zero element for eigen matrices.
122 template <typename T, int N, int M>
123 struct zero_impl< Eigen::Matrix<T, N, M> >
124 {
getamgcl::math::zero_impl125 static Eigen::Matrix<T, N, M> get() {
126 return Eigen::Matrix<T, N, M>::Zero();
127 }
128 };
129
130 /// Specialization of zero element for eigen matrices.
131 template <typename T, int N, int M>
132 struct is_zero_impl< Eigen::Matrix<T, N, M> >
133 {
getamgcl::math::is_zero_impl134 static bool get(const Eigen::Matrix<T, N, M> &x) {
135 return x.isZero();
136 }
137 };
138
139 /// Specialization of identity for eigen matrices.
140 template <typename T, int N>
141 struct identity_impl< Eigen::Matrix<T, N, N> >
142 {
getamgcl::math::identity_impl143 static Eigen::Matrix<T, N, N> get() {
144 return Eigen::Matrix<T, N, N>::Identity();
145 }
146 };
147
148 /// Specialization of constant for eigen matrices.
149 template <typename T, int N, int M>
150 struct constant_impl< Eigen::Matrix<T, N, M> >
151 {
getamgcl::math::constant_impl152 static Eigen::Matrix<T, N, M> get(T c) {
153 return Eigen::Matrix<T, N, M>::Constant(c);
154 }
155 };
156
157 /// Specialization of inversion for eigen matrices.
158 template <typename T, int N>
159 struct inverse_impl< Eigen::Matrix<T, N, N> >
160 {
getamgcl::math::inverse_impl161 static Eigen::Matrix<T, N, N> get(const Eigen::Matrix<T, N, N> &x) {
162 return x.inverse();
163 }
164 };
165
166 } // namespace math
167 } // namespace amgcl
168
169 namespace Eigen {
170
171 template <class A, class B>
operator <(const MatrixBase<A> & a,const MatrixBase<B> & b)172 bool operator<(const MatrixBase<A> &a, const MatrixBase<B> &b) {
173 return a.trace() < b.trace();
174 }
175
176 } // namespace Eigen
177
178 #endif
179