Warning, /DataFormats/Math/test/simpleCholeskyTest.cu is written in an unsupported language. File is not indexed.
0001 #include "DataFormats/Math/interface/choleskyInversion.h"
0002
0003 using namespace math::cholesky;
0004
0005 #include <Eigen/Core>
0006 #include <Eigen/Eigenvalues>
0007
0008 #include <iomanip>
0009 #include <memory>
0010 #include <algorithm>
0011 #include <chrono>
0012 #include <random>
0013
0014 #include <cassert>
0015 #include <iostream>
0016 #include <limits>
0017
0018 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
0019 #include "HeterogeneousCore/CUDAUtilities/interface/requireDevices.h"
0020
0021 template <typename M, int N>
0022 __global__ void invert(M* mm, int n) {
0023 auto i = blockIdx.x * blockDim.x + threadIdx.x;
0024 if (i >= n)
0025 return;
0026
0027 auto& m = mm[i];
0028
0029 printf("before %d %f %f %f\n", N, m(0, 0), m(1, 0), m(1, 1));
0030
0031 invertNN(m, m);
0032
0033 printf("after %d %f %f %f\n", N, m(0, 0), m(1, 0), m(1, 1));
0034 }
0035
0036 template <typename M, int N>
0037 __global__ void invertE(M* mm, int n) {
0038 auto i = blockIdx.x * blockDim.x + threadIdx.x;
0039 if (i >= n)
0040 return;
0041
0042 auto& m = mm[i];
0043
0044 printf("before %d %f %f %f\n", N, m(0, 0), m(1, 0), m(1, 1));
0045
0046 m = m.inverse();
0047
0048 printf("after %d %f %f %f\n", N, m(0, 0), m(1, 0), m(1, 1));
0049 }
0050
0051 // generate matrices
0052 template <class M>
0053 void genMatrix(M& m) {
0054 using T = typename std::remove_reference<decltype(m(0, 0))>::type;
0055 int n = M::ColsAtCompileTime;
0056 std::mt19937 eng;
0057 // std::mt19937 eng2;
0058 std::uniform_real_distribution<T> rgen(0., 1.);
0059
0060 // generate first diagonal elemets
0061 for (int i = 0; i < n; ++i) {
0062 double maxVal = i * 10000 / (n - 1) + 1; // max condition is 10^4
0063 m(i, i) = maxVal * rgen(eng);
0064 }
0065 for (int i = 0; i < n; ++i) {
0066 for (int j = 0; j < i; ++j) {
0067 double v = 0.3 * std::sqrt(m(i, i) * m(j, j)); // this makes the matrix pos defined
0068 m(i, j) = v * rgen(eng);
0069 m(j, i) = m(i, j);
0070 }
0071 }
0072 }
0073
0074 template <int DIM>
0075 using MXN = Eigen::Matrix<double, DIM, DIM>;
0076
0077 int main() {
0078 cms::cudatest::requireDevices();
0079
0080 constexpr int DIM = 6;
0081
0082 using M = MXN<DIM>;
0083
0084 M m;
0085
0086 genMatrix(m);
0087
0088 printf("on CPU before %d %f %f %f\n", DIM, m(0, 0), m(1, 0), m(1, 1));
0089
0090 invertNN(m, m);
0091
0092 printf("on CPU after %d %f %f %f\n", DIM, m(0, 0), m(1, 0), m(1, 1));
0093
0094 double* d;
0095 cudaMalloc(&d, sizeof(M));
0096 cudaMemcpy(d, &m, sizeof(M), cudaMemcpyHostToDevice);
0097 invert<M, DIM><<<1, 1>>>((M*)d, 1);
0098 cudaCheck(cudaDeviceSynchronize());
0099 invertE<M, DIM><<<1, 1>>>((M*)d, 1);
0100 cudaCheck(cudaDeviceSynchronize());
0101
0102 return 0;
0103 }
