product_trmv.cpp 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990
  1. // This file is triangularView of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
  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. #include "main.h"
  10. template<typename MatrixType> void trmv(const MatrixType& m)
  11. {
  12. typedef typename MatrixType::Scalar Scalar;
  13. typedef typename NumTraits<Scalar>::Real RealScalar;
  14. typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  15. RealScalar largerEps = 10*test_precision<RealScalar>();
  16. Index rows = m.rows();
  17. Index cols = m.cols();
  18. MatrixType m1 = MatrixType::Random(rows, cols),
  19. m3(rows, cols);
  20. VectorType v1 = VectorType::Random(rows);
  21. Scalar s1 = internal::random<Scalar>();
  22. m1 = MatrixType::Random(rows, cols);
  23. // check with a column-major matrix
  24. m3 = m1.template triangularView<Eigen::Lower>();
  25. VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Lower>() * v1, largerEps));
  26. m3 = m1.template triangularView<Eigen::Upper>();
  27. VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Upper>() * v1, largerEps));
  28. m3 = m1.template triangularView<Eigen::UnitLower>();
  29. VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLower>() * v1, largerEps));
  30. m3 = m1.template triangularView<Eigen::UnitUpper>();
  31. VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpper>() * v1, largerEps));
  32. // check conjugated and scalar multiple expressions (col-major)
  33. m3 = m1.template triangularView<Eigen::Lower>();
  34. VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::Lower>() * v1, largerEps));
  35. m3 = m1.template triangularView<Eigen::Upper>();
  36. VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::Upper>() * v1.conjugate(), largerEps));
  37. // check with a row-major matrix
  38. m3 = m1.template triangularView<Eigen::Upper>();
  39. VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Lower>() * v1, largerEps));
  40. m3 = m1.template triangularView<Eigen::Lower>();
  41. VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Upper>() * v1, largerEps));
  42. m3 = m1.template triangularView<Eigen::UnitUpper>();
  43. VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLower>() * v1, largerEps));
  44. m3 = m1.template triangularView<Eigen::UnitLower>();
  45. VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpper>() * v1, largerEps));
  46. // check conjugated and scalar multiple expressions (row-major)
  47. m3 = m1.template triangularView<Eigen::Upper>();
  48. VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::Lower>() * v1, largerEps));
  49. m3 = m1.template triangularView<Eigen::Lower>();
  50. VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::Upper>() * (s1*v1.conjugate()), largerEps));
  51. m3 = m1.template triangularView<Eigen::UnitUpper>();
  52. // check transposed cases:
  53. m3 = m1.template triangularView<Eigen::Lower>();
  54. VERIFY((v1.transpose() * m3).isApprox(v1.transpose() * m1.template triangularView<Eigen::Lower>(), largerEps));
  55. VERIFY((v1.adjoint() * m3).isApprox(v1.adjoint() * m1.template triangularView<Eigen::Lower>(), largerEps));
  56. VERIFY((v1.adjoint() * m3.adjoint()).isApprox(v1.adjoint() * m1.template triangularView<Eigen::Lower>().adjoint(), largerEps));
  57. // TODO check with sub-matrices
  58. }
  59. EIGEN_DECLARE_TEST(product_trmv)
  60. {
  61. int s = 0;
  62. for(int i = 0; i < g_repeat ; i++) {
  63. CALL_SUBTEST_1( trmv(Matrix<float, 1, 1>()) );
  64. CALL_SUBTEST_2( trmv(Matrix<float, 2, 2>()) );
  65. CALL_SUBTEST_3( trmv(Matrix3d()) );
  66. s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
  67. CALL_SUBTEST_4( trmv(MatrixXcf(s,s)) );
  68. CALL_SUBTEST_5( trmv(MatrixXcd(s,s)) );
  69. TEST_SET_BUT_UNUSED_VARIABLE(s)
  70. s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
  71. CALL_SUBTEST_6( trmv(Matrix<float,Dynamic,Dynamic,RowMajor>(s, s)) );
  72. TEST_SET_BUT_UNUSED_VARIABLE(s)
  73. }
  74. }