text:
#pragma once

#include <cblas.h>
#include <stdlib.h>
#include <immintrin.h>

#include <faiss/IndexIVF.h>

#include <faiss/utils/Heap.h>

#ifdef OPT_DTYPE_UTILS

namespace faiss {

//==================================Convertion================================

inline const float* convert_x_T_impl (size_t, const float* x, float*) {
    return x;
}

template <typename T>
const T* convert_x_T_impl (size_t d, const float* x, T*) {
    T* conv_x = new T[d];
    for (size_t i = 0; i < d; i++) {
        conv_x[i] = static_cast<T> (x[i]);
    }
    return conv_x;
}

template <typename T>
inline const T* convert_x_T (size_t d, const float* x) {
    return convert_x_T_impl (d, x, (T*)nullptr);
}

inline void del_converted_x_T (size_t, const float*) {
}

template <typename T>
inline void del_converted_x_T (size_t, const T* conv_x) {
    delete[] conv_x;
}

template <typename T>
struct Converter_T {

    const size_t d;
    const T* const x;

    Converter_T (size_t d, const float* x):
            d (d), x (convert_x_T<T> (d, x)) {
    }

    ~Converter_T () {
        del_converted_x_T (d, x);
    }

};

//==============================Distance Function=============================

template <typename Tdis, typename T>
Tdis vec_IP_ref_T (const T* x, const T* y, size_t d, Tdis sum = 0) {
    for (size_t i = 0; i < d; i++) {
        sum += static_cast<Tdis> (x[i]) * static_cast<Tdis> (y[i]);
    }
    return sum;
}

inline float vec_IP_ref_T (const float* x, const float* y, size_t d) {
    return vec_IP_ref_T<float> (x, y, d);
}

inline float vec_IP_ref_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_IP_ref_T<float> (x, y, d);
}

template <typename Tdis, typename T>
Tdis vec_L2Sqr_ref_T (const T* x, const T* y, size_t d, Tdis sum = 0) {
    for (size_t i = 0; i < d; i++) {
        Tdis diff = static_cast<Tdis> (x[i]) - static_cast<Tdis> (y[i]);
        sum += diff * diff;
    }
    return sum;
}

inline float vec_L2Sqr_ref_T (const float* x, const float* y, size_t d) {
    return vec_L2Sqr_ref_T<float> (x, y, d);
}

inline float vec_L2Sqr_ref_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_L2Sqr_ref_T<float> (x, y, d);
}

#ifdef __SSE4_1__

#define USE_SIMD_128

inline __m128 _mm_loadu_ps_T (const float* x) {
    return _mm_loadu_ps (x);
}

inline __m128 _mm_loadu_ps_T (const bfp16_t* x) {
    return _mm_castsi128_ps (_mm_unpacklo_epi16 (
            _mm_setzero_si128 (),
            _mm_loadl_epi64 ((const __m128i*)x)));
}

template <typename T>
float vec_IP_fp_128b_T (const T* x, const T* y, size_t d,
        __m128 msum = _mm_setzero_ps ()) {
    while (d >= 4) {
        __m128 mx = _mm_loadu_ps_T (x);
        x += 4;
        __m128 my = _mm_loadu_ps_T (y);
        y += 4;
        msum = _mm_add_ps (msum, _mm_mul_ps (mx, my));
        d -= 4;
    }
    msum = _mm_hadd_ps (msum, msum);
    msum = _mm_hadd_ps (msum, msum);
    float sum = _mm_cvtss_f32 (msum);
    return d == 0 ? sum : vec_IP_ref_T<float> (x, y, d, sum);
}

inline float vec_IP_128b_T (const float* x, const float* y, size_t d) {
    return vec_IP_fp_128b_T (x, y, d);
}

inline float vec_IP_128b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_IP_fp_128b_T (x, y, d);
}

template <typename T>
float vec_L2Sqr_fp_128b_T (const T* x, const T* y, size_t d,
        __m128 msum = _mm_setzero_ps ()) {
    while (d >= 4) {
        __m128 mx = _mm_loadu_ps_T (x);
        x += 4;
        __m128 my = _mm_loadu_ps_T (y);
        y += 4;
        __m128 mdiff = _mm_sub_ps (mx, my);
        msum = _mm_add_ps (msum, _mm_mul_ps (mdiff, mdiff));
        d -= 4;
    }
    msum = _mm_hadd_ps (msum, msum);
    msum = _mm_hadd_ps (msum, msum);
    float sum = _mm_cvtss_f32 (msum);
    return d == 0 ? sum : vec_L2Sqr_ref_T<float> (x, y, d, sum);
}

inline float vec_L2Sqr_128b_T (const float* x, const float* y, size_t d) {
    return vec_L2Sqr_fp_128b_T (x, y, d);
}

inline float vec_L2Sqr_128b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_L2Sqr_fp_128b_T (x, y, d);
}

#endif

#ifdef __AVX2__

#ifndef USE_SIMD_128
#error "SIMD 256 must have SIMD 128 enabled"
#endif

#define USE_SIMD_256

inline __m256 _mm256_loadu_ps_T (const float* x) {
    return _mm256_loadu_ps (x);
}

inline __m256 _mm256_loadu_ps_T (const bfp16_t* x) {
    return
        _mm256_castsi256_ps (
            _mm256_unpacklo_epi16 (
                _mm256_setzero_si256 (),
                _mm256_insertf128_si256 (
                    _mm256_castsi128_si256 (
                        _mm_loadl_epi64 ((const __m128i*)x)),
                    _mm_loadl_epi64 ((const __m128i*)(x + 4)),
                    1)));
}

template <typename T>
float vec_IP_fp_256b_T (const T* x, const T* y, size_t d,
        __m256 msum = _mm256_setzero_ps ()) {
    while (d >= 8) {
        __m256 mx = _mm256_loadu_ps_T (x);
        x += 8;
        __m256 my = _mm256_loadu_ps_T (y);
        y += 8;
        msum = _mm256_add_ps (msum, _mm256_mul_ps (mx, my));
        d -= 8;
    }
    __m128 msum2 = _mm256_extractf128_ps (msum, 1);
    msum2 = _mm_add_ps (msum2, _mm256_extractf128_ps (msum, 0));
    return vec_IP_fp_128b_T (x, y, d, msum2);
}

inline float vec_IP_256b_T (const float* x, const float* y, size_t d) {
    return vec_IP_fp_256b_T (x, y, d);
}

inline float vec_IP_256b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_IP_fp_256b_T (x, y, d);
}

template <typename T>
float vec_L2Sqr_fp_256b_T (const T* x, const T* y, size_t d,
        __m256 msum = _mm256_setzero_ps ()) {
    while (d >= 8) {
        __m256 mx = _mm256_loadu_ps_T (x);
        x += 8;
        __m256 my = _mm256_loadu_ps_T (y);
        y += 8;
        __m256 mdiff = _mm256_sub_ps (mx, my);
        msum = _mm256_add_ps (msum, _mm256_mul_ps (mdiff, mdiff));
        d -= 8;
    }
    __m128 msum2 = _mm256_extractf128_ps (msum, 1);
    msum2 = _mm_add_ps (msum2, _mm256_extractf128_ps (msum, 0));
    return vec_L2Sqr_fp_128b_T (x, y, d, msum2);
}

inline float vec_L2Sqr_256b_T (const float* x, const float* y, size_t d) {
    return vec_L2Sqr_fp_256b_T (x, y, d);
}

inline float vec_L2Sqr_256b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_L2Sqr_fp_256b_T (x, y, d);
}

#endif

#if defined (__AVX512F__) && defined(__AVX512DQ__) && defined(__AVX512BW__) \
        && defined(__AVX512VL__)

