ground/gcs/src/libs/eigen/test/cuda_basic.cu

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2015-2016 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 // workaround issue between gcc >= 4.7 and cuda 5.5
  11 #if (defined __GNUC__) && (__GNUC__>4 || __GNUC_MINOR__>=7)
  12   #undef _GLIBCXX_ATOMIC_BUILTINS
  13   #undef _GLIBCXX_USE_INT128
  14 #endif
  15
  16 #define EIGEN_TEST_NO_LONGDOUBLE
  17 #define EIGEN_TEST_NO_COMPLEX
  18 #define EIGEN_TEST_FUNC cuda_basic
  19 #define EIGEN_DEFAULT_DENSE_INDEX_TYPE int
  20
  21 #include <math_constants.h>
  22 #include <cuda.h>
  23 #if defined __CUDACC_VER__ && __CUDACC_VER__ >= 70500
  24 #include <cuda_fp16.h>
  25 #endif
  26 #include "main.h"
  27 #include "cuda_common.h"
  28
  29 // Check that dense modules can be properly parsed by nvcc
  30 #include <Eigen/Dense>
  31
  32 // struct Foo{
  33 //   EIGEN_DEVICE_FUNC
  34 //   void operator()(int i, const float* mats, float* vecs) const {
  35 //     using namespace Eigen;
  36 //   //   Matrix3f M(data);
  37 //   //   Vector3f x(data+9);
  38 //   //   Map<Vector3f>(data+9) = M.inverse() * x;
  39 //     Matrix3f M(mats+i/16);
  40 //     Vector3f x(vecs+i*3);
  41 //   //   using std::min;
  42 //   //   using std::sqrt;
  43 //     Map<Vector3f>(vecs+i*3) << x.minCoeff(), 1, 2;// / x.dot(x);//(M.inverse() *  x) / x.x();
  44 //     //x = x*2 + x.y() * x + x * x.maxCoeff() - x / x.sum();
  45 //   }
  46 // };
  47
  48 template<typename T>
  49 struct coeff_wise {
  50   EIGEN_DEVICE_FUNC
  51   void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const
  52   {
  53     using namespace Eigen;
  54     T x1(in+i);
  55     T x2(in+i+1);
  56     T x3(in+i+2);
  57     Map<T> res(out+i*T::MaxSizeAtCompileTime);
  58
  59     res.array() += (in[0] * x1 + x2).array() * x3.array();
  60   }
  61 };
  62
  63 template<typename T>
  64 struct replicate {
  65   EIGEN_DEVICE_FUNC
  66   void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const
  67   {
  68     using namespace Eigen;
  69     T x1(in+i);
  70     int step   = x1.size() * 4;
  71     int stride = 3 * step;
  72
  73     typedef Map<Array<typename T::Scalar,Dynamic,Dynamic> > MapType;
  74     MapType(out+i*stride+0*step, x1.rows()*2, x1.cols()*2) = x1.replicate(2,2);
  75     MapType(out+i*stride+1*step, x1.rows()*3, x1.cols()) = in[i] * x1.colwise().replicate(3);
  76     MapType(out+i*stride+2*step, x1.rows(), x1.cols()*3) = in[i] * x1.rowwise().replicate(3);
  77   }
  78 };
  79
  80 template<typename T>
  81 struct redux {
  82   EIGEN_DEVICE_FUNC
  83   void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const
  84   {
  85     using namespace Eigen;
  86     int N = 10;
  87     T x1(in+i);
  88     out[i*N+0] = x1.minCoeff();
  89     out[i*N+1] = x1.maxCoeff();
  90     out[i*N+2] = x1.sum();
  91     out[i*N+3] = x1.prod();
  92     out[i*N+4] = x1.matrix().squaredNorm();
  93     out[i*N+5] = x1.matrix().norm();
  94     out[i*N+6] = x1.colwise().sum().maxCoeff();
  95     out[i*N+7] = x1.rowwise().maxCoeff().sum();
  96     out[i*N+8] = x1.matrix().colwise().squaredNorm().sum();
  97   }
  98 };
  99
 100 template<typename T1, typename T2>
 101 struct prod_test {
 102   EIGEN_DEVICE_FUNC
 103   void operator()(int i, const typename T1::Scalar* in, typename T1::Scalar* out) const
 104   {
 105     using namespace Eigen;
 106     typedef Matrix<typename T1::Scalar, T1::RowsAtCompileTime, T2::ColsAtCompileTime> T3;
 107     T1 x1(in+i);
 108     T2 x2(in+i+1);
 109     Map<T3> res(out+i*T3::MaxSizeAtCompileTime);
 110     res += in[i] * x1 * x2;
 111   }
 112 };
 113
 114 template<typename T1, typename T2>
 115 struct diagonal {
 116   EIGEN_DEVICE_FUNC
 117   void operator()(int i, const typename T1::Scalar* in, typename T1::Scalar* out) const
 118   {
 119     using namespace Eigen;
 120     T1 x1(in+i);
 121     Map<T2> res(out+i*T2::MaxSizeAtCompileTime);
 122     res += x1.diagonal();
 123   }
 124 };
 125
 126 template<typename T>
 127 struct eigenvalues {
 128   EIGEN_DEVICE_FUNC
 129   void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const
 130   {
 131     using namespace Eigen;
 132     typedef Matrix<typename T::Scalar, T::RowsAtCompileTime, 1> Vec;
 133     T M(in+i);
 134     Map<Vec> res(out+i*Vec::MaxSizeAtCompileTime);
 135     T A = M*M.adjoint();
 136     SelfAdjointEigenSolver<T> eig;
 137     eig.computeDirect(M);
 138     res = eig.eigenvalues();
 139   }
 140 };
 141
 142 void test_cuda_basic()
 143 {
 144   ei_test_init_cuda();
 145
 146   int nthreads = 100;
 147   Eigen::VectorXf in, out;
 148
 149   #ifndef __CUDA_ARCH__
 150   int data_size = nthreads * 512;
 151   in.setRandom(data_size);
 152   out.setRandom(data_size);
 153   #endif
 154
 155   CALL_SUBTEST( run_and_compare_to_cuda(coeff_wise<Vector3f>(), nthreads, in, out) );
 156   CALL_SUBTEST( run_and_compare_to_cuda(coeff_wise<Array44f>(), nthreads, in, out) );
 157
 158   CALL_SUBTEST( run_and_compare_to_cuda(replicate<Array4f>(), nthreads, in, out) );
 159   CALL_SUBTEST( run_and_compare_to_cuda(replicate<Array33f>(), nthreads, in, out) );
 160
 161   CALL_SUBTEST( run_and_compare_to_cuda(redux<Array4f>(), nthreads, in, out) );
 162   CALL_SUBTEST( run_and_compare_to_cuda(redux<Matrix3f>(), nthreads, in, out) );
 163
 164   CALL_SUBTEST( run_and_compare_to_cuda(prod_test<Matrix3f,Matrix3f>(), nthreads, in, out) );
 165   CALL_SUBTEST( run_and_compare_to_cuda(prod_test<Matrix4f,Vector4f>(), nthreads, in, out) );
 166
 167   CALL_SUBTEST( run_and_compare_to_cuda(diagonal<Matrix3f,Vector3f>(), nthreads, in, out) );
 168   CALL_SUBTEST( run_and_compare_to_cuda(diagonal<Matrix4f,Vector4f>(), nthreads, in, out) );
 169
 170   CALL_SUBTEST( run_and_compare_to_cuda(eigenvalues<Matrix3f>(), nthreads, in, out) );
 171   CALL_SUBTEST( run_and_compare_to_cuda(eigenvalues<Matrix2f>(), nthreads, in, out) );
 172
 173 }