| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #include "main.h" |
| |
| template <typename T> |
| int get_random_size() { |
| const int factor = NumTraits<T>::ReadCost; |
| const int max_test_size = EIGEN_TEST_MAX_SIZE > 2 * factor ? EIGEN_TEST_MAX_SIZE / factor : EIGEN_TEST_MAX_SIZE; |
| return internal::random<int>(1, max_test_size); |
| } |
| |
| template <typename Scalar, int Mode, int TriOrder, int OtherOrder, int ResOrder, int OtherCols> |
| void trmm(int rows = get_random_size<Scalar>(), int cols = get_random_size<Scalar>(), |
| int otherCols = OtherCols == Dynamic ? get_random_size<Scalar>() : OtherCols) { |
| typedef Matrix<Scalar, Dynamic, Dynamic, TriOrder> TriMatrix; |
| typedef Matrix<Scalar, Dynamic, OtherCols, OtherCols == 1 ? ColMajor : OtherOrder> OnTheRight; |
| typedef Matrix<Scalar, OtherCols, Dynamic, OtherCols == 1 ? RowMajor : OtherOrder> OnTheLeft; |
| |
| typedef Matrix<Scalar, Dynamic, OtherCols, OtherCols == 1 ? ColMajor : ResOrder> ResXS; |
| typedef Matrix<Scalar, OtherCols, Dynamic, OtherCols == 1 ? RowMajor : ResOrder> ResSX; |
| |
| TriMatrix mat(rows, cols), tri(rows, cols), triTr(cols, rows), s1tri(rows, cols), s1triTr(cols, rows); |
| |
| OnTheRight ge_right(cols, otherCols); |
| OnTheLeft ge_left(otherCols, rows); |
| ResSX ge_sx, ge_sx_save; |
| ResXS ge_xs, ge_xs_save; |
| |
| Scalar s1 = internal::random<Scalar>(), s2 = internal::random<Scalar>(); |
| |
| mat.setRandom(); |
| tri = mat.template triangularView<Mode>(); |
| triTr = mat.transpose().template triangularView<Mode>(); |
| s1tri = (s1 * mat).template triangularView<Mode>(); |
| s1triTr = (s1 * mat).transpose().template triangularView<Mode>(); |
| ge_right.setRandom(); |
| ge_left.setRandom(); |
| |
| VERIFY_IS_APPROX(ge_xs = mat.template triangularView<Mode>() * ge_right, tri * ge_right); |
| VERIFY_IS_APPROX(ge_sx = ge_left * mat.template triangularView<Mode>(), ge_left * tri); |
| |
| VERIFY_IS_APPROX(ge_xs.noalias() = mat.template triangularView<Mode>() * ge_right, tri * ge_right); |
| VERIFY_IS_APPROX(ge_sx.noalias() = ge_left * mat.template triangularView<Mode>(), ge_left * tri); |
| |
| if ((Mode & UnitDiag) == 0) |
| VERIFY_IS_APPROX( |
| ge_xs.noalias() = (s1 * mat.adjoint()).template triangularView<Mode>() * (s2 * ge_left.transpose()), |
| s1 * triTr.conjugate() * (s2 * ge_left.transpose())); |
| |
| VERIFY_IS_APPROX( |
| ge_xs.noalias() = (s1 * mat.transpose()).template triangularView<Mode>() * (s2 * ge_left.transpose()), |
| s1triTr * (s2 * ge_left.transpose())); |
| VERIFY_IS_APPROX(ge_sx.noalias() = (s2 * ge_left) * (s1 * mat).template triangularView<Mode>(), |
| (s2 * ge_left) * s1tri); |
| |
| VERIFY_IS_APPROX(ge_sx.noalias() = ge_right.transpose() * mat.adjoint().template triangularView<Mode>(), |
| ge_right.transpose() * triTr.conjugate()); |
| VERIFY_IS_APPROX(ge_sx.noalias() = ge_right.adjoint() * mat.adjoint().template triangularView<Mode>(), |
| ge_right.adjoint() * triTr.conjugate()); |
| |
| ge_xs_save = ge_xs; |
| if ((Mode & UnitDiag) == 0) |
| VERIFY_IS_APPROX( |
| (ge_xs_save + s1 * triTr.conjugate() * (s2 * ge_left.adjoint())).eval(), |
| ge_xs.noalias() += (s1 * mat.adjoint()).template triangularView<Mode>() * (s2 * ge_left.adjoint())); |
| ge_xs_save = ge_xs; |
| VERIFY_IS_APPROX( |
| (ge_xs_save + s1triTr * (s2 * ge_left.adjoint())).eval(), |
| ge_xs.noalias() += (s1 * mat.transpose()).template triangularView<Mode>() * (s2 * ge_left.adjoint())); |
| ge_sx.setRandom(); |
| ge_sx_save = ge_sx; |
| if ((Mode & UnitDiag) == 0) |
| VERIFY_IS_APPROX( |
| ge_sx_save - (ge_right.adjoint() * (-s1 * triTr).conjugate()).eval(), |
| ge_sx.noalias() -= (ge_right.adjoint() * (-s1 * mat).adjoint().template triangularView<Mode>()).eval()); |
| |
| if ((Mode & UnitDiag) == 0) |
| VERIFY_IS_APPROX(ge_xs = (s1 * mat).adjoint().template triangularView<Mode>() * ge_left.adjoint(), |
| numext::conj(s1) * triTr.conjugate() * ge_left.adjoint()); |
| VERIFY_IS_APPROX(ge_xs = (s1 * mat).transpose().template triangularView<Mode>() * ge_left.adjoint(), |
| s1triTr * ge_left.adjoint()); |
| |
| // TODO check with sub-matrix expressions ? |
| |
| // destination with a non-default inner-stride |
| // see bug 1741 |
| { |
| VERIFY_IS_APPROX(ge_xs.noalias() = mat.template triangularView<Mode>() * ge_right, tri * ge_right); |
| typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX; |
| MatrixX buffer(2 * ge_xs.rows(), 2 * ge_xs.cols()); |
| Map<ResXS, 0, Stride<Dynamic, 2> > map1(buffer.data(), ge_xs.rows(), ge_xs.cols(), |
| Stride<Dynamic, 2>(2 * ge_xs.outerStride(), 2)); |
| buffer.setZero(); |
| VERIFY_IS_APPROX(map1.noalias() = mat.template triangularView<Mode>() * ge_right, tri * ge_right); |
| } |
| } |
| |
| template <typename Scalar, int Mode, int TriOrder> |
| void trmv(int rows = get_random_size<Scalar>(), int cols = get_random_size<Scalar>()) { |
| trmm<Scalar, Mode, TriOrder, ColMajor, ColMajor, 1>(rows, cols, 1); |
| } |
| |
| template <typename Scalar, int Mode, int TriOrder, int OtherOrder, int ResOrder> |
| void trmm(int rows = get_random_size<Scalar>(), int cols = get_random_size<Scalar>(), |
| int otherCols = get_random_size<Scalar>()) { |
| trmm<Scalar, Mode, TriOrder, OtherOrder, ResOrder, Dynamic>(rows, cols, otherCols); |
| } |
| |
| #define CALL_ALL_ORDERS(NB, SCALAR, MODE) \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, ColMajor, ColMajor, ColMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, ColMajor, ColMajor, RowMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, ColMajor, RowMajor, ColMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, ColMajor, RowMajor, RowMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, RowMajor, ColMajor, ColMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, RowMajor, ColMajor, RowMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, RowMajor, RowMajor, ColMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_, NB)((trmm<SCALAR, MODE, RowMajor, RowMajor, RowMajor>())); \ |
| \ |
| EIGEN_CAT(CALL_SUBTEST_1, NB)((trmv<SCALAR, MODE, ColMajor>())); \ |
| EIGEN_CAT(CALL_SUBTEST_1, NB)((trmv<SCALAR, MODE, RowMajor>())); |
| |
| #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) |
| |
| EIGEN_DECLARE_TEST(product_trmm) { |
| for (int i = 0; i < g_repeat; i++) { |
| 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<float>); // EIGEN_SUFFIXES;13;113;23;123;33;133 |
| CALL_ALL(4, std::complex<double>); // EIGEN_SUFFIXES;14;114;24;124;34;134 |
| } |
| } |