#ifndef USE_SIMD_256
#error "SIMD 512 must have SIMD 256 enabled"
#endif
#define USE_SIMD_512

inline __m512 _mm512_loadu_ps_T (const float* x) {
    return _mm512_loadu_ps (x);
}

inline __m512 _mm512_loadu_ps_T (const bfp16_t* x) {
    return
    _mm512_castsi512_ps (
        _mm512_unpacklo_epi16 (
            _mm512_setzero_si512 (),
            _mm512_inserti64x4 (
                _mm512_castsi256_si512 (
                    _mm256_inserti32x4 (
                        _mm256_castsi128_si256 (
                            _mm_loadl_epi64 ((const __m128i*)x)),
                        _mm_loadl_epi64 ((const __m128i*)(x + 4)),
                        1)),
                _mm256_inserti32x4 (
                    _mm256_castsi128_si256 (
                        _mm_loadl_epi64 ((const __m128i*)(x + 8))),
                    _mm_loadl_epi64 ((const __m128i*)(x + 12)),
                    1),
                1)));
}

template <typename T>
float vec_IP_fp_512b_T (const T* x, const T* y, size_t d,
        __m512 msum = _mm512_setzero_ps ()) {
    while (d >= 16) {
        __m512 mx = _mm512_loadu_ps_T (x);
        x += 16;
        __m512 my = _mm512_loadu_ps_T (y);
        y += 16;
        msum = _mm512_add_ps (msum, _mm512_mul_ps (mx, my));
        d -= 16;
    }
    __m256 msum2 = _mm512_extractf32x8_ps (msum, 1);
    msum2 = _mm256_add_ps (msum2, _mm512_extractf32x8_ps (msum, 0));
    return vec_IP_fp_256b_T (x, y, d, msum2);
}

inline float vec_IP_512b_T (const float* x, const float* y, size_t d) {
    return vec_IP_fp_512b_T (x, y, d);
}

inline float vec_IP_512b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_IP_fp_512b_T (x, y, d);
}

template <typename T>
float vec_L2Sqr_fp_512b_T (const T* x, const T* y, size_t d,
        __m512 msum = _mm512_setzero_ps ()) {
    while (d >= 16) {
        __m512 mx = _mm512_loadu_ps_T (x);
        x += 16;
        __m512 my = _mm512_loadu_ps_T (y);
        y += 16;
        __m512 mdiff = _mm512_sub_ps (mx, my);
        msum = _mm512_add_ps (msum, _mm512_mul_ps (mdiff, mdiff));
        d -= 16;
    }
    __m256 msum2 = _mm512_extractf32x8_ps (msum, 1);
    msum2 = _mm256_add_ps (msum2, _mm512_extractf32x8_ps (msum, 0));
    return vec_L2Sqr_fp_256b_T (x, y, d, msum2);
}

inline float vec_L2Sqr_512b_T (const float* x, const float* y, size_t d) {
    return vec_L2Sqr_fp_512b_T (x, y, d);
}

inline float vec_L2Sqr_512b_T (const bfp16_t* x, const bfp16_t* y,
        size_t d) {
    return vec_L2Sqr_fp_512b_T (x, y, d);
}

#endif

#if defined (USE_SIMD_512)

template <typename T>
inline float vec_IP_T (const T* x, const T* y, size_t d) {
    return vec_IP_512b_T (x, y, d);
}

template <typename T>
inline float vec_L2Sqr_T (const T* x, const T* y, size_t d) {
    return vec_L2Sqr_512b_T (x, y, d);
}

#elif defined (USE_SIMD_256)

template <typename T>
inline float vec_IP_T (const T* x, const T* y, size_t d) {
    return vec_IP_256b_T (x, y, d);
}

template <typename T>
inline float vec_L2Sqr_T (const T* x, const T* y, size_t d) {
    return vec_L2Sqr_256b_T (x, y, d);
}

#elif defined (USE_SIMD_128)

template <typename T>
inline float vec_IP_T (const T* x, const T* y, size_t d) {
    return vec_IP_128b_T (x, y, d);
}

template <typename T>
inline float vec_L2Sqr_T (const T* x, const T* y, size_t d) {
    return vec_L2Sqr_128b_T (x, y, d);
}

#else

template <typename T>
inline float vec_IP_T (const T* x, const T* y, size_t d) {
    return vec_IP_ref_T (x, y, d);
}

template <typename T>
inline float vec_L2Sqr_T (const T* x, const T* y, size_t d) {
    return vec_L2Sqr_ref_T (x, y, d);
}

#endif

}

#endif

#ifdef OPT_FLAT_DTYPE

#define FLAT_BATCH_THRESHOLD    4

namespace faiss {

//=================================KNN Routine================================

template <typename T, typename D>
void knn_less_better_alone_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res, D& distance) {
    size_t k = res->k;
    size_t check_period = InterruptCallback::get_period_hint (ny * d);
    check_period *= omp_get_max_threads ();
    for (size_t i0 = 0; i0 < nx; i0 += check_period) {
        size_t i1 = std::min (i0 + check_period, nx);
#pragma omp parallel for
        for (size_t i = i0; i < i1; i++) {
            const T* x_i = x + i * d;
            const T* y_j = y;
            float* simi = res->get_val (i);
            int64_t* idxi = res->get_ids (i);
            maxheap_heapify (k, simi, idxi);
            for (size_t j = 0; j < ny; j++) {
                float dis = distance (i, j, x_i, y_j, d);
                if (dis < simi[0]) {
                    maxheap_pop (k, simi, idxi);
                    maxheap_push (k, simi, idxi, dis, j);
                }
                y_j += d;
            }
            maxheap_reorder (k, simi, idxi);
        }
        InterruptCallback::check ();
    }
}

template <typename T, typename D>
void knn_greater_better_alone_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_minheap_array_t* res, D& distance) {
    size_t k = res->k;
    size_t check_period = InterruptCallback::get_period_hint (ny * d);
    check_period *= omp_get_max_threads ();
    for (size_t i0 = 0; i0 < nx; i0 += check_period) {
        size_t i1 = std::min (i0 + check_period, nx);
#pragma omp parallel for
        for (size_t i = i0; i < i1; i++) {
            const T* x_i = x + i * d;
            const T* y_j = y;
            float* simi = res->get_val (i);
            int64_t* idxi = res->get_ids (i);
            minheap_heapify (k, simi, idxi);
            for (size_t j = 0; j < ny; j++) {
                float dis = distance (i, j, x_i, y_j, d);
                if (dis > simi[0]) {
                    minheap_pop (k, simi, idxi);
                    minheap_push (k, simi, idxi, dis, j);
                }
                y_j += d;
            }
            minheap_reorder (k, simi, idxi);
        }
        InterruptCallback::check ();
    }
}

template <typename T>
inline void knn_inner_product_alone_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_minheap_array_t* res) {

    struct IP {

        inline float operator () (size_t /*ix*/, size_t /*jy*/,
                const T* xi, const T* yj, size_t d) const {
            return vec_IP_T (xi, yj, d);
        }

    }
    distance;
    knn_greater_better_alone_T (x, y, d, nx, ny, res, distance);
}

template <typename T>
inline void knn_L2Sqr_alone_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res) {

    struct L2Sqr {

        inline float operator () (size_t /*ix*/, size_t /*jy*/,
                const T* xi, const T* yj, size_t d) const {
            return vec_L2Sqr_T (xi, yj, d);
        }

    }
    distance;
    knn_less_better_alone_T (x, y, d, nx, ny, res, distance);
}

