diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp index e117f6931..9516dd07c 100644 --- a/test/product_trmm.cpp +++ b/test/product_trmm.cpp @@ -24,70 +24,97 @@ #include "main.h" -template void trmm(int size,int /*othersize*/) +template +void trmm(int rows=internal::random(1,320), int cols=internal::random(1,320), int otherCols = OtherCols==Dynamic?internal::random(1,320):OtherCols) { + typedef typename NumTraits::Real RealScalar; - typedef Matrix MatrixColMaj; - typedef Matrix MatrixRowMaj; + typedef Matrix TriMatrix; + typedef Matrix OnTheRight; + typedef Matrix OnTheLeft; + + typedef Matrix ResXS; + typedef Matrix ResSX; - DenseIndex rows = size; - DenseIndex cols = internal::random(1,size); - - MatrixColMaj triV(rows,cols), triH(cols,rows), upTri(cols,rows), loTri(rows,cols), - unitUpTri(cols,rows), unitLoTri(rows,cols), strictlyUpTri(cols,rows), strictlyLoTri(rows,cols); - MatrixColMaj ge1(rows,cols), ge2(cols,rows), ge3; - MatrixRowMaj rge3; + TriMatrix mat(rows,cols), tri(rows,cols), triTr(cols,rows); + + OnTheRight ge_right(cols,otherCols); + OnTheLeft ge_left(otherCols,rows); + ResSX ge_sx; + ResXS ge_xs; Scalar s1 = internal::random(), s2 = internal::random(); - triV.setRandom(); - triH.setRandom(); - loTri = triV.template triangularView(); - upTri = triH.template triangularView(); - unitLoTri = triV.template triangularView(); - unitUpTri = triH.template triangularView(); - strictlyLoTri = triV.template triangularView(); - strictlyUpTri = triH.template triangularView(); - ge1.setRandom(); - ge2.setRandom(); + mat.setRandom(); + tri = mat.template triangularView(); + triTr = mat.transpose().template triangularView(); + ge_right.setRandom(); + ge_left.setRandom(); - VERIFY_IS_APPROX( ge3 = triV.template triangularView() * ge2, loTri * ge2); - VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView(), ge2 * loTri); - VERIFY_IS_APPROX( ge3 = triH.template triangularView() * ge1, upTri * ge1); - VERIFY_IS_APPROX( ge3 = ge1 * triH.template triangularView(), ge1 * upTri); - VERIFY_IS_APPROX( ge3 = (s1*triV.adjoint()).template triangularView() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1)); - VERIFY_IS_APPROX( ge3 = ge1 * triV.adjoint().template triangularView(), ge1 * loTri.adjoint()); - VERIFY_IS_APPROX( ge3 = triH.adjoint().template triangularView() * ge2, upTri.adjoint() * ge2); - VERIFY_IS_APPROX( ge3 = ge2 * triH.adjoint().template triangularView(), ge2 * upTri.adjoint()); - VERIFY_IS_APPROX( ge3 = triV.template triangularView() * ge1.adjoint(), loTri * ge1.adjoint()); - VERIFY_IS_APPROX( ge3 = ge1.adjoint() * triV.template triangularView(), ge1.adjoint() * loTri); - VERIFY_IS_APPROX( ge3 = triH.template triangularView() * ge2.adjoint(), upTri * ge2.adjoint()); - VERIFY_IS_APPROX(rge3.noalias() = triH.template triangularView() * ge2.adjoint(), upTri * ge2.adjoint()); - VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView() * ge2.adjoint(), internal::conj(s1) * loTri.adjoint() * ge2.adjoint()); - VERIFY_IS_APPROX(rge3.noalias() = triV.adjoint().template triangularView() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); - VERIFY_IS_APPROX( ge3 = triH.adjoint().template triangularView() * ge1.adjoint(), upTri.adjoint() * ge1.adjoint()); - VERIFY_IS_APPROX(rge3.noalias() = triH.adjoint().template triangularView() * ge1.adjoint(), upTri.adjoint() * ge1.adjoint()); - - VERIFY_IS_APPROX( ge3 = triV.template triangularView() * ge2, unitLoTri * ge2); - VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView(), ge2 * unitLoTri); - VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView(), ge2 * unitLoTri); - VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView() * ge2.adjoint(), internal::conj(s1) * unitLoTri.adjoint() * ge2.adjoint()); - - VERIFY_IS_APPROX( ge3 = triV.template triangularView() * ge2, strictlyLoTri * ge2); - VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView(), ge2 * strictlyLoTri); - VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView(), ge2 * strictlyLoTri); - VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView() * ge2.adjoint(), internal::conj(s1) * strictlyLoTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX( ge_xs = mat.template triangularView() * ge_right, tri * ge_right); + VERIFY_IS_APPROX( ge_sx = ge_left * mat.template triangularView(), ge_left * tri); + + VERIFY_IS_APPROX( ge_xs.noalias() = mat.template triangularView() * ge_right, tri * ge_right); + VERIFY_IS_APPROX( ge_sx.noalias() = ge_left * mat.template triangularView(), ge_left * tri); + + VERIFY_IS_APPROX( ge_xs.noalias() = (s1*mat.adjoint()).template triangularView() * (s2*ge_left.transpose()), s1*triTr.conjugate() * (s2*ge_left.transpose())); + VERIFY_IS_APPROX( ge_sx.noalias() = ge_right.transpose() * mat.adjoint().template triangularView(), ge_right.transpose() * triTr.conjugate()); + + VERIFY_IS_APPROX( ge_xs.noalias() = (s1*mat.adjoint()).template triangularView() * (s2*ge_left.adjoint()), s1*triTr.conjugate() * (s2*ge_left.adjoint())); + VERIFY_IS_APPROX( ge_sx.noalias() = ge_right.adjoint() * mat.adjoint().template triangularView(), ge_right.adjoint() * triTr.conjugate()); + + VERIFY_IS_APPROX( (ge_xs+s1*triTr.conjugate() * (s2*ge_left.adjoint())).eval(), ge_xs.noalias() += (s1*mat.adjoint()).template triangularView() * (s2*ge_left.adjoint()) ); + VERIFY_IS_APPROX( (ge_sx-ge_right.adjoint() * triTr.conjugate()).eval(), ge_sx.noalias() -= ge_right.adjoint() * mat.adjoint().template triangularView()); + + VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView() * ge_left.adjoint(), internal::conj(s1) * triTr.conjugate() * ge_left.adjoint()); + + // TODO check with sub-matrix expressions ? } +template +void trmv(int rows=internal::random(1,320), int cols=internal::random(1,320)) +{ + trmm(rows,cols,1); +} + +template +void trmm(int rows=internal::random(1,320), int cols=internal::random(1,320), int otherCols = internal::random(1,320)) +{ + trmm(rows,cols,otherCols); +} + +#define CALL_ALL_ORDERS(NB,SCALAR,MODE) \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + EIGEN_CAT(CALL_SUBTEST_,NB)((trmm())); \ + \ + EIGEN_CAT(CALL_SUBTEST_1,NB)((trmv())); \ + EIGEN_CAT(CALL_SUBTEST_1,NB)((trmv())); + + +#define CALL_ALL(NB,SCALAR) \ + CALL_ALL_ORDERS(EIGEN_CAT(1,NB),SCALAR,Upper) \ + CALL_ALL_ORDERS(EIGEN_CAT(2,NB),SCALAR,UnitUpper) \ + CALL_ALL_ORDERS(EIGEN_CAT(3,NB),SCALAR,StrictlyUpper) \ + CALL_ALL_ORDERS(EIGEN_CAT(1,NB),SCALAR,Lower) \ + CALL_ALL_ORDERS(EIGEN_CAT(2,NB),SCALAR,UnitLower) \ + CALL_ALL_ORDERS(EIGEN_CAT(3,NB),SCALAR,StrictlyLower) + + void test_product_trmm() { for(int i = 0; i < g_repeat ; i++) { - CALL_SUBTEST_1((trmm(internal::random(1,320),internal::random(1,320)))); - CALL_SUBTEST_2((trmm(internal::random(1,320),internal::random(1,320)))); - CALL_SUBTEST_3((trmm >(internal::random(1,200),internal::random(1,200)))); - CALL_SUBTEST_4((trmm >(internal::random(1,200),internal::random(1,200)))); + CALL_ALL(1,float); // EIGEN_SUFFIXES;11;111;21;121;31;131 + CALL_ALL(2,double); // EIGEN_SUFFIXES;12;112;22;122;32;132 + CALL_ALL(3,std::complex); // EIGEN_SUFFIXES;13;113;23;123;33;133 + CALL_ALL(4,std::complex); // EIGEN_SUFFIXES;14;114;24;124;34;134 } }