ground/gcs/src/libs/eigen/test/product_syrk.cpp

   1 // This file is part 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
  10 #include "main.h"
  11
  12 template<typename MatrixType> void syrk(const MatrixType& m)
  13 {
  14   typedef typename MatrixType::Index Index;
  15   typedef typename MatrixType::Scalar Scalar;
  16   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, RowMajor> RMatrixType;
  17   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1;
  18   typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2;
  19   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> Rhs3;
  20
  21   Index rows = m.rows();
  22   Index cols = m.cols();
  23
  24   MatrixType m1 = MatrixType::Random(rows, cols),
  25              m2 = MatrixType::Random(rows, cols),
  26              m3 = MatrixType::Random(rows, cols);
  27   RMatrixType rm2 = MatrixType::Random(rows, cols);
  28
  29   Rhs1 rhs1 = Rhs1::Random(internal::random<int>(1,320), cols); Rhs1 rhs11 = Rhs1::Random(rhs1.rows(), cols);
  30   Rhs2 rhs2 = Rhs2::Random(rows, internal::random<int>(1,320)); Rhs2 rhs22 = Rhs2::Random(rows, rhs2.cols());
  31   Rhs3 rhs3 = Rhs3::Random(internal::random<int>(1,320), rows);
  32
  33   Scalar s1 = internal::random<Scalar>();
  34
  35   Index c = internal::random<Index>(0,cols-1);
  36
  37   m2.setZero();
  38   VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()),
  39                    ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
  40   m2.setZero();
  41   VERIFY_IS_APPROX(((m2.template triangularView<Lower>() += s1 * rhs2  * rhs22.adjoint()).nestedExpression()),
  42                    ((s1 * rhs2 * rhs22.adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
  43
  44
  45   m2.setZero();
  46   VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(),
  47                    (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix());
  48   m2.setZero();
  49   VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * rhs22 * rhs2.adjoint()).nestedExpression(),
  50                    (s1 * rhs22 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix());
  51
  52
  53   m2.setZero();
  54   VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(),
  55                    (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix());
  56   m2.setZero();
  57   VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * rhs11.adjoint() * rhs1).nestedExpression(),
  58                    (s1 * rhs11.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix());
  59
  60
  61   m2.setZero();
  62   VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(),
  63                    (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix());
  64   VERIFY_IS_APPROX((m2.template triangularView<Upper>() = s1 * rhs1.adjoint() * rhs11).nestedExpression(),
  65                    (s1 * rhs1.adjoint() * rhs11).eval().template triangularView<Upper>().toDenseMatrix());
  66
  67
  68   m2.setZero();
  69   VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(),
  70                    (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix());
  71
  72   m2.setZero();
  73   VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs3.adjoint(),s1)._expression(),
  74                    (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Upper>().toDenseMatrix());
  75
  76   m2.setZero();
  77   VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c),s1)._expression()),
  78                    ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
  79
  80   m2.setZero();
  81   VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()),
  82                    ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
  83   rm2.setZero();
  84   VERIFY_IS_APPROX((rm2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()),
  85                    ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
  86   m2.setZero();
  87   VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * m3.col(c) * m1.col(c).adjoint()).nestedExpression(),
  88                    ((s1 * m3.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
  89   rm2.setZero();
  90   VERIFY_IS_APPROX((rm2.template triangularView<Upper>() += s1 * m1.col(c) * m3.col(c).adjoint()).nestedExpression(),
  91                    ((s1 * m1.col(c) * m3.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
  92
  93   m2.setZero();
  94   VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c).conjugate(),s1)._expression()),
  95                    ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
  96
  97   m2.setZero();
  98   VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c).conjugate(),s1)._expression()),
  99                    ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
 100
 101
 102   m2.setZero();
 103   VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()),
 104                    ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
 105   rm2.setZero();
 106   VERIFY_IS_APPROX((rm2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()),
 107                    ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
 108   m2.setZero();
 109   VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(),
 110                    ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
 111   rm2.setZero();
 112   VERIFY_IS_APPROX((rm2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(),
 113                    ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
 114
 115
 116   m2.setZero();
 117   VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.row(c).adjoint(),s1)._expression()),
 118                    ((s1 * m1.row(c).adjoint() * m1.row(c).adjoint().adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
 119 }
 120
 121 void test_product_syrk()
 122 {
 123   for(int i = 0; i < g_repeat ; i++)
 124   {
 125     int s;
 126     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
 127     CALL_SUBTEST_1( syrk(MatrixXf(s, s)) );
 128     CALL_SUBTEST_2( syrk(MatrixXd(s, s)) );
 129     TEST_SET_BUT_UNUSED_VARIABLE(s)
 130
 131     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
 132     CALL_SUBTEST_3( syrk(MatrixXcf(s, s)) );
 133     CALL_SUBTEST_4( syrk(MatrixXcd(s, s)) );
 134     TEST_SET_BUT_UNUSED_VARIABLE(s)
 135   }
 136 }