template <typename T>
inline void knn_L2Sqr_expand_alone_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res,
        const float* y_norm) {

    struct L2SqrExpand {

        const float* y_norm_sqr;

        inline float operator () (size_t /*ix*/, size_t jy,
                const T* xi, const T* yj, size_t d) const {
            return y_norm_sqr[jy] - 2 * vec_IP_T (xi, yj, d);
        }

    }
    distance = {
        .y_norm_sqr = y_norm,
    };
    knn_less_better_alone_T (x, y, d, nx, ny, res, distance);
}

template <typename T>
inline void knn_inner_product_batch_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_minheap_array_t* res) {
    knn_inner_product_alone_T (x, y, d, nx, ny, res);
}

template <typename T>
inline void knn_L2Sqr_batch_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res) {
    knn_L2Sqr_alone_T (x, y, d, nx, ny, res);
}

template <typename T>
inline void knn_L2Sqr_expand_batch_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res,
        const float* y_norm) {
    knn_L2Sqr_expand_alone_T (x, y, d, nx, ny, res, y_norm);
}

template <typename H, typename D>
void knn_batch_T (const float* x, const float* y,
        size_t d, size_t nx, size_t ny, H* heap, D& distance) {
    heap->heapify ();
    if (nx == 0 || ny == 0) {
        return;
    }
    float* distances = new float [nx * ny];
    distance (x, y, d, nx, ny, distances);
    heap->addn (ny, distances, 0, 0, nx);
    delete[] distances;
    InterruptCallback::check ();
    heap->reorder ();
}

inline void knn_inner_product_batch_T (const float* x, const float* y,
        size_t d, size_t nx, size_t ny, float_minheap_array_t* res) {

    struct IP {

        inline void operator () (const float* x, const float* y,
                size_t d, size_t nx, size_t ny, float* distances) const {
            cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasTrans, nx, ny, d,
                    1.0f, x, d, y, d, 0.0f, distances, ny);
        }

    }
    distance;
    knn_batch_T (x, y, d, nx, ny, res, distance);
}

inline void knn_L2Sqr_expand_batch_T (const float* x, const float* y,
        size_t d, size_t nx, size_t ny, float_maxheap_array_t* res,
        const float* y_norm) {

    struct L2SqrExpand {

        const float* y_norm;

        inline void operator () (const float* x, const float* y,
                size_t d, size_t nx, size_t ny, float* distances) const {
            float* distances_i = distances;
            size_t step = ny * sizeof(float);
            for (size_t i = 0; i < nx; i++) {
                memcpy (distances_i, y_norm, step);
                distances_i += ny;
            }
            cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasTrans, nx, ny, d,
                    -2.0f, x, d, y, d, 1.0f, distances, ny);
            }
    }
    distance = {
        .y_norm = y_norm,
    };
    knn_batch_T (x, y, d, nx, ny, res, distance);
}

template <typename T>
inline void knn_inner_product_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_minheap_array_t* res) {
    if (nx < FLAT_BATCH_THRESHOLD) {
        knn_inner_product_alone_T (x, y, d, nx, ny, res);
    }
    else {
        knn_inner_product_batch_T (x, y, d, nx, ny, res);
    }
}

template <typename T>
inline void knn_L2Sqr_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res) {
    if (nx < FLAT_BATCH_THRESHOLD) {
        knn_L2Sqr_alone_T (x, y, d, nx, ny, res);
    }
    else {
        knn_L2Sqr_batch_T (x, y, d, nx, ny, res);
    }
}

template <typename T>
inline void knn_L2Sqr_expand_T (const T* x, const T* y, size_t d,
        size_t nx, size_t ny, float_maxheap_array_t* res,
        const float* y_norm) {
    if (nx < FLAT_BATCH_THRESHOLD) {
        knn_L2Sqr_expand_alone_T (x, y, d, nx, ny, res, y_norm);
    }
    else {
        knn_L2Sqr_expand_batch_T (x, y, d, nx, ny, res, y_norm);
    }
}

}

#endif

#ifdef OPT_IVFFLAT_DTYPE

#define SCANNER_USE_BATCH       false

namespace faiss {

//===========================Inverted List Scanner============================

template <typename T>
class InvertedListScanner_T : public InvertedListScanner {

    using idx_t = InvertedListScanner::idx_t;

protected:
    size_t d;
    size_t code_size;
    bool store_pairs;
    const T* converted_x;
    idx_t list_no;

public:
    InvertedListScanner_T (size_t d, bool store_pairs):
            d (d), code_size (sizeof(T) * d), store_pairs (store_pairs),
            converted_x(nullptr), list_no(-1) {
    }

    virtual ~InvertedListScanner_T () {
        if (converted_x) {
            del_converted_x_T (d, converted_x);
        }
    }

    virtual void set_query (const float* query) override {
        if (converted_x) {
            del_converted_x_T (d, converted_x);
        }
        converted_x = convert_x_T<T> (d, query);
    }

    virtual void set_list (idx_t lidx, float) override {
        list_no = lidx;
    }

    virtual float distance_to_code (const uint8_t*) const override {
        FAISS_THROW_MSG ("not implemented");
    }

    virtual size_t scan_codes (size_t list_size, const uint8_t* codes,
            const idx_t* ids, float* simi, idx_t* idxi, size_t k)
            const = 0;

};

template <typename T, typename C, typename D>
class AloneInvertedListScanner_T : public InvertedListScanner_T<T> {

    using idx_t = InvertedListScanner::idx_t;
    using Scanner = InvertedListScanner_T<T>;

private:
    D* distance;

public:
    AloneInvertedListScanner_T (size_t d, bool store_pairs, D* distance):
            Scanner (d, store_pairs), distance (distance) {
    }

    virtual ~AloneInvertedListScanner_T () {
        delete distance;
    }

    virtual size_t scan_codes (size_t list_size, const uint8_t* codes,
            const idx_t* ids, float* simi, idx_t* idxi,
            size_t k) const override {
        size_t nup = 0;
        for (size_t i = 0; i < list_size; i++) {
            float dis = (*distance) (Scanner::list_no, i,Scanner::converted_x,
                    (const T*)codes, Scanner::d);
            codes += Scanner::code_size;
            if (C::cmp (simi[0], dis)) {
                heap_pop<C> (k, simi, idxi);
                int64_t id = Scanner::store_pairs ?
                        lo_build (Scanner::list_no, i) :
                        ids[i];
                heap_push<C> (k, simi, idxi, dis, id);
                nup++;
            }
        }
        return nup;
    }

};

template <typename T>
InvertedListScanner* get_IP_alone_scanner_T (size_t d, bool store_pairs) {

    struct IP {

        inline float operator () (size_t /*ilist*/, size_t /*jy*/,
                const T* x, const T* yj, size_t d) const {
            return vec_IP_T (x, yj, d);
        }

    }
    *distance = new IP;
    return new AloneInvertedListScanner_T<T, CMin<float, int64_t>, IP> (d,
            store_pairs, distance);
}

template <typename T>
InvertedListScanner* get_L2Sqr_alone_scanner_T (size_t d, bool store_pairs) {

    struct L2Sqr {

        inline float operator () (size_t /*ilist*/, size_t /*jy*/,
                const T* x, const T* yj, size_t d) const {
            return vec_L2Sqr_T (x, yj, d);
        }

    }
    *distance = new L2Sqr;
    return new AloneInvertedListScanner_T<T, CMax<float, int64_t>, L2Sqr> (d,
            store_pairs, distance);
}

template <typename T, typename TNorm>
InvertedListScanner* get_L2Sqr_expand_alone_scanner_T (size_t d,
        bool store_pairs, const TNorm y_norm) {

    struct L2SqrExpand {

        const TNorm y_norm;

        inline float operator () (size_t ilist, size_t jy,
                const T* x, const T* yj, size_t d) {
            return y_norm [ilist] [jy] - 2 * vec_IP_T (x, yj, d);
        }

    }
    *distance = new L2SqrExpand {
        .y_norm = y_norm,
    };
    return new AloneInvertedListScanner_T<T, CMax<float, int64_t>,
            L2SqrExpand> (d, store_pairs, distance);
}

template <typename T, typename C, typename D>
class BatchInvertedListScanner_T : public InvertedListScanner_T<T> {

