product_small.cpp 12 KB

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