Search code examples
cudagputhrust

calling Thrust device_vector from a device function


I have a struct Cap which inside I have a thrust::device_vector of another structure. When I compile the code, I get an error which complains about calling a host function (thrust::device_vector<FloatIntPair>) from a device function SphericalFaceManager::makeCaps. When I add __host__ __device__ instead of only __device__ to the member functions and constructors the code then compiles but I receive a warning same as aforementioned error and I think it copies data between host and device. My question is how can I access to device vectors in my classes avoiding any data transfer between CPU and GPU?

Hereafter you can find the code:

struct ParticleID {

Int solver;
Int ngb;
Int oldNgb;
LLInt no;
LLInt masterNo;

__device__ ParticleID() {
    solver = -8;
    ngb = 0;
    oldNgb = 0;
    no = 0;
    masterNo = -1;
}
};


struct BaseParticle {

Float h;
Float3 pos;
ParticleID id;

__device__ BaseParticle(const Float3& _pos, const Float& _h, const ParticleID& _id) :
     h(_h), pos(_pos), id(_id) { }

};


struct FloatIntPair{

Float first;
Int second;

__device__ FloatIntPair(const Float& _first, Int _second) : first(_first), second(_second) { }
__device__ FloatIntPair(const FloatIntPair& sample) : first(sample.first), second(sample.second) { }

static struct {
    __device__ bool operator()(const FloatIntPair& a, const FloatIntPair& b) {  return a.first < b.first; }
} LessOp;
};


struct Cap {

Float3 eX;
Float3 eY;
Float radius;
Float height;

Float3 center;
Float3 normal;

BaseParticle* aP;
BaseParticle* bP;

thrust::device_vector<FloatIntPair> vertices; // The ordered list of vertices generated from intersections by other circles

__device__ inline Float findAngle(const Float3& vertex) const {

    Float result;
    Float3 r = (vertex - center);
    result = atan2(r|eY,r|eX);
    return result += (result < 0.0) * (2.0 * _PI);
}

__device__ void insertVertex(const Float3& vertex, Int id) {

    Float theta;
    if (!vertices.empty())
        theta = findAngle(vertex);
    else {
        eX = normalVec(vertex - center);
        eY = normal ^ eX;
        theta = 0.0;
    }
    vertices.push_back(FloatIntPair(theta,id));
}

__device__ Cap(BaseParticle* _aP, BaseParticle* _bP) : aP(_aP), bP(_bP) {

    //Compute normal, center, radius
    Float d = mag(bP->pos - aP->pos);
    if(d == 0.0){
        normal = Vector1(0.0);
        center = aP->pos;
        radius = height = 0.0;
    } else {
        normal = (bP->pos - aP->pos) / d;
        Float x = (d * d - bP->h * bP->h + aP->h * aP->h) / (2.0 * d);
        center = aP->pos + normal * x;
        if (x >= aP->h) {
            radius = height = 0.0;
            return;
        }
        radius = sqrt(aP->h * aP->h - x * x);
        height = min(2.0 * aP->h, aP->h - x);

        Float3 vec001 = Vector(0.0,0.0,1.0);
            Float3 vec011 = Vector(0.0,1.0,1.0);

        eX = normalVec(vec001 ^ normal);
        if (mag2(eX) < geoEps()) {
            eX = eX = normalVec(vec011 ^ normal);
        }

        eY = normal ^ eX;
    }
}
};

class SphericalFaceManager {
BaseParticle* particle;
Int baseSigma;
public:
thrust::device_vector<Cap> caps;  
thrust::device_vector<Float3> vertexPool;    
__device__ void makeCaps();
};


__device__ void SphericalFaceManager::makeCaps() {

BaseParticle* aP;
BaseParticle* bP;
Cap aCap(aP,bP);
}

Solution

  • You cannot use thrust vectors (or std::vector) directly in device code. This is mentioned in various other SO questions such as here

    If you want to use the data in a thrust::device_vector in device code, you should pass a pointer to the data as a functor initializing parameter. Various other SO questions give examples of this, such as here

    Likewise, you cannot use vector methods, e.g. .empty() or .push_back() in device code.

    You will need to replace these with ordinary C-style allocators and C-style indexed data access.

    For a multi-threaded implementation of push_back in device code, I would recommend something like this. That is a fully worked example that demonstrates how to allocate space for the vector and how each thread can use it for insertVertex for example.