    using idx_t = InvertedListScanner::idx_t;
    using Scanner = InvertedListScanner_T<T>;

private:
    D* distance;

public:
    BatchInvertedListScanner_T (size_t d, bool store_pairs, D* distance):
            Scanner (d, store_pairs), distance (distance) {
    }

    virtual ~BatchInvertedListScanner_T () {
        delete distance;
    }

    virtual size_t scan_codes (size_t list_size, const uint8_t* codes,
            const idx_t* ids, float* simi, idx_t* idxi,
            size_t k) const override {
        float* distances = new float [list_size];
        (*distance) (Scanner::converted_x, Scanner::list_no, list_size,
                (const T*)codes, Scanner::d, distances);
        size_t nup = 0;
        for (size_t i = 0; i < list_size; i++) {
            float dis = distances [i];
            if (C::cmp (simi[0], dis)) {
                heap_pop<C> (k, simi, idxi);
                int64_t id = Scanner::store_pairs ?
                        lo_build (Scanner::list_no, i) : ids[i];
                heap_push<C> (k, simi, idxi, dis, id);
                nup++;
            }
        }
        delete[] distances;
        return nup;
    }

};

template <typename T>
inline InvertedListScanner* get_IP_batch_scanner_T (size_t d, bool store_pairs,
        T*) {
    return get_IP_alone_scanner_T<T> (d, store_pairs);
}

inline InvertedListScanner* get_IP_batch_scanner_T (size_t d, bool store_pairs,
        float*) {

    struct IP {

        inline void operator () (const float* x, size_t /*ilist*/,
                size_t list_size, const float* y, size_t d, float* distances)
                const {
            cblas_sgemv (CblasRowMajor, CblasNoTrans, list_size, d, 1.0f,
                    y, d, x, 1, 0.0f, distances, 1);
        }

    }
    *distance = new IP;
    return new BatchInvertedListScanner_T<float, CMin<float, int64_t>, IP> (d,
            store_pairs, distance);
}

template <typename T>
inline InvertedListScanner* get_L2Sqr_batch_scanner_T (size_t d,
        bool store_pairs, T*) {
    return get_L2Sqr_alone_scanner_T<T> (d, store_pairs);
}

template <typename T, typename TNorm>
inline InvertedListScanner* get_L2Sqr_expand_batch_scanner_T (size_t d,
        bool store_pairs, const TNorm y_norm, T*) {
    return get_L2Sqr_expand_alone_scanner_T<T> (d, store_pairs,
            y_norm);
}

template <typename TNorm>
inline InvertedListScanner* get_L2Sqr_expand_batch_scanner_T (size_t d,
        bool store_pairs, const TNorm y_norm, float*) {

    struct L2SqrExpand {

        const TNorm y_norm;

        inline void operator () (const float* x, size_t ilist,
                size_t list_size, const float* y, size_t d, float* distances)
                const {
            memcpy (distances, &(y_norm [ilist] [0]),
                    list_size * sizeof(float));
            cblas_sgemv (CblasRowMajor, CblasNoTrans, list_size, d, -2.0f,
                    y, d, x, 1, 1.0f, distances, 1);
        }

    }
    *distance = new L2SqrExpand {
        .y_norm = y_norm,
    };
    return new BatchInvertedListScanner_T<float, CMax<float, int64_t>,
            L2SqrExpand> (d, store_pairs, distance);
}

template <typename T>
inline InvertedListScanner* get_IP_scanner_T (size_t d, bool store_pairs) {
    if (!SCANNER_USE_BATCH) {
        return get_IP_alone_scanner_T<T> (d, store_pairs);
    }
    else {
        return get_IP_batch_scanner_T (d, store_pairs, (T*)nullptr);
    }
}

template <typename T>
inline InvertedListScanner* get_L2Sqr_scanner_T (size_t d, bool store_pairs) {
    if (!SCANNER_USE_BATCH) {
        return get_L2Sqr_alone_scanner_T<T> (d, store_pairs);
    }
    else {
        return get_L2Sqr_batch_scanner_T (d, store_pairs, (T*)nullptr);
    }
}

template <typename T, typename TNorm>
inline InvertedListScanner* get_L2Sqr_expand_scanner_T (size_t d,
        bool store_pairs, const TNorm y_norm) {
    if (!SCANNER_USE_BATCH) {
        return get_L2Sqr_expand_alone_scanner_T<T> (d, store_pairs, y_norm);
    }
    else {
        return get_L2Sqr_expand_batch_scanner_T (d, store_pairs, y_norm,
                (T*)nullptr);
    }
}

}

#endif

max_stars_repo_name: zhou-yuxin/faiss

max_issues_repo_name: zhou-yuxin/faiss

max_forks_repo_name: zhou-yuxin/faiss


