File: batched_determinant.cpp

package info (click to toggle)
boost1.62 1.62.0%2Bdfsg-4
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 686,420 kB
  • sloc: cpp: 2,609,004; xml: 972,558; ansic: 53,674; python: 32,437; sh: 8,829; asm: 3,071; cs: 2,121; makefile: 964; perl: 859; yacc: 472; php: 132; ruby: 94; f90: 55; sql: 13; csh: 6
file content (96 lines) | stat: -rw-r--r-- 3,477 bytes parent folder | download | duplicates (14)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
//---------------------------------------------------------------------------//
// Copyright (c) 2013-2014 Kyle Lutz <kyle.r.lutz@gmail.com>
//
// Distributed under the Boost Software License, Version 1.0
// See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt
//
// See http://boostorg.github.com/compute for more information.
//---------------------------------------------------------------------------//

#include <iostream>

#include <Eigen/Core>
#include <Eigen/LU>

#include <boost/compute/function.hpp>
#include <boost/compute/system.hpp>
#include <boost/compute/algorithm/transform.hpp>
#include <boost/compute/container/vector.hpp>
#include <boost/compute/types/fundamental.hpp>

namespace compute = boost::compute;

// this example shows how to compute the determinant of many 4x4 matrices
// using a determinant function and the transform() algorithm. in OpenCL the
// float16 type can be used to store a 4x4 matrix and the components are laid
// out in the following order:
//
// M = [ s0 s4 s8 sc ]
//     [ s1 s5 s9 sd ]
//     [ s2 s6 sa se ]
//     [ s3 s7 sb sf ]
//
// the input matrices are created using eigen's random matrix and then
// used again at the end to verify the results of the determinant function.
int main()
{
    // get default device and setup context
    compute::device gpu = compute::system::default_device();
    compute::context context(gpu);
    compute::command_queue queue(context, gpu);
    std::cout << "device: " << gpu.name() << std::endl;

    size_t n = 1000;

    // create random 4x4 matrices on the host
    std::vector<Eigen::Matrix4f> matrices(n);
    for(size_t i = 0; i < n; i++){
        matrices[i] = Eigen::Matrix4f::Random();
    }

    // copy matrices to the device
    using compute::float16_;
    compute::vector<float16_> input(n, context);
    compute::copy(
        matrices.begin(), matrices.end(), input.begin(), queue
    );

    // function returning the determinant of a 4x4 matrix.
    BOOST_COMPUTE_FUNCTION(float, determinant4x4, (const float16_ m),
    {
        return m.s0*m.s5*m.sa*m.sf + m.s0*m.s6*m.sb*m.sd + m.s0*m.s7*m.s9*m.se +
               m.s1*m.s4*m.sb*m.se + m.s1*m.s6*m.s8*m.sf + m.s1*m.s7*m.sa*m.sc +
               m.s2*m.s4*m.s9*m.sf + m.s2*m.s5*m.sb*m.sc + m.s2*m.s7*m.s8*m.sd +
               m.s3*m.s4*m.sa*m.sd + m.s3*m.s5*m.s8*m.se + m.s3*m.s6*m.s9*m.sc -
               m.s0*m.s5*m.sb*m.se - m.s0*m.s6*m.s9*m.sf - m.s0*m.s7*m.sa*m.sd -
               m.s1*m.s4*m.sa*m.sf - m.s1*m.s6*m.sb*m.sc - m.s1*m.s7*m.s8*m.se -
               m.s2*m.s4*m.sb*m.sd - m.s2*m.s5*m.s8*m.sf - m.s2*m.s7*m.s9*m.sc -
               m.s3*m.s4*m.s9*m.se - m.s3*m.s5*m.sa*m.sc - m.s3*m.s6*m.s8*m.sd;
    });

    // calculate determinants on the gpu
    compute::vector<float> determinants(n, context);
    compute::transform(
        input.begin(), input.end(), determinants.begin(), determinant4x4, queue
    );

    // check determinants
    std::vector<float> host_determinants(n);
    compute::copy(
        determinants.begin(), determinants.end(), host_determinants.begin(), queue
    );

    for(size_t i = 0; i < n; i++){
        float det = matrices[i].determinant();

        if(std::abs(det - host_determinants[i]) > 1e-6){
          std::cerr << "error: wrong determinant at " << i << " ("
                    << host_determinants[i] << " != " << det << ")"
                    << std::endl;
            return -1;
        }
    }

    return 0;
}