Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
| 10 | #include "main.h" |
| 11 | #include <Eigen/QR> |
| 12 | |
| 13 | template<typename Derived1, typename Derived2> |
| 14 | bool areNotApprox(const MatrixBase<Derived1>& m1, const MatrixBase<Derived2>& m2, typename Derived1::RealScalar epsilon = NumTraits<typename Derived1::RealScalar>::dummy_precision()) |
| 15 | { |
| 16 | return !((m1-m2).cwiseAbs2().maxCoeff() < epsilon * epsilon |
| 17 | * (std::max)(m1.cwiseAbs2().maxCoeff(), m2.cwiseAbs2().maxCoeff())); |
| 18 | } |
| 19 | |
| 20 | template<typename MatrixType> void product(const MatrixType& m) |
| 21 | { |
| 22 | /* this test covers the following files: |
| 23 | Identity.h Product.h |
| 24 | */ |
| 25 | typedef typename MatrixType::Index Index; |
| 26 | typedef typename MatrixType::Scalar Scalar; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 27 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RowVectorType; |
| 28 | typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ColVectorType; |
| 29 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType; |
| 30 | typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> ColSquareMatrixType; |
| 31 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, |
| 32 | MatrixType::Flags&RowMajorBit?ColMajor:RowMajor> OtherMajorMatrixType; |
| 33 | |
| 34 | Index rows = m.rows(); |
| 35 | Index cols = m.cols(); |
| 36 | |
| 37 | // this test relies a lot on Random.h, and there's not much more that we can do |
| 38 | // to test it, hence I consider that we will have tested Random.h |
| 39 | MatrixType m1 = MatrixType::Random(rows, cols), |
| 40 | m2 = MatrixType::Random(rows, cols), |
| 41 | m3(rows, cols); |
| 42 | RowSquareMatrixType |
| 43 | identity = RowSquareMatrixType::Identity(rows, rows), |
| 44 | square = RowSquareMatrixType::Random(rows, rows), |
| 45 | res = RowSquareMatrixType::Random(rows, rows); |
| 46 | ColSquareMatrixType |
| 47 | square2 = ColSquareMatrixType::Random(cols, cols), |
| 48 | res2 = ColSquareMatrixType::Random(cols, cols); |
| 49 | RowVectorType v1 = RowVectorType::Random(rows); |
| 50 | ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols); |
| 51 | OtherMajorMatrixType tm1 = m1; |
| 52 | |
| 53 | Scalar s1 = internal::random<Scalar>(); |
| 54 | |
| 55 | Index r = internal::random<Index>(0, rows-1), |
| 56 | c = internal::random<Index>(0, cols-1), |
| 57 | c2 = internal::random<Index>(0, cols-1); |
| 58 | |
| 59 | // begin testing Product.h: only associativity for now |
| 60 | // (we use Transpose.h but this doesn't count as a test for it) |
| 61 | VERIFY_IS_APPROX((m1*m1.transpose())*m2, m1*(m1.transpose()*m2)); |
| 62 | m3 = m1; |
| 63 | m3 *= m1.transpose() * m2; |
| 64 | VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); |
| 65 | VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); |
| 66 | |
| 67 | // continue testing Product.h: distributivity |
| 68 | VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2); |
| 69 | VERIFY_IS_APPROX(square*(m1 - m2), square*m1-square*m2); |
| 70 | |
| 71 | // continue testing Product.h: compatibility with ScalarMultiple.h |
| 72 | VERIFY_IS_APPROX(s1*(square*m1), (s1*square)*m1); |
| 73 | VERIFY_IS_APPROX(s1*(square*m1), square*(m1*s1)); |
| 74 | |
| 75 | // test Product.h together with Identity.h |
| 76 | VERIFY_IS_APPROX(v1, identity*v1); |
| 77 | VERIFY_IS_APPROX(v1.transpose(), v1.transpose() * identity); |
| 78 | // again, test operator() to check const-qualification |
| 79 | VERIFY_IS_APPROX(MatrixType::Identity(rows, cols)(r,c), static_cast<Scalar>(r==c)); |
| 80 | |
| 81 | if (rows!=cols) |
| 82 | VERIFY_RAISES_ASSERT(m3 = m1*m1); |
| 83 | |
| 84 | // test the previous tests were not screwed up because operator* returns 0 |
| 85 | // (we use the more accurate default epsilon) |
| 86 | if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1) |
| 87 | { |
| 88 | VERIFY(areNotApprox(m1.transpose()*m2,m2.transpose()*m1)); |
| 89 | } |
| 90 | |
| 91 | // test optimized operator+= path |
| 92 | res = square; |
| 93 | res.noalias() += m1 * m2.transpose(); |
| 94 | VERIFY_IS_APPROX(res, square + m1 * m2.transpose()); |
| 95 | if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1) |
| 96 | { |
| 97 | VERIFY(areNotApprox(res,square + m2 * m1.transpose())); |
| 98 | } |
| 99 | vcres = vc2; |
| 100 | vcres.noalias() += m1.transpose() * v1; |
| 101 | VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1); |
| 102 | |
| 103 | // test optimized operator-= path |
| 104 | res = square; |
| 105 | res.noalias() -= m1 * m2.transpose(); |
| 106 | VERIFY_IS_APPROX(res, square - (m1 * m2.transpose())); |
| 107 | if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1) |
| 108 | { |
| 109 | VERIFY(areNotApprox(res,square - m2 * m1.transpose())); |
| 110 | } |
| 111 | vcres = vc2; |
| 112 | vcres.noalias() -= m1.transpose() * v1; |
| 113 | VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1); |
| 114 | |
| 115 | tm1 = m1; |
| 116 | VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1); |
| 117 | VERIFY_IS_APPROX(v1.transpose() * tm1, v1.transpose() * m1); |
| 118 | |
| 119 | // test submatrix and matrix/vector product |
| 120 | for (int i=0; i<rows; ++i) |
| 121 | res.row(i) = m1.row(i) * m2.transpose(); |
| 122 | VERIFY_IS_APPROX(res, m1 * m2.transpose()); |
| 123 | // the other way round: |
| 124 | for (int i=0; i<rows; ++i) |
| 125 | res.col(i) = m1 * m2.transpose().col(i); |
| 126 | VERIFY_IS_APPROX(res, m1 * m2.transpose()); |
| 127 | |
| 128 | res2 = square2; |
| 129 | res2.noalias() += m1.transpose() * m2; |
| 130 | VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2); |
| 131 | if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1) |
| 132 | { |
| 133 | VERIFY(areNotApprox(res2,square2 + m2.transpose() * m1)); |
| 134 | } |
| 135 | |
| 136 | VERIFY_IS_APPROX(res.col(r).noalias() = square.adjoint() * square.col(r), (square.adjoint() * square.col(r)).eval()); |
| 137 | VERIFY_IS_APPROX(res.col(r).noalias() = square * square.col(r), (square * square.col(r)).eval()); |
| 138 | |
| 139 | // inner product |
| 140 | Scalar x = square2.row(c) * square2.col(c2); |
| 141 | VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum()); |
| 142 | } |