text:
/* specfunc/test_mathieu.c
 *
 * Copyright (C) 2003 Lowell Johnson
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

#define NVAL 100

static double c[NVAL];

int test_mathieu(void)
{
  gsl_sf_result r;
  gsl_sf_mathieu_workspace *work;
  int s = 0;
  int sa;

  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 0.0, 0.0, &r),
          0.7071067811865475, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 0.0, M_PI_2, &r),
          0.7071067811865475, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 5.0, 0.0, &r),
          0.04480018165188902, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 5.0, M_PI_2, &r),
          1.334848674698019, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 10.0, 0.0, &r),
          0.007626517570935782, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 10.0, M_PI_2, &r),
          1.468660470712856, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 15.0, 0.0, &r),
          0.001932508315204592, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 15.0, M_PI_2, &r),
          1.550108146686649, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 20.0, 0.0, &r),
          0.0006037438292242197, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 20.0, M_PI_2, &r),
          1.609890857395926, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 25.0, 0.0, &r),
          0.0002158630184146612, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (0, 25.0, M_PI_2, &r),
          1.657510298323475, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 0.0, 0.0, &r),
          1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 5.0, 0.0, &r),
          0.2565428793223637, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 10.0, 0.0, &r),
          0.05359874774717657, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 15.0, 0.0, &r),
          0.01504006645382623, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 20.0, 0.0, &r),
          0.005051813764712904, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (1, 25.0, 0.0, &r),
          0.001911051506657645, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 0.0, M_PI_2, &r),
          1.0000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 5.0, M_PI_2, &r),
          1.337433887022345, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 10.0, M_PI_2, &r),
          1.468755664102938, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 15.0, M_PI_2, &r),
          1.550115074357552, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 20.0, M_PI_2, &r),
          1.609891592603772, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (1, 25.0, M_PI_2, &r),
          1.657510398374516, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 0.0, 0.0, &r),
          1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 0.0, M_PI_2, &r),
          -1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 5.0, 0.0, &r),
          0.7352943084006845, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 5.0, M_PI_2, &r),
          -0.7244881519676682, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 10.0, 0.0, &r),
          0.2458883492913189, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 10.0, M_PI_2, &r),
          -0.9267592641263211, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 15.0, 0.0, &r),
          0.07879282784639313, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 15.0, M_PI_2, &r),
          -1.019966226030262, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 20.0, 0.0, &r),
          0.02864894314707431, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 20.0, M_PI_2, &r),
          -1.075293228779687, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 25.0, 0.0, &r),
          0.0115128663308875, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (2, 25.0, M_PI_2, &r),
          -1.116278953295253, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 0.0, 0.0, &r),
          1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 5.0, 0.0, &r),
          1.12480725063848, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 10.0, 0.0, &r),
          1.258019941308287, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 15.0, 0.0, &r),
          1.193432230413072, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 20.0, 0.0, &r),
          0.9365755314226215, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (5, 25.0, 0.0, &r),
          0.6106943100506986, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 0.0, M_PI_2, &r),
          1.0000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 5.0, M_PI_2, &r),
          0.9060779302023551, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 10.0, M_PI_2, &r),
          0.8460384335355106, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 15.0, M_PI_2, &r),
          0.837949340012484, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 20.0, M_PI_2, &r),
          0.8635431218533667, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (5, 25.0, M_PI_2, &r),
          0.8992683245108413, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 0.0, 0.0, &r),
          1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 0.0, M_PI_2, &r),
          -1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 5.0, 0.0, &r),
          1.025995027089438, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 5.0, M_PI_2, &r),
          -0.975347487235964, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 10.0, 0.0, &r),
          1.053815992100935, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 10.0, M_PI_2, &r),
          -0.9516453181789554, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 15.0, 0.0, &r),
          1.084106311839221, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 15.0, M_PI_2, &r),
          -0.9285480638845388, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 20.0, 0.0, &r),
          1.117788631259397, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 20.0, M_PI_2, &r),
          -0.9057107845940974, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 25.0, 0.0, &r),
          1.156239918632239, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (10, 25.0, M_PI_2, &r),
          -0.8826919105636903, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 0.0, 0.0, &r),
          1.00000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 5.0, 0.0, &r),
          1.011293732529566, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 10.0, 0.0, &r),
          1.022878282438181, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 15.0, 0.0, &r),
          1.034793652236873, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 20.0, 0.0, &r),
          1.047084344162887, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_ce_e, (15, 25.0, 0.0, &r),
          1.059800441813937, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 0.0, M_PI_2, &r),
          -1.0000000, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 5.0, M_PI_2, &r),
          -0.9889607027406357, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 10.0, M_PI_2, &r),
          -0.9781423471832157, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 15.0, M_PI_2, &r),
          -0.9675137031854538, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 20.0, M_PI_2, &r),
          -0.9570452540612817, TEST_SNGL, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_mathieu_se_e, (15, 25.0, M_PI_2, &r),
          -0.9467086958780897, TEST_SNGL, GSL_SUCCESS);

  work = gsl_sf_mathieu_alloc(NVAL, 20.0);
  sa = 0;
  gsl_sf_mathieu_ce_array(0, 5, 0.0, M_PI_2, work, c);
  sa += (test_sf_frac_diff(c[0], 0.7071067811865475) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[2], -1.0) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[4], 1.0) > TEST_SNGL);
  gsl_test(sa, "gsl_sf_mathieu_ce_array");
  s += sa;

  sa = 0;
  gsl_sf_mathieu_ce_array(0, 15, 20.0, 0.0, work, c);
  sa += (test_sf_frac_diff(c[0], 0.0006037438292242197) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[1], 0.005051813764712904) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[2], 0.02864894314707431) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[5], 0.9365755314226215) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[10], 1.117788631259397) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[15], 1.047084344162887) > TEST_SNGL);
  gsl_test(sa, "gsl_sf_mathieu_ce_array");
  s += sa;

  sa = 0;
  gsl_sf_mathieu_se_array(1, 15, 20.0, M_PI_2, work, c);
  sa += (test_sf_frac_diff(c[0], 1.609891592603772) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[4], 0.8635431218533667) > TEST_SNGL);
  sa += (test_sf_frac_diff(c[14], -0.9570452540612817) > TEST_SNGL);
  gsl_test(sa, "gsl_sf_mathieu_se_array");
  s += sa;

  gsl_sf_mathieu_free(work);
  return s;
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/**
 *
 * @file qwrapper_zlaswp.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaswp(Quark *quark, Quark_Task_Flags *task_flags,
                       int n, PLASMA_Complex64_t *A, int lda,
                       int i1,  int i2, const int *ipiv, int inc)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_zlaswp_quark, task_flags,
        sizeof(int),                      &n,    VALUE,
        sizeof(PLASMA_Complex64_t)*lda*n,  A,        INOUT | LOCALITY,
        sizeof(int),                      &lda,  VALUE,
        sizeof(int),                      &i1,   VALUE,
        sizeof(int),                      &i2,   VALUE,
        sizeof(int)*n,                     ipiv,     INPUT,
        sizeof(int),                      &inc,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaswp_quark = PCORE_zlaswp_quark
#define CORE_zlaswp_quark PCORE_zlaswp_quark
#endif
void CORE_zlaswp_quark(Quark *quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    PLASMA_Complex64_t *A;

    quark_unpack_args_7(quark, n, A, lda, i1, i2, ipiv, inc);
    LAPACKE_zlaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaswp_f2(Quark *quark, Quark_Task_Flags *task_flags,
                          int n, PLASMA_Complex64_t *A, int lda,
                          int i1,  int i2, const int *ipiv, int inc,
                          PLASMA_Complex64_t *fake1, int szefake1, int flag1,
                          PLASMA_Complex64_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_zlaswp_f2_quark, task_flags,
        sizeof(int),                        &n,     VALUE,
        sizeof(PLASMA_Complex64_t)*lda*n,    A,         INOUT | LOCALITY,
        sizeof(int),                        &lda,   VALUE,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*n,                       ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(PLASMA_Complex64_t)*szefake1, fake1,     flag1,
        sizeof(PLASMA_Complex64_t)*szefake2, fake2,     flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaswp_f2_quark = PCORE_zlaswp_f2_quark
#define CORE_zlaswp_f2_quark PCORE_zlaswp_f2_quark
#endif
void CORE_zlaswp_f2_quark(Quark* quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    PLASMA_Complex64_t *A;
    void *fake1, *fake2;

    quark_unpack_args_9(quark, n, A, lda, i1, i2, ipiv, inc, fake1, fake2);
    LAPACKE_zlaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaswp_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex64_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc, PLASMA_Complex64_t *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_zlaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex64_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex64_t)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_zlaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex64_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex64_t)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaswp_ontile_quark = PCORE_zlaswp_ontile_quark
#define CORE_zlaswp_ontile_quark PCORE_zlaswp_ontile_quark
#endif
void CORE_zlaswp_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex64_t *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_zlaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaswp_ontile_f2(Quark *quark, Quark_Task_Flags *task_flags,
                                 PLASMA_desc descA, PLASMA_Complex64_t *Aij,
                                 int i1,  int i2, const int *ipiv, int inc,
                                 PLASMA_Complex64_t *fake1, int szefake1, int flag1,
                                 PLASMA_Complex64_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_zlaswp_ontile_f2_quark, task_flags,
        sizeof(PLASMA_desc),                &descA, VALUE,
        sizeof(PLASMA_Complex64_t)*1,        Aij,       INOUT | LOCALITY,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),      ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(PLASMA_Complex64_t)*szefake1, fake1, flag1,
        sizeof(PLASMA_Complex64_t)*szefake2, fake2, flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaswp_ontile_f2_quark = PCORE_zlaswp_ontile_f2_quark
#define CORE_zlaswp_ontile_f2_quark PCORE_zlaswp_ontile_f2_quark
#endif
void CORE_zlaswp_ontile_f2_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex64_t *A;
    PLASMA_desc descA;
    void *fake1, *fake2;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, fake1, fake2);
    CORE_zlaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zswptr_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex64_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc,
                              const PLASMA_Complex64_t *Akk, int ldak)
{
    DAG_CORE_TRSM;
    QUARK_Insert_Task(
        quark, CORE_zswptr_ontile_quark, task_flags,
        sizeof(PLASMA_desc),              &descA, VALUE,
        sizeof(PLASMA_Complex64_t)*1,      Aij,       INOUT | LOCALITY,
        sizeof(int),                      &i1,    VALUE,
        sizeof(int),                      &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),    ipiv,      INPUT,
        sizeof(int),                      &inc,   VALUE,
        sizeof(PLASMA_Complex64_t)*ldak,   Akk,       INPUT,
        sizeof(int),                      &ldak,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zswptr_ontile_quark = PCORE_zswptr_ontile_quark
#define CORE_zswptr_ontile_quark PCORE_zswptr_ontile_quark
#endif
void CORE_zswptr_ontile_quark(Quark *quark)
{
    int i1, i2, inc, ldak;
    int *ipiv;
    PLASMA_Complex64_t *A, *Akk;
    PLASMA_desc descA;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, Akk, ldak);
    CORE_zswptr_ontile(descA, i1, i2, ipiv, inc, Akk, ldak);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaswpc_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex64_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc, PLASMA_Complex64_t *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_zlaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex64_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex64_t)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_zlaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex64_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex64_t)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaswpc_ontile_quark = PCORE_zlaswpc_ontile_quark
#define CORE_zlaswpc_ontile_quark PCORE_zlaswpc_ontile_quark
#endif
void CORE_zlaswpc_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex64_t *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_zlaswpc_ontile(descA, i1, i2, ipiv, inc);
}


text:
//====---- Sudoku/Value.h                                             ----====//
//
// Strongly typed Value type.
// Unsigned for use as a bit-mask.
//====--------------------------------------------------------------------====//
#pragma once

#include "Size.h"

#include <gsl/gsl>

#include <vector>

#include <algorithm> // all_of
#include <limits>
#include <stdexcept>
#include <type_traits> // is_same

#include <cstddef> // size_t

namespace Sudoku
{
class Value
{
public:
	Value() noexcept = default;
	explicit constexpr Value(size_t val) noexcept : value_(val) {}

	explicit constexpr operator size_t() const noexcept { return value_; }
	explicit constexpr operator bool() const noexcept { return value_ != 0; }
	// note: static_assert performs an explicit conversion to bool

	friend constexpr bool operator==(const Value& a, const Value& b) noexcept
	{ // Friend definitions: https://abseil.io/tips/99
		return a.value_ == b.value_;
	}
	friend constexpr bool operator<(const Value& a, const Value& b) noexcept
	{
		return a.value_ < b.value_;
	}

	constexpr Value& operator++() noexcept;
	constexpr Value operator++(int) noexcept;

private:
	size_t value_{0};
};

//====--------------------------------------------------------------------====//
// Free-function declarations
constexpr bool operator!=(const Value&, const Value&) noexcept;
constexpr bool operator>(const Value&, const Value&) noexcept;
constexpr bool operator<=(const Value&, const Value&) noexcept;
constexpr bool operator>=(const Value&, const Value&) noexcept;

template<typename T, int N>
constexpr Value to_Value(T val);

template<int N>
constexpr bool is_valid(const Value&) noexcept;
template<int E>
inline constexpr bool is_valid_option(const Value&) noexcept;
template<int N>
bool is_valid(const std::vector<Value>&) noexcept;

//====--------------------------------------------------------------------====//

// Test input value
template<int N>
inline constexpr bool is_valid(const Value& value) noexcept
{
	return is_valid_option<elem_size<N>>(value);
}

// Test input value, for use with Options<E>
template<int E>
inline constexpr bool is_valid_option(const Value& value) noexcept
{
	return value > Value{0} && value <= Value{E};
}

// Test input values
template<int N>
bool is_valid(const std::vector<Value>& values) noexcept
{
	return (
		!values.empty() &&
		std::all_of(values.cbegin(), values.cend(), [](Value i) {
			return is_valid<N>(i);
		}));
}

//====--------------------------------------------------------------------====//
// Checked conversion to Value
template<int N, typename T>
inline constexpr Value to_Value(T val)
{
	if constexpr (std::is_same_v<Value, T>)
	{
		if (val > static_cast<Value>(elem_size<N>))
		{
			throw std::domain_error("Value input too large");
		}
		return val;
	}
	else if constexpr (std::is_unsigned_v<T>)
	{
		static_assert(elem_size<N> < std::numeric_limits<T>::max());
		return to_Value<N>(static_cast<Value>(val));
	}
	else
	{
		static_assert(std::is_integral_v<T>);
		static_assert(std::numeric_limits<T>::max() > elem_size<N>);
		if (val < 0)
		{
			throw std::domain_error("Value can not be negative");
		}
		return to_Value<N>(gsl::narrow_cast<size_t>(val));
	}
}

//====--------------------------------------------------------------------====//
constexpr Value& Value::operator++() noexcept
{ // pre-increment
	++value_;
	return *this;
}
// NOLINTNEXTLINE(readability/casting) false-positive cpplint
constexpr Value Value::operator++(int) noexcept
{ // post-increment
	const Value pre{*this};
	operator++();
	return pre;
}

inline constexpr bool operator!=(const Value& left, const Value& right) noexcept
{
	return !(left == right);
}
inline constexpr bool operator<=(const Value& left, const Value& right) noexcept
{
	return !(right < left);
}
inline constexpr bool operator>(const Value& left, const Value& right) noexcept
{
	return right < left;
}
inline constexpr bool operator>=(const Value& left, const Value& right) noexcept
{
	return !(left < right);
}
static_assert(Value{4} == Value{4}); //-V501 //-V112
static_assert(Value{1} != Value{0});
static_assert(Value{3} > Value{2});
// NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
static_assert(Value{10} >= Value{7}); // NOLINT(readability-magic-numbers)
// NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
static_assert(Value{8} <= Value{9}); // NOLINT(readability-magic-numbers)
// NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
static_assert(Value{5} < Value{6}); // NOLINT(readability-magic-numbers)

} // namespace Sudoku

ext: h

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
#include <cblas.h>

#include "tasks.h"

void gemm_task_seq(void *ptr)
{
    struct gemm_task_arg *arg = (struct gemm_task_arg*) ptr;

    int m        = arg->m;
    int n        = arg->n;
    int k        = arg->k;
    double *A21  = arg->A21;
    double *A21T = arg->A21T;
    double *A22  = arg->A22;
    int ldA      = arg->ldA;

    // Compute A22 := A22 - A21 * A21T'.
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
                m, n, k,
                -1.0, A21,  ldA,
                      A21T, ldA,
                 1.0, A22,  ldA);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/*--------------------------------------------------------------------
 * $Id$
 *
 * This file is part of libRadtran.
 * Copyright (c) 1997-2012 by Arve Kylling, Bernhard Mayer,
 *                            Claudia Emde, Robert Buras
 *
 * ######### Contact info: http://www.libradtran.org #########
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 *--------------------------------------------------------------------*/

