convex_hull_bench.cpp

英特尔&#174 线程构建模块（英特尔&#174 TBB）是一个屡获殊荣的 C++ 运行时库

    FindXExtremum(const pointVec_t& points_, extremumType exType_)
        : points(points_), exType(exType_), extrXPoint(points[0]) {}

    FindXExtremum(const FindXExtremum& fxex, tbb::split)
        : points(fxex.points), exType(fxex.exType), extrXPoint(fxex.extrXPoint) {}

    void operator()(const range_t& range) {
        const size_t i_end = range.end();
        if(!range.empty()) {
            for(size_t i = range.begin(); i != i_end; ++i) {
                if(closerToExtremum(points[i])) {
                    extrXPoint = points[i];
                }
            }
        }
    }

    void join(const FindXExtremum &rhs) {
        if(closerToExtremum(rhs.extrXPoint)) {
            extrXPoint = rhs.extrXPoint;
        }
    }

    point_t extremeXPoint() {
        return extrXPoint;
    }

private:
    const pointVec_t    &points;
    const extremumType   exType;
    point_t              extrXPoint;
    bool closerToExtremum(const point_t &p) const {
        switch(exType) {
        case minX:
            return p.x<extrXPoint.x; break;
        case maxX:
            return p.x>extrXPoint.x; break;
        }
        return false; // avoid warning
    }
};

template <FindXExtremum::extremumType type>
point_t extremum(const pointVec_t &P) {
    FindXExtremum fxBody(P, type);
    tbb::parallel_reduce(range_t(0, P.size(), FindXExtremum::grainSize), fxBody);
    return fxBody.extremeXPoint();
}

class SplitByCP {
    const pointVec_t    &initialSet;
    pointVec_t          &reducedSet;
    point_t              p1, p2;
    point_t              farPoint;
    double               howFar;
public:
    static const size_t grainSize = cfg::DIVIDE_GS;
#if !USECONCVEC
    static mutex_t      pushBackMutex;
#endif // USECONCVEC

    SplitByCP( point_t _p1, point_t _p2,
        const pointVec_t &_initialSet, pointVec_t &_reducedSet)
        : p1(_p1), p2(_p2),
        initialSet(_initialSet), reducedSet(_reducedSet),
        howFar(0), farPoint(p1) {
    }

    SplitByCP( SplitByCP& sbcp, tbb::split )
        : p1(sbcp.p1), p2(sbcp.p2),
        initialSet(sbcp.initialSet), reducedSet(sbcp.reducedSet),
        howFar(0), farPoint(p1) {}

    void operator()( const range_t& range ) {
        const size_t i_end = range.end();
        double cp;
        for(size_t i=range.begin(); i!=i_end; ++i) {
            if( (initialSet[i] != p1) && (initialSet[i] != p2) ) {
                cp = util::cross_product(p1, p2, initialSet[i]);
                if(cp>0) {
#if USECONCVEC
                    reducedSet.push_back(initialSet[i]);
#else // Locked push_back to a not thread-safe STD::VECTOR
                    {
                        mutex_t::scoped_lock lock(pushBackMutex);
                        reducedSet.push_back(initialSet[i]);
                    }
#endif // USECONCVEC
                    if(cp>howFar) {
                        farPoint = initialSet[i];
                        howFar   = cp;
                    }
                }
            }
        }
    }

    void join(const SplitByCP& rhs) {
        if(rhs.howFar>howFar) {
            howFar   = rhs.howFar;
            farPoint = rhs.farPoint;
        }
    }

    point_t farthestPoint() const {
        return farPoint;
    }
};

class SplitByCP_buf {
    const pointVec_t    &initialSet;
    pointVec_t          &reducedSet;
    point_t              p1, p2;
    point_t              farPoint;
    double               howFar;
public:
    static const size_t  grainSize = cfg::DIVIDE_GS;
#if !USECONCVEC
    static mutex_t       insertMutex;
#endif // USECONCVEC

    SplitByCP_buf( point_t _p1, point_t _p2,
        const pointVec_t &_initialSet, pointVec_t &_reducedSet)
        : p1(_p1), p2(_p2),
        initialSet(_initialSet), reducedSet(_reducedSet),
        howFar(0), farPoint(p1) {}

    SplitByCP_buf(SplitByCP_buf& sbcp, tbb::split)
        : p1(sbcp.p1), p2(sbcp.p2),
        initialSet(sbcp.initialSet), reducedSet(sbcp.reducedSet),
        howFar(0), farPoint(p1) {}

    void operator()(const range_t& range) {
        const size_t i_end = range.end();
        size_t j = 0;
        double cp;
        point_t tmp_vec[grainSize];
        for(size_t i = range.begin(); i != i_end; ++i) {
            if( (initialSet[i] != p1) && (initialSet[i] != p2) ) {            
                cp = util::cross_product(p1, p2, initialSet[i]);
                if(cp>0) {
                    tmp_vec[j++] = initialSet[i];
                    if(cp>howFar) {
                        farPoint = initialSet[i];
                        howFar   = cp;
                    }
                }
            }
        }

#if USECONCVEC
        appendVector(tmp_vec, j, reducedSet);
#else // USE STD::VECTOR
        appendVector(insertMutex, tmp_vec, j, reducedSet);
#endif // USECONCVEC
    }

