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
operator ()coeff_wise51 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
operator ()replicate66 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
operator ()redux83 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
operator ()prod_test103 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
operator ()diagonal117 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
operator ()eigenvalues129 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
test_cuda_basic()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