#include <unistd.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include "errors.h"
#include "uvspecrandom.h"

#if HAVE_LIBGSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#endif

/* global variable: random number generator */
#if HAVE_LIBGSL
gsl_rng* uvspecrng;
#endif

/*****************************************************************/
/* uvspec random number generator; either the highly-recommended */
/* MT19937 if GSL is available, or random() otherwise.           */
/*****************************************************************/

double uvspec_random() {
#if HAVE_LIBGSL
  return gsl_rng_uniform (uvspecrng);
#else
  return (double)random();
#endif
}

double uvspec_random_gauss (double sigma) {
#if HAVE_LIBGSL
  return gsl_ran_gaussian (uvspecrng, sigma);
#else
  return (double)random();
#endif
}

int init_uvspec_random (int* randomseed, int need_mt19937, int quiet) {
  int rseed = 0;

/* global variable, defined in uvspec_lex.l */
#if HAVE_LIBGSL
  /*  extern gsl_rng *uvspecrng; */ /* global random number generator */
  const gsl_rng_type* T;
#endif

  /* always initialize standard rng */
  srandom (rseed);

  if (*randomseed > 0)
    rseed = *randomseed;
  else
    rseed = (int)time (NULL) + (int)getpid();

#if HAVE_LIBGSL
  T                    = gsl_rng_mt19937;
  gsl_rng_default_seed = rseed;
  uvspecrng            = gsl_rng_alloc (T);
  if (!quiet && need_mt19937)
    fprintf (stderr, " ... using %s random number generator\n", gsl_rng_name (uvspecrng));
#else
  if (need_mt19937) {
    fprintf (stderr, "Error, GSL not found. The MT19937 random number generator is\n");
    fprintf (stderr, "required for MYSTIC simulations. Please install GSL and recompile!\n");
    return -1;
  }
#endif

  return rseed;
}

