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 #include "product.h"
12 #include <Eigen/LU>
13
14 // regression test for bug 447
15 template<int>
product1x1()16 void product1x1()
17 {
18 Matrix<float,1,3> matAstatic;
19 Matrix<float,3,1> matBstatic;
20 matAstatic.setRandom();
21 matBstatic.setRandom();
22 VERIFY_IS_APPROX( (matAstatic * matBstatic).coeff(0,0),
23 matAstatic.cwiseProduct(matBstatic.transpose()).sum() );
24
25 MatrixXf matAdynamic(1,3);
26 MatrixXf matBdynamic(3,1);
27 matAdynamic.setRandom();
28 matBdynamic.setRandom();
29 VERIFY_IS_APPROX( (matAdynamic * matBdynamic).coeff(0,0),
30 matAdynamic.cwiseProduct(matBdynamic.transpose()).sum() );
31 }
32
33 template<typename TC, typename TA, typename TB>
ref_prod(TC & C,const TA & A,const TB & B)34 const TC& ref_prod(TC &C, const TA &A, const TB &B)
35 {
36 for(Index i=0;i<C.rows();++i)
37 for(Index j=0;j<C.cols();++j)
38 for(Index k=0;k<A.cols();++k)
39 C.coeffRef(i,j) += A.coeff(i,k) * B.coeff(k,j);
40 return C;
41 }
42
43 template<typename T, int Rows, int Cols, int Depth, int OC, int OA, int OB>
44 typename internal::enable_if<! ( (Rows ==1&&Depth!=1&&OA==ColMajor)
45 || (Depth==1&&Rows !=1&&OA==RowMajor)
46 || (Cols ==1&&Depth!=1&&OB==RowMajor)
47 || (Depth==1&&Cols !=1&&OB==ColMajor)
48 || (Rows ==1&&Cols !=1&&OC==ColMajor)
49 || (Cols ==1&&Rows !=1&&OC==RowMajor)),void>::type
test_lazy_single(int rows,int cols,int depth)50 test_lazy_single(int rows, int cols, int depth)
51 {
52 Matrix<T,Rows,Depth,OA> A(rows,depth); A.setRandom();
53 Matrix<T,Depth,Cols,OB> B(depth,cols); B.setRandom();
54 Matrix<T,Rows,Cols,OC> C(rows,cols); C.setRandom();
55 Matrix<T,Rows,Cols,OC> D(C);
56 VERIFY_IS_APPROX(C+=A.lazyProduct(B), ref_prod(D,A,B));
57 }
58
test_dynamic_bool()59 void test_dynamic_bool()
60 {
61 int rows = internal::random<int>(1,64);
62 int cols = internal::random<int>(1,64);
63 int depth = internal::random<int>(1,65);
64
65 typedef Matrix<bool,Dynamic,Dynamic> MatrixX;
66 MatrixX A(rows,depth); A.setRandom();
67 MatrixX B(depth,cols); B.setRandom();
68 MatrixX C(rows,cols); C.setRandom();
69 MatrixX D(C);
70 for(Index i=0;i<C.rows();++i)
71 for(Index j=0;j<C.cols();++j)
72 for(Index k=0;k<A.cols();++k)
73 D.coeffRef(i,j) |= A.coeff(i,k) & B.coeff(k,j);
74 C += A * B;
75 VERIFY_IS_EQUAL(C, D);
76
77 MatrixX E = B.transpose();
78 for(Index i=0;i<B.rows();++i)
79 for(Index j=0;j<B.cols();++j)
80 VERIFY_IS_EQUAL(B(i,j), E(j,i));
81 }
82
83 template<typename T, int Rows, int Cols, int Depth, int OC, int OA, int OB>
84 typename internal::enable_if< ( (Rows ==1&&Depth!=1&&OA==ColMajor)
85 || (Depth==1&&Rows !=1&&OA==RowMajor)
86 || (Cols ==1&&Depth!=1&&OB==RowMajor)
87 || (Depth==1&&Cols !=1&&OB==ColMajor)
88 || (Rows ==1&&Cols !=1&&OC==ColMajor)
89 || (Cols ==1&&Rows !=1&&OC==RowMajor)),void>::type
test_lazy_single(int,int,int)90 test_lazy_single(int, int, int)
91 {
92 }
93
94 template<typename T, int Rows, int Cols, int Depth>
test_lazy_all_layout(int rows=Rows,int cols=Cols,int depth=Depth)95 void test_lazy_all_layout(int rows=Rows, int cols=Cols, int depth=Depth)
96 {
97 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,ColMajor,ColMajor,ColMajor>(rows,cols,depth) ));
98 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,RowMajor,ColMajor,ColMajor>(rows,cols,depth) ));
99 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,ColMajor,RowMajor,ColMajor>(rows,cols,depth) ));
100 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,RowMajor,RowMajor,ColMajor>(rows,cols,depth) ));
101 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,ColMajor,ColMajor,RowMajor>(rows,cols,depth) ));
102 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,RowMajor,ColMajor,RowMajor>(rows,cols,depth) ));
103 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,ColMajor,RowMajor,RowMajor>(rows,cols,depth) ));
104 CALL_SUBTEST(( test_lazy_single<T,Rows,Cols,Depth,RowMajor,RowMajor,RowMajor>(rows,cols,depth) ));
105 }
106
107 template<typename T>
test_lazy_l1()108 void test_lazy_l1()
109 {
110 int rows = internal::random<int>(1,12);
111 int cols = internal::random<int>(1,12);
112 int depth = internal::random<int>(1,12);
113
114 // Inner
115 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,1>() ));
116 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,2>() ));
117 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,3>() ));
118 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,8>() ));
119 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,9>() ));
120 CALL_SUBTEST(( test_lazy_all_layout<T,1,1,-1>(1,1,depth) ));
121
122 // Outer
123 CALL_SUBTEST(( test_lazy_all_layout<T,2,1,1>() ));
124 CALL_SUBTEST(( test_lazy_all_layout<T,1,2,1>() ));
125 CALL_SUBTEST(( test_lazy_all_layout<T,2,2,1>() ));
126 CALL_SUBTEST(( test_lazy_all_layout<T,3,3,1>() ));
127 CALL_SUBTEST(( test_lazy_all_layout<T,4,4,1>() ));
128 CALL_SUBTEST(( test_lazy_all_layout<T,4,8,1>() ));
129 CALL_SUBTEST(( test_lazy_all_layout<T,4,-1,1>(4,cols) ));
130 CALL_SUBTEST(( test_lazy_all_layout<T,7,-1,1>(7,cols) ));
131 CALL_SUBTEST(( test_lazy_all_layout<T,-1,8,1>(rows) ));
132 CALL_SUBTEST(( test_lazy_all_layout<T,-1,3,1>(rows) ));
133 CALL_SUBTEST(( test_lazy_all_layout<T,-1,-1,1>(rows,cols) ));
134 }
135
136 template<typename T>
test_lazy_l2()137 void test_lazy_l2()
138 {
139 int rows = internal::random<int>(1,12);
140 int cols = internal::random<int>(1,12);
141 int depth = internal::random<int>(1,12);
142
143 // mat-vec
144 CALL_SUBTEST(( test_lazy_all_layout<T,2,1,2>() ));
145 CALL_SUBTEST(( test_lazy_all_layout<T,2,1,4>() ));
146 CALL_SUBTEST(( test_lazy_all_layout<T,4,1,2>() ));
147 CALL_SUBTEST(( test_lazy_all_layout<T,4,1,4>() ));
148 CALL_SUBTEST(( test_lazy_all_layout<T,5,1,4>() ));
149 CALL_SUBTEST(( test_lazy_all_layout<T,4,1,5>() ));
150 CALL_SUBTEST(( test_lazy_all_layout<T,4,1,6>() ));
151 CALL_SUBTEST(( test_lazy_all_layout<T,6,1,4>() ));
152 CALL_SUBTEST(( test_lazy_all_layout<T,8,1,8>() ));
153 CALL_SUBTEST(( test_lazy_all_layout<T,-1,1,4>(rows) ));
154 CALL_SUBTEST(( test_lazy_all_layout<T,4,1,-1>(4,1,depth) ));
155 CALL_SUBTEST(( test_lazy_all_layout<T,-1,1,-1>(rows,1,depth) ));
156
157 // vec-mat
158 CALL_SUBTEST(( test_lazy_all_layout<T,1,2,2>() ));
159 CALL_SUBTEST(( test_lazy_all_layout<T,1,2,4>() ));
160 CALL_SUBTEST(( test_lazy_all_layout<T,1,4,2>() ));
161 CALL_SUBTEST(( test_lazy_all_layout<T,1,4,4>() ));
162 CALL_SUBTEST(( test_lazy_all_layout<T,1,5,4>() ));
163 CALL_SUBTEST(( test_lazy_all_layout<T,1,4,5>() ));
164 CALL_SUBTEST(( test_lazy_all_layout<T,1,4,6>() ));
165 CALL_SUBTEST(( test_lazy_all_layout<T,1,6,4>() ));
166 CALL_SUBTEST(( test_lazy_all_layout<T,1,8,8>() ));
167 CALL_SUBTEST(( test_lazy_all_layout<T,1,-1, 4>(1,cols) ));
168 CALL_SUBTEST(( test_lazy_all_layout<T,1, 4,-1>(1,4,depth) ));
169 CALL_SUBTEST(( test_lazy_all_layout<T,1,-1,-1>(1,cols,depth) ));
170 }
171
172 template<typename T>
test_lazy_l3()173 void test_lazy_l3()
174 {
175 int rows = internal::random<int>(1,12);
176 int cols = internal::random<int>(1,12);
177 int depth = internal::random<int>(1,12);
178 // mat-mat
179 CALL_SUBTEST(( test_lazy_all_layout<T,2,4,2>() ));
180 CALL_SUBTEST(( test_lazy_all_layout<T,2,6,4>() ));
181 CALL_SUBTEST(( test_lazy_all_layout<T,4,3,2>() ));
182 CALL_SUBTEST(( test_lazy_all_layout<T,4,8,4>() ));
183 CALL_SUBTEST(( test_lazy_all_layout<T,5,6,4>() ));
184 CALL_SUBTEST(( test_lazy_all_layout<T,4,2,5>() ));
185 CALL_SUBTEST(( test_lazy_all_layout<T,4,7,6>() ));
186 CALL_SUBTEST(( test_lazy_all_layout<T,6,8,4>() ));
187 CALL_SUBTEST(( test_lazy_all_layout<T,8,3,8>() ));
188 CALL_SUBTEST(( test_lazy_all_layout<T,-1,6,4>(rows) ));
189 CALL_SUBTEST(( test_lazy_all_layout<T,4,3,-1>(4,3,depth) ));
190 CALL_SUBTEST(( test_lazy_all_layout<T,-1,6,-1>(rows,6,depth) ));
191 CALL_SUBTEST(( test_lazy_all_layout<T,8,2,2>() ));
192 CALL_SUBTEST(( test_lazy_all_layout<T,5,2,4>() ));
193 CALL_SUBTEST(( test_lazy_all_layout<T,4,4,2>() ));
194 CALL_SUBTEST(( test_lazy_all_layout<T,8,4,4>() ));
195 CALL_SUBTEST(( test_lazy_all_layout<T,6,5,4>() ));
196 CALL_SUBTEST(( test_lazy_all_layout<T,4,4,5>() ));
197 CALL_SUBTEST(( test_lazy_all_layout<T,3,4,6>() ));
198 CALL_SUBTEST(( test_lazy_all_layout<T,2,6,4>() ));
199 CALL_SUBTEST(( test_lazy_all_layout<T,7,8,8>() ));
200 CALL_SUBTEST(( test_lazy_all_layout<T,8,-1, 4>(8,cols) ));
201 CALL_SUBTEST(( test_lazy_all_layout<T,3, 4,-1>(3,4,depth) ));
202 CALL_SUBTEST(( test_lazy_all_layout<T,4,-1,-1>(4,cols,depth) ));
203 }
204
205 template<typename T,int N,int M,int K>
test_linear_but_not_vectorizable()206 void test_linear_but_not_vectorizable()
207 {
208 // Check tricky cases for which the result of the product is a vector and thus must exhibit the LinearBit flag,
209 // but is not vectorizable along the linear dimension.
210 Index n = N==Dynamic ? internal::random<Index>(1,32) : N;
211 Index m = M==Dynamic ? internal::random<Index>(1,32) : M;
212 Index k = K==Dynamic ? internal::random<Index>(1,32) : K;
213
214 {
215 Matrix<T,N,M+1> A; A.setRandom(n,m+1);
216 Matrix<T,M*2,K> B; B.setRandom(m*2,k);
217 Matrix<T,1,K> C;
218 Matrix<T,1,K> R;
219
220 C.noalias() = A.template topLeftCorner<1,M>() * (B.template topRows<M>()+B.template bottomRows<M>());
221 R.noalias() = A.template topLeftCorner<1,M>() * (B.template topRows<M>()+B.template bottomRows<M>()).eval();
222 VERIFY_IS_APPROX(C,R);
223 }
224
225 {
226 Matrix<T,M+1,N,RowMajor> A; A.setRandom(m+1,n);
227 Matrix<T,K,M*2,RowMajor> B; B.setRandom(k,m*2);
228 Matrix<T,K,1> C;
229 Matrix<T,K,1> R;
230
231 C.noalias() = (B.template leftCols<M>()+B.template rightCols<M>()) * A.template topLeftCorner<M,1>();
232 R.noalias() = (B.template leftCols<M>()+B.template rightCols<M>()).eval() * A.template topLeftCorner<M,1>();
233 VERIFY_IS_APPROX(C,R);
234 }
235 }
236
237 template<int Rows>
bug_1311()238 void bug_1311()
239 {
240 Matrix< double, Rows, 2 > A; A.setRandom();
241 Vector2d b = Vector2d::Random() ;
242 Matrix<double,Rows,1> res;
243 res.noalias() = 1. * (A * b);
244 VERIFY_IS_APPROX(res, A*b);
245 res.noalias() = 1.*A * b;
246 VERIFY_IS_APPROX(res, A*b);
247 res.noalias() = (1.*A).lazyProduct(b);
248 VERIFY_IS_APPROX(res, A*b);
249 res.noalias() = (1.*A).lazyProduct(1.*b);
250 VERIFY_IS_APPROX(res, A*b);
251 res.noalias() = (A).lazyProduct(1.*b);
252 VERIFY_IS_APPROX(res, A*b);
253 }
254
255 template<int>
product_small_regressions()256 void product_small_regressions()
257 {
258 {
259 // test compilation of (outer_product) * vector
260 Vector3f v = Vector3f::Random();
261 VERIFY_IS_APPROX( (v * v.transpose()) * v, (v * v.transpose()).eval() * v);
262 }
263
264 {
265 // regression test for pull-request #93
266 Eigen::Matrix<double, 1, 1> A; A.setRandom();
267 Eigen::Matrix<double, 18, 1> B; B.setRandom();
268 Eigen::Matrix<double, 1, 18> C; C.setRandom();
269 VERIFY_IS_APPROX(B * A.inverse(), B * A.inverse()[0]);
270 VERIFY_IS_APPROX(A.inverse() * C, A.inverse()[0] * C);
271 }
272
273 {
274 Eigen::Matrix<double, 10, 10> A, B, C;
275 A.setRandom();
276 C = A;
277 for(int k=0; k<79; ++k)
278 C = C * A;
279 B.noalias() = (((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)) * ((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)))
280 * (((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)) * ((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)));
281 VERIFY_IS_APPROX(B,C);
282 }
283 }
284
EIGEN_DECLARE_TEST(product_small)285 EIGEN_DECLARE_TEST(product_small)
286 {
287 for(int i = 0; i < g_repeat; i++) {
288 CALL_SUBTEST_1( product(Matrix<float, 3, 2>()) );
289 CALL_SUBTEST_2( product(Matrix<int, 3, 17>()) );
290 CALL_SUBTEST_8( product(Matrix<double, 3, 17>()) );
291 CALL_SUBTEST_3( product(Matrix3d()) );
292 CALL_SUBTEST_4( product(Matrix4d()) );
293 CALL_SUBTEST_5( product(Matrix4f()) );
294 CALL_SUBTEST_6( product1x1<0>() );
295
296 CALL_SUBTEST_11( test_lazy_l1<float>() );
297 CALL_SUBTEST_12( test_lazy_l2<float>() );
298 CALL_SUBTEST_13( test_lazy_l3<float>() );
299
300 CALL_SUBTEST_21( test_lazy_l1<double>() );
301 CALL_SUBTEST_22( test_lazy_l2<double>() );
302 CALL_SUBTEST_23( test_lazy_l3<double>() );
303
304 CALL_SUBTEST_31( test_lazy_l1<std::complex<float> >() );
305 CALL_SUBTEST_32( test_lazy_l2<std::complex<float> >() );
306 CALL_SUBTEST_33( test_lazy_l3<std::complex<float> >() );
307
308 CALL_SUBTEST_41( test_lazy_l1<std::complex<double> >() );
309 CALL_SUBTEST_42( test_lazy_l2<std::complex<double> >() );
310 CALL_SUBTEST_43( test_lazy_l3<std::complex<double> >() );
311
312 CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,2,1,Dynamic>() ));
313 CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,3,1,Dynamic>() ));
314 CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,2,1,16>() ));
315
316 CALL_SUBTEST_6( bug_1311<3>() );
317 CALL_SUBTEST_6( bug_1311<5>() );
318
319 CALL_SUBTEST_9( test_dynamic_bool() );
320 }
321
322 CALL_SUBTEST_6( product_small_regressions<0>() );
323 }
324