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 (a] 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 } 174