/***********************************************************************************/
/* Function: read_and_write_random_status                                 @62_30i@ */
/* Description:                                                                    */
/*  check whether to read/write random status.                                     */
/* Parameters:                                                                     */
/* Return value:                                                                   */
/* Example:                                                                        */
/* Files:                                                                          */
/* Known bugs:                                                                     */
/* Author: RPB                                                                     */
/*                                                                        @i62_30@ */
/***********************************************************************************/

#if HAVE_LIBGSL
int read_and_write_random_status (char* randomstatusfilename,
                                  int*  ili,
                                  int*  backward_is,
                                  int*  backward_js,
                                  int   photoncounter,
                                  int*  readrandomstatus) {
  FILE* randomstatusfile = NULL;
  int   status           = 0;

  /************************************+**************************/
  /* starting with the randomseed is nice, but starting with the */
  /* actual photon that produced the problem is even nicer...    */
  /* here we read and write the status of the random number      */
  /* generator                                                   */
  /* this only works with the gsl !!!                            */
  /************************************+**************************/

  /* read randomstatusfile if asked for */
  if (*readrandomstatus && photoncounter == 1) {
    if ((randomstatusfile = fopen (randomstatusfilename, "r")) != NULL) {
      status = gsl_rng_fread (randomstatusfile, uvspecrng);
      if (status)
        return err_out ("Error %d returned by gsl_rng_fread()\n", status);
      if (!fscanf (randomstatusfile, "%d %d %d", ili, backward_is, backward_js)) {
        fprintf (stderr, "Error reading file %s\n", randomstatusfilename);
        return -1;
      }
      (void)fclose (randomstatusfile);
      fprintf (stderr, "reading random status for photon...\n");
    } else
      return err_out ("Error %d, could not open randomstatusfile\n", -1);
  }

#ifdef WRITERANDOMSTATUS
  /* write status every 10 photons, starting with photon 1 */
  if (photoncounter % 10 == 1 && !(*readrandomstatus)) {
    if ((randomstatusfile = fopen (randomstatusfilename, "w")) != NULL) {
      status = gsl_rng_fwrite (randomstatusfile, uvspecrng);
      if (status)
        return err_out ("Error %d returned by gsl_rng_fwrite()\n", status);
      fprintf (randomstatusfile, "%d %d %d\n", *ili, *backward_is, *backward_js);
      (void)fclose (randomstatusfile);
    } else
      return err_out ("Error %d, could not open randomstatusfile\n", -1);
  }
#endif
  *readrandomstatus = 0;

  return 0;
}
#endif


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>

#include "cosmocalc.h"

static double deltaTm[9] = {200.0, 300.0, 400.0, 600.0, 800.0, 1200.0, 1600.0, 2400.0, 3200.0};
static double ATm[9] = {0.186, 0.2, 0.212, 0.218, 0.248, 0.255, 0.260, 0.260, 0.260};
static double aTm[9] = {1.47, 1.52, 1.56, 1.61, 1.87, 2.13, 2.30, 2.53, 2.66};
static double bTm[9] = {2.57, 2.25, 2.05, 1.87, 1.59, 1.51, 1.46, 1.44, 1.41};
static double cTm[9] = {1.19, 1.27, 1.34, 1.45, 1.58, 1.80, 1.97, 2.24, 2.44};

double mass_function(double m, double a)
{
  return tinker2008_mass_function(m,a,cosmoData.delta);
}

double tinker2008_mass_function(double m, double a, double delta)
{
  double A,am,b,c,alpha;

  static int init = 1;
  static gsl_spline *ATm_spline,*aTm_spline,*bTm_spline,*cTm_spline;
  static gsl_interp_accel *ATm_acc,*aTm_acc,*bTm_acc,*cTm_acc;

  if(init)
    {
#define TINKERCODEPARMS
#ifdef TINKERCODEPARMS
      //first parameter
      ATm[0] = 1.858659e-01 ;
      ATm[1] = 1.995973e-01 ;
      ATm[2] = 2.115659e-01 ;
      ATm[3] = 2.184113e-01 ;
      ATm[4] = 2.480968e-01 ;
      ATm[5] = 2.546053e-01 ;
      ATm[6] = 2.600000e-01 ;
      ATm[7] = 2.600000e-01 ;
      ATm[8] = 2.600000e-01 ;

      //second parameter
      aTm[0] = 1.466904e+00 ;
      aTm[1] = 1.521782e+00 ;
      aTm[2] = 1.559186e+00 ;
      aTm[3] = 1.614585e+00 ;
      aTm[4] = 1.869936e+00 ;
      aTm[5] = 2.128056e+00 ;
      aTm[6] = 2.301275e+00 ;
      aTm[7] = 2.529241e+00 ;
      aTm[8] = 2.661983e+00 ;

      //third parameter
      bTm[0] = 2.571104e+00 ;
      bTm[1] = 2.254217e+00 ;
      bTm[2] = 2.048674e+00 ;
      bTm[3] = 1.869559e+00 ;
      bTm[4] = 1.588649e+00 ;
      bTm[5] = 1.507134e+00 ;
      bTm[6] = 1.464374e+00 ;
      bTm[7] = 1.436827e+00 ;
      bTm[8] = 1.405210e+00 ;

      //fourth parameter
      cTm[0] = 1.193958e+00;
      cTm[1] = 1.270316e+00;
      cTm[2] = 1.335191e+00;
      cTm[3] = 1.446266e+00;
      cTm[4] = 1.581345e+00;
      cTm[5] = 1.795050e+00;
      cTm[6] = 1.965613e+00;
      cTm[7] = 2.237466e+00;
      cTm[8] = 2.439729e+00;
#endif

      ATm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      aTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      bTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      cTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);

      gsl_spline_init(ATm_spline,deltaTm,ATm,9);
      gsl_spline_init(aTm_spline,deltaTm,aTm,9);
      gsl_spline_init(bTm_spline,deltaTm,bTm,9);
      gsl_spline_init(cTm_spline,deltaTm,cTm,9);

      ATm_acc = gsl_interp_accel_alloc();
      aTm_acc = gsl_interp_accel_alloc();
      bTm_acc = gsl_interp_accel_alloc();
      cTm_acc = gsl_interp_accel_alloc();

      init = 0;
    }

  A = gsl_spline_eval(ATm_spline,delta,ATm_acc);
#ifdef TINKERCODEPARMS
  if(delta > 1600.0)
    A = 0.26;
#endif
  am = gsl_spline_eval(aTm_spline,delta,aTm_acc);
  b = gsl_spline_eval(bTm_spline,delta,bTm_acc);
  c = gsl_spline_eval(cTm_spline,delta,cTm_acc);

  double onepz = 1.0/a;
