Skip to content

Commit aafc9dc

Browse files
Eggy115Eggy115
authored
Add files via upload
1 parent ad1267d commit aafc9dc

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

55 files changed

+8637
-0
lines changed

‎adjoint.cpp

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
// This file is part of Eigen, a lightweight C++ template library
2+
// for linear algebra.
3+
//
4+
// Copyright (C) 2006-2008 Benoit Jacob <[email protected]>
5+
//
6+
// This Source Code Form is subject to the terms of the Mozilla
7+
// Public License v. 2.0. If a copy of the MPL was not distributed
8+
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9+
10+
#define EIGEN_NO_STATIC_ASSERT
11+
12+
#include "main.h"
13+
14+
template<bool IsInteger> struct adjoint_specific;
15+
16+
template<> struct adjoint_specific<true> {
17+
template<typename Vec, typename Mat, typename Scalar>
18+
static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) {
19+
VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), numext::conj(s1) * v1.dot(v3) + numext::conj(s2) * v2.dot(v3), 0));
20+
VERIFY(test_isApproxWithRef(v3.dot(s1 * v1 + s2 * v2), s1*v3.dot(v1)+s2*v3.dot(v2), 0));
21+
22+
// check compatibility of dot and adjoint
23+
VERIFY(test_isApproxWithRef(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), 0));
24+
}
25+
};
26+
27+
template<> struct adjoint_specific<false> {
28+
template<typename Vec, typename Mat, typename Scalar>
29+
static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) {
30+
typedef typename NumTraits<Scalar>::Real RealScalar;
31+
using std::abs;
32+
33+
RealScalar ref = NumTraits<Scalar>::IsInteger ? RealScalar(0) : (std::max)((s1 * v1 + s2 * v2).norm(),v3.norm());
34+
VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), numext::conj(s1) * v1.dot(v3) + numext::conj(s2) * v2.dot(v3), ref));
35+
VERIFY(test_isApproxWithRef(v3.dot(s1 * v1 + s2 * v2), s1*v3.dot(v1)+s2*v3.dot(v2), ref));
36+
37+
VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm());
38+
// check normalized() and normalize()
39+
VERIFY_IS_APPROX(v1, v1.norm() * v1.normalized());
40+
v3 = v1;
41+
v3.normalize();
42+
VERIFY_IS_APPROX(v1, v1.norm() * v3);
43+
VERIFY_IS_APPROX(v3, v1.normalized());
44+
VERIFY_IS_APPROX(v3.norm(), RealScalar(1));
45+
46+
// check null inputs
47+
VERIFY_IS_APPROX((v1*0).normalized(), (v1*0));
48+
#if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE)
49+
RealScalar very_small = (std::numeric_limits<RealScalar>::min)();
50+
VERIFY( (v1*very_small).norm() == 0 );
51+
VERIFY_IS_APPROX((v1*very_small).normalized(), (v1*very_small));
52+
v3 = v1*very_small;
53+
v3.normalize();
54+
VERIFY_IS_APPROX(v3, (v1*very_small));
55+
#endif
56+
57+
// check compatibility of dot and adjoint
58+
ref = NumTraits<Scalar>::IsInteger ? 0 : (std::max)((std::max)(v1.norm(),v2.norm()),(std::max)((square * v2).norm(),(square.adjoint() * v1).norm()));
59+
VERIFY(internal::isMuchSmallerThan(abs(v1.dot(square * v2) - (square.adjoint() * v1).dot(v2)), ref, test_precision<Scalar>()));
60+
61+
// check that Random().normalized() works: tricky as the random xpr must be evaluated by
62+
// normalized() in order to produce a consistent result.
63+
VERIFY_IS_APPROX(Vec::Random(v1.size()).normalized().norm(), RealScalar(1));
64+
}
65+
};
66+
67+
template<typename MatrixType> void adjoint(const MatrixType& m)
68+
{
69+
/* this test covers the following files:
70+
Transpose.h Conjugate.h Dot.h
71+
*/
72+
using std::abs;
73+
typedef typename MatrixType::Index Index;
74+
typedef typename MatrixType::Scalar Scalar;
75+
typedef typename NumTraits<Scalar>::Real RealScalar;
76+
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
77+
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
78+
const Index PacketSize = internal::packet_traits<Scalar>::size;
79+
80+
Index rows = m.rows();
81+
Index cols = m.cols();
82+
83+
MatrixType m1 = MatrixType::Random(rows, cols),
84+
m2 = MatrixType::Random(rows, cols),
85+
m3(rows, cols),
86+
square = SquareMatrixType::Random(rows, rows);
87+
VectorType v1 = VectorType::Random(rows),
88+
v2 = VectorType::Random(rows),
89+
v3 = VectorType::Random(rows),
90+
vzero = VectorType::Zero(rows);
91+
92+
Scalar s1 = internal::random<Scalar>(),
93+
s2 = internal::random<Scalar>();
94+
95+
// check basic compatibility of adjoint, transpose, conjugate
96+
VERIFY_IS_APPROX(m1.transpose().conjugate().adjoint(), m1);
97+
VERIFY_IS_APPROX(m1.adjoint().conjugate().transpose(), m1);
98+
99+
// check multiplicative behavior
100+
VERIFY_IS_APPROX((m1.adjoint() * m2).adjoint(), m2.adjoint() * m1);
101+
VERIFY_IS_APPROX((s1 * m1).adjoint(), numext::conj(s1) * m1.adjoint());
102+
103+
// check basic properties of dot, squaredNorm
104+
VERIFY_IS_APPROX(numext::conj(v1.dot(v2)), v2.dot(v1));
105+
VERIFY_IS_APPROX(numext::real(v1.dot(v1)), v1.squaredNorm());
106+
107+
adjoint_specific<NumTraits<Scalar>::IsInteger>::run(v1, v2, v3, square, s1, s2);
108+
109+
VERIFY_IS_MUCH_SMALLER_THAN(abs(vzero.dot(v1)), static_cast<RealScalar>(1));
110+
111+
// like in testBasicStuff, test operator() to check const-qualification
112+
Index r = internal::random<Index>(0, rows-1),
113+
c = internal::random<Index>(0, cols-1);
114+
VERIFY_IS_APPROX(m1.conjugate()(r,c), numext::conj(m1(r,c)));
115+
VERIFY_IS_APPROX(m1.adjoint()(c,r), numext::conj(m1(r,c)));
116+
117+
// check inplace transpose
118+
m3 = m1;
119+
m3.transposeInPlace();
120+
VERIFY_IS_APPROX(m3,m1.transpose());
121+
m3.transposeInPlace();
122+
VERIFY_IS_APPROX(m3,m1);
123+
124+
if(PacketSize<m3.rows() && PacketSize<m3.cols())
125+
{
126+
m3 = m1;
127+
Index i = internal::random<Index>(0,m3.rows()-PacketSize);
128+
Index j = internal::random<Index>(0,m3.cols()-PacketSize);
129+
m3.template block<PacketSize,PacketSize>(i,j).transposeInPlace();
130+
VERIFY_IS_APPROX( (m3.template block<PacketSize,PacketSize>(i,j)), (m1.template block<PacketSize,PacketSize>(i,j).transpose()) );
131+
m3.template block<PacketSize,PacketSize>(i,j).transposeInPlace();
132+
VERIFY_IS_APPROX(m3,m1);
133+
}
134+
135+
// check inplace adjoint
136+
m3 = m1;
137+
m3.adjointInPlace();
138+
VERIFY_IS_APPROX(m3,m1.adjoint());
139+
m3.transposeInPlace();
140+
VERIFY_IS_APPROX(m3,m1.conjugate());
141+
142+
// check mixed dot product
143+
typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
144+
RealVectorType rv1 = RealVectorType::Random(rows);
145+
VERIFY_IS_APPROX(v1.dot(rv1.template cast<Scalar>()), v1.dot(rv1));
146+
VERIFY_IS_APPROX(rv1.template cast<Scalar>().dot(v1), rv1.dot(v1));
147+
}
148+
149+
void test_adjoint()
150+
{
151+
for(int i = 0; i < g_repeat; i++) {
152+
CALL_SUBTEST_1( adjoint(Matrix<float, 1, 1>()) );
153+
CALL_SUBTEST_2( adjoint(Matrix3d()) );
154+
CALL_SUBTEST_3( adjoint(Matrix4f()) );
155+
156+
CALL_SUBTEST_4( adjoint(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
157+
CALL_SUBTEST_5( adjoint(MatrixXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
158+
CALL_SUBTEST_6( adjoint(MatrixXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
159+
160+
// Complement for 128 bits vectorization:
161+
CALL_SUBTEST_8( adjoint(Matrix2d()) );
162+
CALL_SUBTEST_9( adjoint(Matrix<int,4,4>()) );
163+
164+
// 256 bits vectorization:
165+
CALL_SUBTEST_10( adjoint(Matrix<float,8,8>()) );
166+
CALL_SUBTEST_11( adjoint(Matrix<double,4,4>()) );
167+
CALL_SUBTEST_12( adjoint(Matrix<int,8,8>()) );
168+
}
169+
// test a large static matrix only once
170+
CALL_SUBTEST_7( adjoint(Matrix<float, 100, 100>()) );
171+
172+
#ifdef EIGEN_TEST_PART_13
173+
{
174+
MatrixXcf a(10,10), b(10,10);
175+
VERIFY_RAISES_ASSERT(a = a.transpose());
176+
VERIFY_RAISES_ASSERT(a = a.transpose() + b);
177+
VERIFY_RAISES_ASSERT(a = b + a.transpose());
178+
VERIFY_RAISES_ASSERT(a = a.conjugate().transpose());
179+
VERIFY_RAISES_ASSERT(a = a.adjoint());
180+
VERIFY_RAISES_ASSERT(a = a.adjoint() + b);
181+
VERIFY_RAISES_ASSERT(a = b + a.adjoint());
182+
183+
// no assertion should be triggered for these cases:
184+
a.transpose() = a.transpose();
185+
a.transpose() += a.transpose();
186+
a.transpose() += a.transpose() + b;
187+
a.transpose() = a.adjoint();
188+
a.transpose() += a.adjoint();
189+
a.transpose() += a.adjoint() + b;
190+
191+
// regression tests for check_for_aliasing
192+
MatrixXd c(10,10);
193+
c = 1.0 * MatrixXd::Ones(10,10) + c;
194+
c = MatrixXd::Ones(10,10) * 1.0 + c;
195+
c = c + MatrixXd::Ones(10,10) .cwiseProduct( MatrixXd::Zero(10,10) );
196+
c = MatrixXd::Ones(10,10) * MatrixXd::Zero(10,10);
197+
}
198+
#endif
199+
}
200+

0 commit comments

Comments
 (0)