    void join(const SplitByCP_buf& rhs) {
        if(rhs.howFar>howFar) {
            howFar   = rhs.howFar;
            farPoint = rhs.farPoint;
        }
    }

    point_t farthestPoint() const {
        return farPoint;
    }
};

#if !USECONCVEC
mutex_t SplitByCP::pushBackMutex   = mutex_t();
mutex_t SplitByCP_buf::insertMutex = mutex_t();
#endif

template <typename BodyType>
point_t divide(const pointVec_t &P, pointVec_t &P_reduced,
              const point_t &p1, const point_t &p2) {
    BodyType body(p1, p2, P, P_reduced);
    tbb::parallel_reduce(range_t(0, P.size(), BodyType::grainSize), body);

    if(util::VERBOSE) {
        std::stringstream ss;
        ss << P.size() << " nodes in bucket"<< ", "
            << "dividing by: [ " << p1 << ", " << p2 << " ], "
            << "farthest node: " << body.farthestPoint();
        util::OUTPUT.push_back(ss.str());
    }

    return body.farthestPoint();
}

void divide_and_conquer(const pointVec_t &P, pointVec_t &H,
                        point_t p1, point_t p2, bool buffered) {
    if (P.size()<2) {
        H.push_back(p1);
#if USECONCVEC
        appendVector(P, H);
#else // insert into STD::VECTOR
        H.insert(H.end(), P.begin(), P.end());
#endif
    }
    else {
        pointVec_t P_reduced;
        pointVec_t H1, H2;
        point_t p_far;

        if(buffered) {
            p_far = divide<SplitByCP_buf>(P, P_reduced, p1, p2);
        } else {
            p_far = divide<SplitByCP>(P, P_reduced, p1, p2);
        }

        divide_and_conquer(P_reduced, H1, p1, p_far, buffered);
        divide_and_conquer(P_reduced, H2, p_far, p2, buffered);

#if USECONCVEC
        appendVector(H1, H);
        appendVector(H2, H);
#else // insert into STD::VECTOR
        H.insert(H.end(), H1.begin(), H1.end());
        H.insert(H.end(), H2.begin(), H2.end());
#endif
    }
}

void quickhull(const pointVec_t &points, pointVec_t &hull, bool buffered) {
    hull.clear();

    point_t p_maxx = extremum<FindXExtremum::maxX>(points);
    point_t p_minx = extremum<FindXExtremum::minX>(points);

    pointVec_t H;

    divide_and_conquer(points, hull, p_maxx, p_minx, buffered);
    divide_and_conquer(points, H, p_minx, p_maxx, buffered);
#if USECONCVEC
    appendVector(H, hull);
#else // STD::VECTOR
    hull.insert(hull.end(), H.begin(), H.end());
#endif // USECONCVEC
}

int main(int argc, char* argv[]) {
    util::ParseInputArgs(argc, argv);

    pointVec_t      points;
    pointVec_t      hull;
    size_t          nthreads;
    util::my_time_t tm_init, tm_start, tm_end;
    pointVec_t      tmp_points;

#if USECONCVEC
    std::cout << "Starting TBB unbufferred push_back version of QUICK HULL algorithm" << std::endl;
#else
    std::cout << "Starting STL locked unbufferred push_back version of QUICK HULL algorithm" << std::endl;
#endif // USECONCVEC

    for(nthreads=cfg::NUM_THREADS_START; nthreads<=cfg::NUM_THREADS_END;
        ++nthreads) {
        tbb::task_scheduler_init init(nthreads);
#if INIT_ONCE
        if(nthreads==cfg::NUM_THREADS_START) {
            tm_init = util::gettime();
            initialize<FillRNDPointsVector>(points);
        }
        else /* timing generation for stats, but use original data set */ {
            tm_init = util::gettime();
            initialize<FillRNDPointsVector>(tmp_points);
        }
#else
        tm_init = util::gettime();
        initialize<FillRNDPointsVector>(points);
#endif // INIT_ONCE
        tm_start = util::gettime();
        quickhull(points, hull, false);
        tm_end = util::gettime();

        util::WriteResults(nthreads, util::time_diff(tm_init, tm_start),
            util::time_diff(tm_start, tm_end));
    }

#if USECONCVEC 
    std::cout << "Starting TBB bufferred version of QUICK HULL algorithm" << std::endl;
#else
    std::cout << "Starting STL locked bufferred version of QUICK HULL algorithm" << std::endl;
#endif

    for(nthreads=cfg::NUM_THREADS_START; nthreads<=cfg::NUM_THREADS_END;
        ++nthreads) {
        tbb::task_scheduler_init init(nthreads);
#if INIT_ONCE
        if(nthreads==cfg::NUM_THREADS_START) {
            tm_init = util::gettime();
            initialize<FillRNDPointsVector_buf>(points);
        }
        else /* timing generation for stats, but use original data set */ {
            tm_init = util::gettime();
            initialize<FillRNDPointsVector_buf>(tmp_points);
        }
#else
        tm_init = util::gettime();
        initialize<FillRNDPointsVector_buf>(points);
#endif // INIT_ONCE
        tm_start = util::gettime();
        quickhull(points, hull, true);
        tm_end = util::gettime();

        util::WriteResults(nthreads, util::time_diff(tm_init, tm_start),
            util::time_diff(tm_start, tm_end));
    }    

    return 0;
}

#endif // USETBB