#ifdef TINKERCODEPARMS
  if(onepz > 4.0)
    onepz = 4.0;
#endif

  A *= pow(onepz,-0.14);
  am *= pow(onepz,-0.06);
  alpha = pow(10.0,-1.0*pow(0.75/log10(delta/75.0),1.2));
  b *= pow(onepz,-1.0*alpha);

  //fprintf(stderr,"A = %f, a = %f, b = %f, c = %f\n",A,am,b,c);

  double sigma = sigmaMtophat(m,a);
  double fsigma = A*(pow(sigma/b,-1.0*am) + 1.0)*exp(-1.0*c/sigma/sigma);
  double dm = 1e-6*m;
  double dlnsiginvdm = log(sigmaMtophat(m-dm/2.0,a)/sigmaMtophat(m+dm/2.0,a))/dm;

  return fsigma*RHO_CRIT*cosmoData.OmegaM/m*dlnsiginvdm;
}

static double alphaTm[9] = {0.368, 0.363, 0.385, 0.389, 0.393, 0.365, 0.379, 0.355, 0.327};
static double betaTm[9] = {0.589, 0.585, 0.544, 0.543, 0.564, 0.623, 0.637, 0.673, 0.702};
static double gammaTm[9] = {0.864, 0.922, 0.987, 1.09, 1.20, 1.34, 1.50, 1.68, 1.81};
static double phiTm[9] = {-0.729, -0.789, -0.910, -1.05, -1.20, -1.26, -1.45, -1.50, -1.49};
static double etaTm[9] = {-0.243, -0.261, -0.261, -0.273, -0.278, -0.301, -0.301, -0.319, -0.336};

static double tinker2010_mf_norm_integ(double lnnu, void *p)
{
  double *params = (double*) p;
  double nu = exp(lnnu);
  double fnu,bnu;

  double y = log10(params[4]);
  double A = 1.0 + 0.24*y*exp(-1.0*pow(4.0/y,4.0));
  double _a = 0.44*y - 0.88;
  double B = 0.183;
  double b = 1.5;
  double C = 0.019 + 0.107*y + 0.19*exp(-1.0*pow(4.0/y,4.0));
  double c = 2.4;

  if(nu == 0.0)
    return 0.0;
  else
    {
      fnu = (1.0 + pow(params[0]*nu,-2.0*params[1]))*pow(nu,2.0*params[2])*exp(-1.0*params[3]*nu*nu/2.0);
      bnu = 1.0 - A*pow(nu,_a)/(pow(nu,_a) + pow(1.686,_a)) + B*pow(nu,b) + C*pow(nu,c);

      return bnu*fnu*nu;
    }
}

double tinker2010_mass_function(double m, double a, double delta)
{
  double z,gsigma,nu,dlnsiginvdm;
  double dm = 1e-6*m;
  double alpha,beta,gamma,phi,eta;

  static int init = 1;
  static gsl_spline *alphaTm_spline,*betaTm_spline,*gammaTm_spline,*phiTm_spline,*etaTm_spline;
  static gsl_interp_accel *alphaTm_acc,*betaTm_acc,*gammaTm_acc,*phiTm_acc,*etaTm_acc;
  //static gsl_integration_glfixed_table *gltab;
  static gsl_integration_workspace *w;
  static double aprev = -1.0;
  static double alphaprev = -1.0;

  gsl_function f;
  size_t limit = 1000;
  double epsabs,epsrel;
  double abserr;
  double numin,numax;
  //int i,Nnu,dnu;
  double params[5];

  if(init)
    {
      alphaTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      betaTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      gammaTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      phiTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);
      etaTm_spline = gsl_spline_alloc(gsl_interp_cspline,9);

      gsl_spline_init(alphaTm_spline,deltaTm,alphaTm,9);
      gsl_spline_init(betaTm_spline,deltaTm,betaTm,9);
      gsl_spline_init(gammaTm_spline,deltaTm,gammaTm,9);
      gsl_spline_init(phiTm_spline,deltaTm,phiTm,9);
      gsl_spline_init(etaTm_spline,deltaTm,etaTm,9);

      alphaTm_acc = gsl_interp_accel_alloc();
      betaTm_acc = gsl_interp_accel_alloc();
      gammaTm_acc = gsl_interp_accel_alloc();
      phiTm_acc = gsl_interp_accel_alloc();
      etaTm_acc = gsl_interp_accel_alloc();

      //gltab = gsl_integration_glfixed_table_alloc((size_t) 500);
      w = gsl_integration_workspace_alloc(limit);

      init = 0;
    }

  alpha = gsl_spline_eval(alphaTm_spline,delta,alphaTm_acc);
  beta = gsl_spline_eval(betaTm_spline,delta,betaTm_acc);
  gamma = gsl_spline_eval(gammaTm_spline,delta,gammaTm_acc);
  phi = gsl_spline_eval(phiTm_spline,delta,phiTm_acc);
  eta = gsl_spline_eval(etaTm_spline,delta,etaTm_acc);

  nu = 1.686/sigmaMtophat(m,a);
  z = 1.0/a - 1.0;
  if(z > 3.0)
    z = 3.0;

  beta = beta*pow(1.0+z,0.20);
  phi = phi*pow(1.0+z,-0.08);
  eta = eta*pow(1.0+z,0.27);
  gamma = gamma*pow(1.0+z,-0.01);

  //determine alpha by normalization condition
  if(a != aprev)
    {
      aprev = a;

      params[0] = beta;
      params[1] = phi;
      params[2] = eta;
      params[3] = gamma;
      params[4] = delta;
      f.function = &tinker2010_mf_norm_integ;
      f.params = params;
      numin = 1e-10;
      numax = 1e10;
      epsabs = 1e-6;
      epsrel = 1e-6;

      gsl_integration_qag(&f,log(numin),log(numax),epsabs,epsrel,limit,GSL_INTEG_GAUSS61,w,&alphaprev,&abserr);
      //alphaprev = gsl_integration_glfixed(&f,log(numin),log(numax),gltab);

      /*
      alphaprev = 0.0;
      Nnu = 100;

      dnu = log(numax/numin)/Nnu;
      for(i=0;i<Nnu;++i)
	alphaprev += tinker2010_mf_norm_integ(log(numin)+dnu*i,f.params)*dnu;
      */
      //fprintf(stderr,"alpha = %g, norm = %g, alpha_norm = %g\n",alpha,alphaprev,1.0/alphaprev);

      alphaprev = 1.0/alphaprev;
    }

  alpha = alphaprev;

  gsigma = alpha*(1.0 + pow(beta*nu,-2.0*phi))*pow(nu,2.0*eta)*exp(-1.0*gamma*nu*nu/2.0)*nu;
  dlnsiginvdm = log(sigmaMtophat(m-dm/2.0,a)/sigmaMtophat(m+dm/2.0,a))/dm;

  return gsigma*RHO_CRIT*cosmoData.OmegaM/m*dlnsiginvdm;
}

double bias_function(double m, double a)
{
  return tinker2010_bias(m,a,cosmoData.delta);
}

double tinker2010_bias(double m, double a, double delta)
{
  double nu;
  nu =  1.686/sigmaMtophat(m,a);

  double y = log10(delta);
  double A = 1.0 + 0.24*y*exp(-1.0*pow(4.0/y,4.0));
  double _a = 0.44*y - 0.88;
  double B = 0.183;
  double b = 1.5;
  double C = 0.019 + 0.107*y + 0.19*exp(-1.0*pow(4.0/y,4.0));
  double c = 2.4;

  return 1.0 - A*pow(nu,_a)/(pow(nu,_a) + pow(1.686,_a)) + B*pow(nu,b) + C*pow(nu,c);
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc