Functions
product.h File Reference
#include "main.h"
#include <Eigen/QR>

Go to the source code of this file.

Functions

template<typename Derived1 , typename Derived2 >
bool areNotApprox (const MatrixBase< Derived1 > &m1, const MatrixBase< Derived2 > &m2, typename Derived1::RealScalar epsilon=NumTraits< typename Derived1::RealScalar >::dummy_precision())
 
template<typename Type1 , typename Type2 , typename Tol >
bool verifyIsApprox (const Type1 &a, const Type2 &b, Tol tol)
 
template<typename LhsType , typename RhsType >
std::enable_if_t< RhsType::SizeAtCompileTime==Dynamic, void > check_mismatched_product (LhsType &lhs, const RhsType &rhs)
 
template<typename LhsType , typename RhsType >
std::enable_if_t< RhsType::SizeAtCompileTime !=Dynamic, void > check_mismatched_product (LhsType &, const RhsType &)
 
template<typename MatrixType >
void product (const MatrixType &m)
 

Function Documentation

◆ areNotApprox()

template<typename Derived1 , typename Derived2 >
bool areNotApprox ( const MatrixBase< Derived1 > &  m1,
const MatrixBase< Derived2 > &  m2,
typename Derived1::RealScalar  epsilon = NumTraits<typename Derived1::RealScalar>::dummy_precision() 
)
15  {
16  return !((m1 - m2).cwiseAbs2().maxCoeff() <
17  epsilon * epsilon * (std::max)(m1.cwiseAbs2().maxCoeff(), m2.cwiseAbs2().maxCoeff()));
18 }
m1
Matrix3d m1
Definition: IOFormat.cpp:2
m2
MatrixType m2(n_dims)
max
#define max(a, b)
Definition: datatypes.h:23
oomph::SarahBL::epsilon
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43

References oomph::SarahBL::epsilon, m1, m2(), and max.

Referenced by product().

◆ check_mismatched_product() [1/2]

template<typename LhsType , typename RhsType >
std::enable_if_t<RhsType::SizeAtCompileTime != Dynamic, void> check_mismatched_product ( LhsType &  ,
const RhsType &   
)
39  {}

◆ check_mismatched_product() [2/2]

template<typename LhsType , typename RhsType >
std::enable_if_t<RhsType::SizeAtCompileTime == Dynamic, void> check_mismatched_product ( LhsType &  lhs,
const RhsType &  rhs 
)
33  {
34  VERIFY_RAISES_ASSERT(lhs = rhs * rhs);
35 }
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:329

References VERIFY_RAISES_ASSERT.

Referenced by product().

◆ product()

template<typename MatrixType >
void product ( const MatrixType &  m)
42  {
43  /* this test covers the following files:
44  Identity.h Product.h
45  */
46  typedef typename MatrixType::Scalar Scalar;
47  typedef typename MatrixType::RealScalar RealScalar;
48  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RowVectorType;
49  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ColVectorType;
50  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType;
51  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> ColSquareMatrixType;
52  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
53  MatrixType::Flags & RowMajorBit ? ColMajor : RowMajor>
54  OtherMajorMatrixType;
55 
56  // We want a tighter epsilon for not-approx tests. Otherwise, for certain
57  // low-precision types (e.g. bfloat16), the bound ends up being relatively large
58  // (e.g. 0.12), causing flaky tests.
59  RealScalar not_approx_epsilon = RealScalar(0.1) * NumTraits<RealScalar>::dummy_precision();
60 
61  Index rows = m.rows();
62  Index cols = m.cols();
63 
64  // this test relies a lot on Random.h, and there's not much more that we can do
65  // to test it, hence I consider that we will have tested Random.h
66  MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), m3(rows, cols);
67  RowSquareMatrixType identity = RowSquareMatrixType::Identity(rows, rows),
68  square = RowSquareMatrixType::Random(rows, rows), res = RowSquareMatrixType::Random(rows, rows);
69  ColSquareMatrixType square2 = ColSquareMatrixType::Random(cols, cols), res2 = ColSquareMatrixType::Random(cols, cols);
70  RowVectorType v1 = RowVectorType::Random(rows);
71  ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols);
72 
73  // Prevent overflows for integer types.
74  if (Eigen::NumTraits<Scalar>::IsInteger) {
75  Scalar kMaxVal = Scalar(10000);
76  m1.array() = m1.array() - kMaxVal * (m1.array() / kMaxVal);
77  m2.array() = m2.array() - kMaxVal * (m2.array() / kMaxVal);
78  v1.array() = v1.array() - kMaxVal * (v1.array() / kMaxVal);
79  }
80 
81  OtherMajorMatrixType tm1 = m1;
82 
83  Scalar s1 = internal::random<Scalar>();
84 
85  Index r = internal::random<Index>(0, rows - 1), c = internal::random<Index>(0, cols - 1),
86  c2 = internal::random<Index>(0, cols - 1);
87 
88  // begin testing Product.h: only associativity for now
89  // (we use Transpose.h but this doesn't count as a test for it)
90  {
91  // Increase tolerance, since coefficients here can get relatively large.
92  RealScalar tol = RealScalar(2) * get_test_precision(m1);
93  VERIFY(verifyIsApprox((m1 * m1.transpose()) * m2, m1 * (m1.transpose() * m2), tol));
94  }
95  m3 = m1;
96  m3 *= m1.transpose() * m2;
97  VERIFY_IS_APPROX(m3, m1 * (m1.transpose() * m2));
98  VERIFY_IS_APPROX(m3, m1 * (m1.transpose() * m2));
99 
100  // continue testing Product.h: distributivity
101  VERIFY_IS_APPROX(square * (m1 + m2), square * m1 + square * m2);
102  VERIFY_IS_APPROX(square * (m1 - m2), square * m1 - square * m2);
103 
104  // continue testing Product.h: compatibility with ScalarMultiple.h
105  VERIFY_IS_APPROX(s1 * (square * m1), (s1 * square) * m1);
106  VERIFY_IS_APPROX(s1 * (square * m1), square * (m1 * s1));
107 
108  // test Product.h together with Identity.h
109  VERIFY_IS_APPROX(v1, identity * v1);
110  VERIFY_IS_APPROX(v1.transpose(), v1.transpose() * identity);
111  // again, test operator() to check const-qualification
112  VERIFY_IS_APPROX(MatrixType::Identity(rows, cols)(r, c), static_cast<Scalar>(r == c));
113 
114  if (rows != cols) {
115  check_mismatched_product(m3, m1);
116  }
117 
118  // test the previous tests were not screwed up because operator* returns 0
119  // (we use the more accurate default epsilon)
120  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
121  VERIFY(areNotApprox(m1.transpose() * m2, m2.transpose() * m1, not_approx_epsilon));
122  }
123 
124  // test optimized operator+= path
125  res = square;
126  res.noalias() += m1 * m2.transpose();
127  VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
128  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
129  VERIFY(areNotApprox(res, square + m2 * m1.transpose(), not_approx_epsilon));
130  }
131  vcres = vc2;
132  vcres.noalias() += m1.transpose() * v1;
133  VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1);
134 
135  // test optimized operator-= path
136  res = square;
137  res.noalias() -= m1 * m2.transpose();
138  VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
139  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
140  VERIFY(areNotApprox(res, square - m2 * m1.transpose(), not_approx_epsilon));
141  }
142  vcres = vc2;
143  vcres.noalias() -= m1.transpose() * v1;
144  VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1);
145 
146  // test scaled products
147  res = square;
148  res.noalias() = s1 * m1 * m2.transpose();
149  VERIFY_IS_APPROX(res, ((s1 * m1).eval() * m2.transpose()));
150  res = square;
151  res.noalias() += s1 * m1 * m2.transpose();
152  VERIFY_IS_APPROX(res, square + ((s1 * m1).eval() * m2.transpose()));
153  res = square;
154  res.noalias() -= s1 * m1 * m2.transpose();
155  VERIFY_IS_APPROX(res, square - ((s1 * m1).eval() * m2.transpose()));
156 
157  // test d ?= a+b*c rules
158  res.noalias() = square + m1 * m2.transpose();
159  VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
160  res.noalias() += square + m1 * m2.transpose();
161  VERIFY_IS_APPROX(res, Scalar(2) * (square + m1 * m2.transpose()));
162  res.noalias() -= square + m1 * m2.transpose();
163  VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
164 
165  // test d ?= a-b*c rules
166  res.noalias() = square - m1 * m2.transpose();
167  VERIFY_IS_APPROX(res, square - m1 * m2.transpose());
168  res.noalias() += square - m1 * m2.transpose();
169  VERIFY_IS_APPROX(res, Scalar(2) * (square - m1 * m2.transpose()));
170  res.noalias() -= square - m1 * m2.transpose();
171  VERIFY_IS_APPROX(res, square - m1 * m2.transpose());
172 
173  tm1 = m1;
174  VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1);
175  VERIFY_IS_APPROX(v1.transpose() * tm1, v1.transpose() * m1);
176 
177  // test submatrix and matrix/vector product
178  for (int i = 0; i < rows; ++i) res.row(i) = m1.row(i) * m2.transpose();
179  VERIFY_IS_APPROX(res, m1 * m2.transpose());
180  // the other way round:
181  for (int i = 0; i < rows; ++i) res.col(i) = m1 * m2.transpose().col(i);
182  VERIFY_IS_APPROX(res, m1 * m2.transpose());
183 
184  res2 = square2;
185  res2.noalias() += m1.transpose() * m2;
186  VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
187  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
188  VERIFY(areNotApprox(res2, square2 + m2.transpose() * m1, not_approx_epsilon));
189  }
190 
191  res.col(r).noalias() = square.adjoint() * square.col(r);
192  VERIFY_IS_APPROX(res.col(r), (square.adjoint() * square.col(r)).eval());
193  res.col(r).noalias() = square * square.col(r);
194  VERIFY_IS_APPROX(res.col(r), (square * square.col(r)).eval());
195 
196  // vector at runtime (see bug 1166)
197  {
198  RowSquareMatrixType ref(square);
199  ColSquareMatrixType ref2(square2);
200  ref = res = square;
201  VERIFY_IS_APPROX(res.block(0, 0, 1, rows).noalias() = m1.col(0).transpose() * square.transpose(),
202  (ref.row(0) = m1.col(0).transpose() * square.transpose()));
203  VERIFY_IS_APPROX(res.block(0, 0, 1, rows).noalias() = m1.block(0, 0, rows, 1).transpose() * square.transpose(),
204  (ref.row(0) = m1.col(0).transpose() * square.transpose()));
205  VERIFY_IS_APPROX(res.block(0, 0, 1, rows).noalias() = m1.col(0).transpose() * square,
206  (ref.row(0) = m1.col(0).transpose() * square));
207  VERIFY_IS_APPROX(res.block(0, 0, 1, rows).noalias() = m1.block(0, 0, rows, 1).transpose() * square,
208  (ref.row(0) = m1.col(0).transpose() * square));
209  ref2 = res2 = square2;
210  VERIFY_IS_APPROX(res2.block(0, 0, 1, cols).noalias() = m1.row(0) * square2.transpose(),
211  (ref2.row(0) = m1.row(0) * square2.transpose()));
212  VERIFY_IS_APPROX(res2.block(0, 0, 1, cols).noalias() = m1.block(0, 0, 1, cols) * square2.transpose(),
213  (ref2.row(0) = m1.row(0) * square2.transpose()));
214  VERIFY_IS_APPROX(res2.block(0, 0, 1, cols).noalias() = m1.row(0) * square2, (ref2.row(0) = m1.row(0) * square2));
215  VERIFY_IS_APPROX(res2.block(0, 0, 1, cols).noalias() = m1.block(0, 0, 1, cols) * square2,
216  (ref2.row(0) = m1.row(0) * square2));
217  }
218 
219  // vector.block() (see bug 1283)
220  {
221  RowVectorType w1(rows);
222  VERIFY_IS_APPROX(square * v1.block(0, 0, rows, 1), square * v1);
223  VERIFY_IS_APPROX(w1.noalias() = square * v1.block(0, 0, rows, 1), square * v1);
224  VERIFY_IS_APPROX(w1.block(0, 0, rows, 1).noalias() = square * v1.block(0, 0, rows, 1), square * v1);
225 
226  Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> w2(cols);
227  VERIFY_IS_APPROX(vc2.block(0, 0, cols, 1).transpose() * square2, vc2.transpose() * square2);
228  VERIFY_IS_APPROX(w2.noalias() = vc2.block(0, 0, cols, 1).transpose() * square2, vc2.transpose() * square2);
229  VERIFY_IS_APPROX(w2.block(0, 0, 1, cols).noalias() = vc2.block(0, 0, cols, 1).transpose() * square2,
230  vc2.transpose() * square2);
231 
232  vc2 = square2.block(0, 0, 1, cols).transpose();
233  VERIFY_IS_APPROX(square2.block(0, 0, 1, cols) * square2, vc2.transpose() * square2);
234  VERIFY_IS_APPROX(w2.noalias() = square2.block(0, 0, 1, cols) * square2, vc2.transpose() * square2);
235  VERIFY_IS_APPROX(w2.block(0, 0, 1, cols).noalias() = square2.block(0, 0, 1, cols) * square2,
236  vc2.transpose() * square2);
237 
238  vc2 = square2.block(0, 0, cols, 1);
239  VERIFY_IS_APPROX(square2.block(0, 0, cols, 1).transpose() * square2, vc2.transpose() * square2);
240  VERIFY_IS_APPROX(w2.noalias() = square2.block(0, 0, cols, 1).transpose() * square2, vc2.transpose() * square2);
241  VERIFY_IS_APPROX(w2.block(0, 0, 1, cols).noalias() = square2.block(0, 0, cols, 1).transpose() * square2,
242  vc2.transpose() * square2);
243  }
244 
245  // inner product
246  {
247  Scalar x = square2.row(c) * square2.col(c2);
248  VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum());
249  }
250 
251  // outer product
252  {
253  VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0, c, rows, 1) * m1.block(r, 0, 1, cols));
254  VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(),
255  m1.block(r, 0, 1, cols).transpose() * m1.block(0, c, rows, 1).transpose());
256  VERIFY_IS_APPROX(m1.block(0, c, rows, 1) * m1.row(r), m1.block(0, c, rows, 1) * m1.block(r, 0, 1, cols));
257  VERIFY_IS_APPROX(m1.col(c) * m1.block(r, 0, 1, cols), m1.block(0, c, rows, 1) * m1.block(r, 0, 1, cols));
258  VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0, 0, rows, 1) * m1.block(r, 0, 1, cols));
259  VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0, c, rows, 1) * m1.block(0, 0, 1, cols));
260  }
261 
262  // Aliasing
263  {
264  ColVectorType x(cols);
265  x.setRandom();
266  ColVectorType z(x);
267  ColVectorType y(cols);
268  y.setZero();
269  ColSquareMatrixType A(cols, cols);
270  A.setRandom();
271  // CwiseBinaryOp
272  VERIFY_IS_APPROX(x = y + A * x, A * z);
273  x = z;
274  VERIFY_IS_APPROX(x = y - A * x, A * (-z));
275  x = z;
276  // CwiseUnaryOp
277  VERIFY_IS_APPROX(x = Scalar(1.) * (A * x), A * z);
278  }
279 
280  // regression for blas_trais
281  {
282  // Increase test tolerance, since coefficients can get relatively large.
283  RealScalar tol = RealScalar(2) * get_test_precision(square);
284  VERIFY(
285  verifyIsApprox(square * (square * square).transpose(), square * square.transpose() * square.transpose(), tol));
286  VERIFY(verifyIsApprox(square * (-(square * square)), -square * square * square, tol));
287  VERIFY(verifyIsApprox(square * (s1 * (square * square)), s1 * square * square * square, tol));
288  VERIFY(
289  verifyIsApprox(square * (square * square).conjugate(), square * square.conjugate() * square.conjugate(), tol));
290  }
291 
292  // destination with a non-default inner-stride
293  // see bug 1741
294  if (!MatrixType::IsRowMajor) {
295  typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX;
296  MatrixX buffer(2 * rows, 2 * rows);
297  Map<RowSquareMatrixType, 0, Stride<Dynamic, 2> > map1(buffer.data(), rows, rows, Stride<Dynamic, 2>(2 * rows, 2));
298  buffer.setZero();
299  VERIFY_IS_APPROX(map1 = m1 * m2.transpose(), (m1 * m2.transpose()).eval());
300  buffer.setZero();
301  VERIFY_IS_APPROX(map1.noalias() = m1 * m2.transpose(), (m1 * m2.transpose()).eval());
302  buffer.setZero();
303  VERIFY_IS_APPROX(map1.noalias() += m1 * m2.transpose(), (m1 * m2.transpose()).eval());
304  }
305 }
get_test_precision
AnnoyingScalar get_test_precision(const AnnoyingScalar &)
Definition: AnnoyingScalar.h:190
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
v1
M1<< 1, 2, 3, 4, 5, 6, 7, 8, 9;Map< RowVectorXf > v1(M1.data(), M1.size())
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:47
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:46
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
Eigen::PlainObjectBase::setRandom
Derived & setRandom(Index size)
Definition: Random.h:147
Eigen::Stride
Holds strides information for Map.
Definition: Stride.h:55
min
#define min(a, b)
Definition: datatypes.h:22
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:318
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:320
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:70
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:13
y
Scalar * y
Definition: level1_cplx_impl.h:128
m
int * m
Definition: level2_cplx_impl.h:294
VERIFY
#define VERIFY(a)
Definition: main.h:362
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:83
Eigen::square
squared absolute sa ArrayBase::abs2 DOXCOMMA MatrixBase::cwiseAbs2 square(power 2)
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose
void transpose()
Definition: skew_symmetric_matrix3.cpp:135
calibrate.c
int c
Definition: calibrate.py:100
plotDoE.x
list x
Definition: plotDoE.py:28
check_mismatched_product
std::enable_if_t< RhsType::SizeAtCompileTime==Dynamic, void > check_mismatched_product(LhsType &lhs, const RhsType &rhs)
Definition: product.h:32
areNotApprox
bool areNotApprox(const MatrixBase< Derived1 > &m1, const MatrixBase< Derived2 > &m2, typename Derived1::RealScalar epsilon=NumTraits< typename Derived1::RealScalar >::dummy_precision())
Definition: product.h:14
verifyIsApprox
bool verifyIsApprox(const Type1 &a, const Type2 &b, Tol tol)
Definition: product.h:22
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:47
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:217

References areNotApprox(), calibrate::c, check_mismatched_product(), Eigen::ColMajor, cols, eval(), get_test_precision(), i, m, m1, m2(), min, UniformPSDSelfTest::r, res, Eigen::RowMajor, Eigen::RowMajorBit, rows, Eigen::PlainObjectBase< Derived >::setRandom(), Eigen::square(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), v1(), VERIFY, VERIFY_IS_APPROX, verifyIsApprox(), plotDoE::x, and y.

Referenced by EIGEN_DECLARE_TEST(), oomph::AxisymmetricNavierStokesEquations::fill_in_contribution_to_hessian_vector_products(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_hessian_vector_products(), oomph::Problem::get_hessian_vector_products(), oomph::AssemblyHandler::get_hessian_vector_products(), oomph::PitchForkHandler::get_hessian_vector_products(), oomph::GeneralisedElement::get_hessian_vector_products(), and Eigen::SparseVector< Scalar_, Options_, StorageIndex_ >::operator=().

◆ verifyIsApprox()

template<typename Type1 , typename Type2 , typename Tol >
bool verifyIsApprox ( const Type1 &  a,
const Type2 &  b,
Tol  tol 
)
inline
22  {
23  bool ret = a.isApprox(b, tol);
24  if (!ret) {
25  std::cerr << "Difference too large wrt tolerance " << tol << ", relative error is: " << test_relative_error(a, b)
26  << std::endl;
27  }
28  return ret;
29 }
test_relative_error
AnnoyingScalar test_relative_error(const AnnoyingScalar &a, const AnnoyingScalar &b)
Definition: AnnoyingScalar.h:192
b
Scalar * b
Definition: benchVecAdd.cpp:17
ret
Eigen::DenseIndex ret
Definition: level1_cplx_impl.h:43
a
const Scalar * a
Definition: level2_cplx_impl.h:32

References a, b, ret, and test_relative_error().

Referenced by product().