/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *  * Neither the name of NVIDIA CORPORATION nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/* Determine eigenvalues for large symmetric, tridiagonal matrix. First
  step of the computation. */

#ifndef _BISECT_KERNEL_LARGE_H_
#define _BISECT_KERNEL_LARGE_H_
#include <cooperative_groups.h>

namespace cg = cooperative_groups;
// includes, project
#include "config.h"
#include "util.h"

// additional kernel
#include "bisect_util.cu"

// declaration, forward

////////////////////////////////////////////////////////////////////////////////
//! Write data to global memory
////////////////////////////////////////////////////////////////////////////////
__device__ void writeToGmem(const unsigned int tid,
                            const unsigned int tid_2,
                            const unsigned int num_threads_active,
                            const unsigned int num_blocks_mult,
                            float             *g_left_one,
                            float             *g_right_one,
                            unsigned int      *g_pos_one,
                            float             *g_left_mult,
                            float             *g_right_mult,
                            unsigned int      *g_left_count_mult,
                            unsigned int      *g_right_count_mult,
                            float             *s_left,
                            float             *s_right,
                            unsigned short    *s_left_count,
                            unsigned short    *s_right_count,
                            unsigned int      *g_blocks_mult,
                            unsigned int      *g_blocks_mult_sum,
                            unsigned short    *s_compaction_list,
                            unsigned short    *s_cl_helper,
                            unsigned int       offset_mult_lambda);

////////////////////////////////////////////////////////////////////////////////
//! Perform final stream compaction before writing out data
////////////////////////////////////////////////////////////////////////////////
__device__ void compactStreamsFinal(const unsigned int tid,
                                    const unsigned int tid_2,
                                    const unsigned int num_threads_active,
                                    unsigned int      &offset_mult_lambda,
                                    float             *s_left,
                                    float             *s_right,
                                    unsigned short    *s_left_count,
                                    unsigned short    *s_right_count,
                                    unsigned short    *s_cl_one,
                                    unsigned short    *s_cl_mult,
                                    unsigned short    *s_cl_blocking,
                                    unsigned short    *s_cl_helper,
                                    unsigned int       is_one_lambda,
                                    unsigned int       is_one_lambda_2,
                                    float             &left,
                                    float             &right,
                                    float             &left_2,
                                    float             &right_2,
                                    unsigned int      &left_count,
                                    unsigned int      &right_count,
                                    unsigned int      &left_count_2,
                                    unsigned int      &right_count_2,
                                    unsigned int       c_block_iend,
                                    unsigned int       c_sum_block,
                                    unsigned int       c_block_iend_2,
                                    unsigned int       c_sum_block_2,
                                    cg::thread_block   cta);

////////////////////////////////////////////////////////////////////////////////
//! Perform scan to compact list of block start addresses
////////////////////////////////////////////////////////////////////////////////
__device__ void scanCompactBlocksStartAddress(const unsigned int tid,
                                              const unsigned int tid_2,
                                              const unsigned int num_threads_compaction,
                                              unsigned short    *s_cl_blocking,
                                              unsigned short    *s_cl_helper,
                                              cg::thread_block   cta);

////////////////////////////////////////////////////////////////////////////////
//! Perform scan to obtain number of eigenvalues before a specific block
////////////////////////////////////////////////////////////////////////////////
__device__ void scanSumBlocks(const unsigned int tid,
                              const unsigned int tid_2,
                              const unsigned int num_threads_active,
                              const unsigned int num_threads_compaction,
                              unsigned short    *s_cl_blocking,
                              unsigned short    *s_cl_helper,
                              cg::thread_block   cta);

////////////////////////////////////////////////////////////////////////////////
//! Perform initial scan for compaction of intervals containing one and
//! multiple eigenvalues; also do initial scan to build blocks
////////////////////////////////////////////////////////////////////////////////
__device__ void scanInitial(const unsigned int tid,
                            const unsigned int tid_2,
                            const unsigned int num_threads_active,
                            const unsigned int num_threads_compaction,
                            unsigned short    *s_cl_one,
                            unsigned short    *s_cl_mult,
                            unsigned short    *s_cl_blocking,
                            unsigned short    *s_cl_helper,
                            cg::thread_block   cta);

////////////////////////////////////////////////////////////////////////////////
//! Store all non-empty intervals resulting from the subdivision of the interval
//! currently processed by the thread
//! @param  addr  address where to store
////////////////////////////////////////////////////////////////////////////////
__device__ void storeNonEmptyIntervalsLarge(unsigned int         addr,
                                            const unsigned int   num_threads_active,
                                            float               *s_left,
                                            float               *s_right,
                                            unsigned short      *s_left_count,
                                            unsigned short      *s_right_count,
                                            float                left,
                                            float                mid,
                                            float                right,
                                            const unsigned short left_count,
                                            const unsigned short mid_count,
                                            const unsigned short right_count,
                                            float                epsilon,
                                            unsigned int        &compact_second_chunk,
                                            unsigned short      *s_compaction_list,
                                            unsigned int        &is_active_second);

////////////////////////////////////////////////////////////////////////////////
//! Bisection to find eigenvalues of a real, symmetric, and tridiagonal matrix
//! @param  g_d  diagonal elements in global memory
//! @param  g_s  superdiagonal elements in global elements (stored so that the
//!              element *(g_s - 1) can be accessed an equals 0
//! @param  n   size of matrix
//! @param  lg  lower bound of input interval (e.g. Gerschgorin interval)
//! @param  ug  upper bound of input interval (e.g. Gerschgorin interval)
//! @param  lg_eig_count  number of eigenvalues that are smaller than \a lg
//! @param  lu_eig_count  number of eigenvalues that are smaller than \a lu
//! @param  epsilon  desired accuracy of eigenvalues to compute
////////////////////////////////////////////////////////////////////////////////
__global__ void bisectKernelLarge(float             *g_d,
                                  float             *g_s,
                                  const unsigned int n,
                                  const float        lg,
                                  const float        ug,
                                  const unsigned int lg_eig_count,
                                  const unsigned int ug_eig_count,
                                  float              epsilon,
                                  unsigned int      *g_num_one,
                                  unsigned int      *g_num_blocks_mult,
                                  float             *g_left_one,
                                  float             *g_right_one,
                                  unsigned int      *g_pos_one,
                                  float             *g_left_mult,
                                  float             *g_right_mult,
                                  unsigned int      *g_left_count_mult,
                                  unsigned int      *g_right_count_mult,
                                  unsigned int      *g_blocks_mult,
                                  unsigned int      *g_blocks_mult_sum)
{
    // Handle to thread block group
    cg::thread_block   cta = cg::this_thread_block();
    const unsigned int tid = threadIdx.x;

    // intervals (store left and right because the subdivision tree is in general
    // not dense
    __shared__ float s_left[2 * MAX_THREADS_BLOCK + 1];
    __shared__ float s_right[2 * MAX_THREADS_BLOCK + 1];

    // number of eigenvalues that are smaller than s_left / s_right
    // (correspondence is realized via indices)
    __shared__ unsigned short s_left_count[2 * MAX_THREADS_BLOCK + 1];
    __shared__ unsigned short s_right_count[2 * MAX_THREADS_BLOCK + 1];

    // helper for stream compaction
    __shared__ unsigned short s_compaction_list[2 * MAX_THREADS_BLOCK + 1];

    // state variables for whole block
    // if 0 then compaction of second chunk of child intervals is not necessary
    // (because all intervals had exactly one non-dead child)
    __shared__ unsigned int compact_second_chunk;
    // if 1 then all threads are converged
    __shared__ unsigned int all_threads_converged;

    // number of currently active threads
    __shared__ unsigned int num_threads_active;

    // number of threads to use for stream compaction
    __shared__ unsigned int num_threads_compaction;

    // helper for exclusive scan
    unsigned short *s_compaction_list_exc = s_compaction_list + 1;

    // variables for currently processed interval
    // left and right limit of active interval
    float        left        = 0.0f;
    float        right       = 0.0f;
    unsigned int left_count  = 0;
    unsigned int right_count = 0;
    // midpoint of active interval
    float mid = 0.0f;
    // number of eigenvalues smaller then mid
    unsigned int mid_count = 0;
    // helper for stream compaction (tracking of threads generating second child)
    unsigned int is_active_second = 0;

    // initialize lists
    s_compaction_list[tid] = 0;
    s_left[tid]            = 0;
    s_right[tid]           = 0;
    s_left_count[tid]      = 0;
    s_right_count[tid]     = 0;

    cg::sync(cta);

    // set up initial configuration
    if (0 == tid) {
        s_left[0]        = lg;
        s_right[0]       = ug;
        s_left_count[0]  = lg_eig_count;
        s_right_count[0] = ug_eig_count;

        compact_second_chunk = 0;
        num_threads_active   = 1;

        num_threads_compaction = 1;

        all_threads_converged = 1;
    }

    cg::sync(cta);

    // for all active threads read intervals from the last level
    // the number of (worst case) active threads per level l is 2^l
    while (true) {
        subdivideActiveInterval(tid,
                                s_left,
                                s_right,
                                s_left_count,
                                s_right_count,
                                num_threads_active,
                                left,
                                right,
                                left_count,
                                right_count,
                                mid,
                                all_threads_converged);

        cg::sync(cta);

        // check if done
        if (1 == all_threads_converged) {
            break;
        }

        // compute number of eigenvalues smaller than mid
        // use all threads for reading the necessary matrix data from global
        // memory
        // use s_left and s_right as scratch space for diagonal and
        // superdiagonal of matrix
        mid_count = computeNumSmallerEigenvalsLarge(
            g_d, g_s, n, mid, threadIdx.x, num_threads_active, s_left, s_right, (left == right), cta);

        cg::sync(cta);

        // store intervals
        // for all threads store the first child interval in a continuous chunk of
        // memory, and the second child interval -- if it exists -- in a second
        // chunk; it is likely that all threads reach convergence up to
        // \a epsilon at the same level; furthermore, for higher level most / all
        // threads will have only one child, storing the first child compactly will
        // (first) avoid to perform a compaction step on the first chunk, (second)
        // make it for higher levels (when all threads / intervals have
        // exactly one child)  unnecessary to perform a compaction of the second
        // chunk
        if (tid < num_threads_active) {
            if (left != right) {
                // store intervals
                storeNonEmptyIntervalsLarge(tid,
                                            num_threads_active,
                                            s_left,
                                            s_right,
                                            s_left_count,
                                            s_right_count,
                                            left,
                                            mid,
                                            right,
                                            left_count,
                                            mid_count,
                                            right_count,
                                            epsilon,
                                            compact_second_chunk,
                                            s_compaction_list_exc,
                                            is_active_second);
            }
            else {
                // re-write converged interval (has to be stored again because s_left
                // and s_right are used as scratch space for
                // computeNumSmallerEigenvalsLarge()
                s_left[tid]        = left;
                s_right[tid]       = left;
                s_left_count[tid]  = left_count;
                s_right_count[tid] = right_count;

                is_active_second = 0;
            }
        }

        // necessary so that compact_second_chunk is up-to-date
        cg::sync(cta);

        // perform compaction of chunk where second children are stored
        // scan of (num_threads_active / 2) elements, thus at most
        // (num_threads_active / 4) threads are needed
        if (compact_second_chunk > 0) {
            // create indices for compaction
            createIndicesCompaction(s_compaction_list_exc, num_threads_compaction, cta);

            compactIntervals(s_left,
                             s_right,
                             s_left_count,
                             s_right_count,
                             mid,
                             right,
                             mid_count,
                             right_count,
                             s_compaction_list,
                             num_threads_active,
                             is_active_second);
        }

        cg::sync(cta);

        // update state variables
        if (0 == tid) {
            // update number of active threads with result of reduction
            num_threads_active += s_compaction_list[num_threads_active];
            num_threads_compaction = ceilPow2(num_threads_active);

            compact_second_chunk  = 0;
            all_threads_converged = 1;
        }

        cg::sync(cta);

        if (num_threads_compaction > blockDim.x) {
            break;
        }
    }

    cg::sync(cta);

    // generate two lists of intervals; one with intervals that contain one
    // eigenvalue (or are converged), and one with intervals that need further
    // subdivision

    // perform two scans in parallel

    unsigned int left_count_2;
    unsigned int right_count_2;

    unsigned int tid_2 = tid + blockDim.x;

    // cache in per thread registers so that s_left_count and s_right_count
    // can be used for scans
    left_count  = s_left_count[tid];
    right_count = s_right_count[tid];

    // some threads have to cache data for two intervals
    if (tid_2 < num_threads_active) {
        left_count_2  = s_left_count[tid_2];
        right_count_2 = s_right_count[tid_2];
    }

    // compaction list for intervals containing one and multiple eigenvalues
    // do not affect first element for exclusive scan
    unsigned short *s_cl_one  = s_left_count + 1;
    unsigned short *s_cl_mult = s_right_count + 1;

    // compaction list for generating blocks of intervals containing multiple
    // eigenvalues
    unsigned short *s_cl_blocking = s_compaction_list_exc;
    // helper compaction list for generating blocks of intervals
    __shared__ unsigned short s_cl_helper[2 * MAX_THREADS_BLOCK + 1];

    if (0 == tid) {
        // set to 0 for exclusive scan
        s_left_count[0]  = 0;
        s_right_count[0] = 0;
    }

    cg::sync(cta);

    // flag if interval contains one or multiple eigenvalues
    unsigned int is_one_lambda   = 0;
    unsigned int is_one_lambda_2 = 0;

    // number of eigenvalues in the interval
    unsigned int multiplicity = right_count - left_count;
    is_one_lambda             = (1 == multiplicity);

    s_cl_one[tid]  = is_one_lambda;
    s_cl_mult[tid] = (!is_one_lambda);

    // (note: s_cl_blocking is non-zero only where s_cl_mult[] is non-zero)
    s_cl_blocking[tid] = (1 == is_one_lambda) ? 0 : multiplicity;
    s_cl_helper[tid]   = 0;

    if (tid_2 < num_threads_active) {
        unsigned int multiplicity = right_count_2 - left_count_2;
        is_one_lambda_2           = (1 == multiplicity);

        s_cl_one[tid_2]  = is_one_lambda_2;
        s_cl_mult[tid_2] = (!is_one_lambda_2);

        // (note: s_cl_blocking is non-zero only where s_cl_mult[] is non-zero)
        s_cl_blocking[tid_2] = (1 == is_one_lambda_2) ? 0 : multiplicity;
        s_cl_helper[tid_2]   = 0;
    }
    else if (tid_2 < (2 * MAX_THREADS_BLOCK + 1)) {
        // clear
        s_cl_blocking[tid_2] = 0;
        s_cl_helper[tid_2]   = 0;
    }

    scanInitial(
        tid, tid_2, num_threads_active, num_threads_compaction, s_cl_one, s_cl_mult, s_cl_blocking, s_cl_helper, cta);

    scanSumBlocks(tid, tid_2, num_threads_active, num_threads_compaction, s_cl_blocking, s_cl_helper, cta);

    // end down sweep of scan
    cg::sync(cta);

    unsigned int c_block_iend   = 0;
    unsigned int c_block_iend_2 = 0;
    unsigned int c_sum_block    = 0;
    unsigned int c_sum_block_2  = 0;

    // for each thread / interval that corresponds to root node of interval block
    // store start address of block and total number of eigenvalues in all blocks
    // before this block (particular thread is irrelevant, constraint is to
    // have a subset of threads so that one and only one of them is in each
    // interval)
    if (1 == s_cl_helper[tid]) {
        c_block_iend = s_cl_mult[tid] + 1;
        c_sum_block  = s_cl_blocking[tid];
    }

    if (1 == s_cl_helper[tid_2]) {
        c_block_iend_2 = s_cl_mult[tid_2] + 1;
        c_sum_block_2  = s_cl_blocking[tid_2];
    }

    scanCompactBlocksStartAddress(tid, tid_2, num_threads_compaction, s_cl_blocking, s_cl_helper, cta);

    // finished second scan for s_cl_blocking
    cg::sync(cta);

    // determine the global results
    __shared__ unsigned int num_blocks_mult;
    __shared__ unsigned int num_mult;
    __shared__ unsigned int offset_mult_lambda;

    if (0 == tid) {
        num_blocks_mult    = s_cl_blocking[num_threads_active - 1];
        offset_mult_lambda = s_cl_one[num_threads_active - 1];
        num_mult           = s_cl_mult[num_threads_active - 1];

        *g_num_one         = offset_mult_lambda;
        *g_num_blocks_mult = num_blocks_mult;
    }

    cg::sync(cta);

    float left_2, right_2;
    --s_cl_one;
    --s_cl_mult;
    --s_cl_blocking;
    compactStreamsFinal(tid,
                        tid_2,
                        num_threads_active,
                        offset_mult_lambda,
                        s_left,
                        s_right,
                        s_left_count,
                        s_right_count,
                        s_cl_one,
                        s_cl_mult,
                        s_cl_blocking,
                        s_cl_helper,
                        is_one_lambda,
                        is_one_lambda_2,
                        left,
                        right,
                        left_2,
                        right_2,
                        left_count,
                        right_count,
                        left_count_2,
                        right_count_2,
                        c_block_iend,
                        c_sum_block,
                        c_block_iend_2,
                        c_sum_block_2,
                        cta);

    cg::sync(cta);

    // final adjustment before writing out data to global memory
    if (0 == tid) {
        s_cl_blocking[num_blocks_mult] = num_mult;
        s_cl_helper[0]                 = 0;
    }

    cg::sync(cta);

    // write to global memory
    writeToGmem(tid,
                tid_2,
                num_threads_active,
                num_blocks_mult,
                g_left_one,
                g_right_one,
                g_pos_one,
                g_left_mult,
                g_right_mult,
                g_left_count_mult,
                g_right_count_mult,
                s_left,
                s_right,
                s_left_count,
                s_right_count,
                g_blocks_mult,
                g_blocks_mult_sum,
                s_compaction_list,
                s_cl_helper,
                offset_mult_lambda);
}

////////////////////////////////////////////////////////////////////////////////
//! Write data to global memory
////////////////////////////////////////////////////////////////////////////////
__device__ void writeToGmem(const unsigned int tid,
                            const unsigned int tid_2,
                            const unsigned int num_threads_active,
                            const unsigned int num_blocks_mult,
                            float             *g_left_one,
                            float             *g_right_one,
                            unsigned int      *g_pos_one,
                            float             *g_left_mult,
                            float             *g_right_mult,
                            unsigned int      *g_left_count_mult,
                            unsigned int      *g_right_count_mult,
                            float             *s_left,
                            float             *s_right,
                            unsigned short    *s_left_count,
                            unsigned short    *s_right_count,
                            unsigned int      *g_blocks_mult,
                            unsigned int      *g_blocks_mult_sum,
                            unsigned short    *s_compaction_list,
                            unsigned short    *s_cl_helper,
                            unsigned int       offset_mult_lambda)
{
    if (tid < offset_mult_lambda) {
        g_left_one[tid]  = s_left[tid];
        g_right_one[tid] = s_right[tid];
        // right count can be used to order eigenvalues without sorting
        g_pos_one[tid] = s_right_count[tid];
    }
    else {
        g_left_mult[tid - offset_mult_lambda]        = s_left[tid];
        g_right_mult[tid - offset_mult_lambda]       = s_right[tid];
        g_left_count_mult[tid - offset_mult_lambda]  = s_left_count[tid];
        g_right_count_mult[tid - offset_mult_lambda] = s_right_count[tid];
    }

    if (tid_2 < num_threads_active) {
        if (tid_2 < offset_mult_lambda) {
            g_left_one[tid_2]  = s_left[tid_2];
            g_right_one[tid_2] = s_right[tid_2];
            // right count can be used to order eigenvalues without sorting
            g_pos_one[tid_2] = s_right_count[tid_2];
        }
        else {
            g_left_mult[tid_2 - offset_mult_lambda]        = s_left[tid_2];
            g_right_mult[tid_2 - offset_mult_lambda]       = s_right[tid_2];
            g_left_count_mult[tid_2 - offset_mult_lambda]  = s_left_count[tid_2];
            g_right_count_mult[tid_2 - offset_mult_lambda] = s_right_count[tid_2];
        }

    } // end writing out data

    // note that s_cl_blocking = s_compaction_list + 1;, that is by writing out
    // s_compaction_list we write the exclusive scan result
    if (tid <= num_blocks_mult) {
        g_blocks_mult[tid]     = s_compaction_list[tid];
        g_blocks_mult_sum[tid] = s_cl_helper[tid];
    }

    if (tid_2 <= num_blocks_mult) {
        g_blocks_mult[tid_2]     = s_compaction_list[tid_2];
        g_blocks_mult_sum[tid_2] = s_cl_helper[tid_2];
    }
}

////////////////////////////////////////////////////////////////////////////////
//! Perform final stream compaction before writing data to global memory
////////////////////////////////////////////////////////////////////////////////
__device__ void compactStreamsFinal(const unsigned int tid,
                                    const unsigned int tid_2,
                                    const unsigned int num_threads_active,
                                    unsigned int      &offset_mult_lambda,
                                    float             *s_left,
                                    float             *s_right,
                                    unsigned short    *s_left_count,
                                    unsigned short    *s_right_count,
                                    unsigned short    *s_cl_one,
                                    unsigned short    *s_cl_mult,
                                    unsigned short    *s_cl_blocking,
                                    unsigned short    *s_cl_helper,
                                    unsigned int       is_one_lambda,
                                    unsigned int       is_one_lambda_2,
                                    float             &left,
                                    float             &right,
                                    float             &left_2,
                                    float             &right_2,
                                    unsigned int      &left_count,
                                    unsigned int      &right_count,
                                    unsigned int      &left_count_2,
                                    unsigned int      &right_count_2,
                                    unsigned int       c_block_iend,
                                    unsigned int       c_sum_block,
                                    unsigned int       c_block_iend_2,
                                    unsigned int       c_sum_block_2,
                                    cg::thread_block   cta)
{
    // cache data before performing compaction
    left  = s_left[tid];
    right = s_right[tid];

    if (tid_2 < num_threads_active) {
        left_2  = s_left[tid_2];
        right_2 = s_right[tid_2];
    }

    cg::sync(cta);

    // determine addresses for intervals containing multiple eigenvalues and
    // addresses for blocks of intervals
    unsigned int ptr_w            = 0;
    unsigned int ptr_w_2          = 0;
    unsigned int ptr_blocking_w   = 0;
    unsigned int ptr_blocking_w_2 = 0;

    ptr_w = (1 == is_one_lambda) ? s_cl_one[tid] : s_cl_mult[tid] + offset_mult_lambda;

    if (0 != c_block_iend) {
        ptr_blocking_w = s_cl_blocking[tid];
    }

    if (tid_2 < num_threads_active) {
        ptr_w_2 = (1 == is_one_lambda_2) ? s_cl_one[tid_2] : s_cl_mult[tid_2] + offset_mult_lambda;

        if (0 != c_block_iend_2) {
            ptr_blocking_w_2 = s_cl_blocking[tid_2];
        }
    }

    cg::sync(cta);

    // store compactly in shared mem
    s_left[ptr_w]        = left;
    s_right[ptr_w]       = right;
    s_left_count[ptr_w]  = left_count;
    s_right_count[ptr_w] = right_count;

    if (0 != c_block_iend) {
        s_cl_blocking[ptr_blocking_w + 1] = c_block_iend - 1;
        s_cl_helper[ptr_blocking_w + 1]   = c_sum_block;
    }

    if (tid_2 < num_threads_active) {
        // store compactly in shared mem
        s_left[ptr_w_2]        = left_2;
        s_right[ptr_w_2]       = right_2;
        s_left_count[ptr_w_2]  = left_count_2;
        s_right_count[ptr_w_2] = right_count_2;

        if (0 != c_block_iend_2) {
            s_cl_blocking[ptr_blocking_w_2 + 1] = c_block_iend_2 - 1;
            s_cl_helper[ptr_blocking_w_2 + 1]   = c_sum_block_2;
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
//! Compute addresses to obtain compact list of block start addresses
////////////////////////////////////////////////////////////////////////////////
__device__ void scanCompactBlocksStartAddress(const unsigned int tid,
                                              const unsigned int tid_2,
                                              const unsigned int num_threads_compaction,
                                              unsigned short    *s_cl_blocking,
                                              unsigned short    *s_cl_helper,
                                              cg::thread_block   cta)
{
    // prepare for second step of block generation: compaction of the block
    // list itself to efficiently write out these
    s_cl_blocking[tid] = s_cl_helper[tid];

    if (tid_2 < num_threads_compaction) {
        s_cl_blocking[tid_2] = s_cl_helper[tid_2];
    }

    cg::sync(cta);

    // additional scan to compact s_cl_blocking that permits to generate a
    // compact list of eigenvalue blocks each one containing about
    // MAX_THREADS_BLOCK eigenvalues (so that each of these blocks may be
    // processed by one thread block in a subsequent processing step

    unsigned int offset = 1;

    // build scan tree
    for (int d = (num_threads_compaction >> 1); d > 0; d >>= 1) {
        cg::sync(cta);

        if (tid < d) {
            unsigned int ai   = offset * (2 * tid + 1) - 1;
            unsigned int bi   = offset * (2 * tid + 2) - 1;
            s_cl_blocking[bi] = s_cl_blocking[bi] + s_cl_blocking[ai];
        }

        offset <<= 1;
    }

    // traverse down tree: first down to level 2 across
    for (int d = 2; d < num_threads_compaction; d <<= 1) {
        offset >>= 1;
        cg::sync(cta);

        //
        if (tid < (d - 1)) {
            unsigned int ai   = offset * (tid + 1) - 1;
            unsigned int bi   = ai + (offset >> 1);
            s_cl_blocking[bi] = s_cl_blocking[bi] + s_cl_blocking[ai];
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
//! Perform scan to obtain number of eigenvalues before a specific block
////////////////////////////////////////////////////////////////////////////////
__device__ void scanSumBlocks(const unsigned int tid,
                              const unsigned int tid_2,
                              const unsigned int num_threads_active,
                              const unsigned int num_threads_compaction,
                              unsigned short    *s_cl_blocking,
                              unsigned short    *s_cl_helper,
                              cg::thread_block   cta)
{
    unsigned int offset = 1;

    // first step of scan to build the sum of elements within each block
    // build up tree
    for (int d = num_threads_compaction >> 1; d > 0; d >>= 1) {
        cg::sync(cta);

        if (tid < d) {
            unsigned int ai = offset * (2 * tid + 1) - 1;
            unsigned int bi = offset * (2 * tid + 2) - 1;

            s_cl_blocking[bi] += s_cl_blocking[ai];
        }

        offset *= 2;
    }

    // first step of scan to build the sum of elements within each block
    // traverse down tree
    for (int d = 2; d < (num_threads_compaction - 1); d <<= 1) {
        offset >>= 1;
        cg::sync(cta);

        if (tid < (d - 1)) {
            unsigned int ai = offset * (tid + 1) - 1;
            unsigned int bi = ai + (offset >> 1);

            s_cl_blocking[bi] += s_cl_blocking[ai];
        }
    }

    cg::sync(cta);

    if (0 == tid) {
        // move last element of scan to last element that is valid
        // necessary because the number of threads employed for scan is a power
        // of two and not necessarily the number of active threasd
        s_cl_helper[num_threads_active - 1]   = s_cl_helper[num_threads_compaction - 1];
        s_cl_blocking[num_threads_active - 1] = s_cl_blocking[num_threads_compaction - 1];
    }
}

////////////////////////////////////////////////////////////////////////////////
//! Perform initial scan for compaction of intervals containing one and
//! multiple eigenvalues; also do initial scan to build blocks
////////////////////////////////////////////////////////////////////////////////
__device__ void scanInitial(const unsigned int tid,
                            const unsigned int tid_2,
                            const unsigned int num_threads_active,
                            const unsigned int num_threads_compaction,
                            unsigned short    *s_cl_one,
                            unsigned short    *s_cl_mult,
                            unsigned short    *s_cl_blocking,
                            unsigned short    *s_cl_helper,
                            cg::thread_block   cta)
{
    // perform scan to compactly write out the intervals containing one and
    // multiple eigenvalues
    // also generate tree for blocking of intervals containing multiple
    // eigenvalues

    unsigned int offset = 1;

    // build scan tree
    for (int d = (num_threads_compaction >> 1); d > 0; d >>= 1) {
        cg::sync(cta);

        if (tid < d) {
            unsigned int ai = offset * (2 * tid + 1);
            unsigned int bi = offset * (2 * tid + 2) - 1;

            s_cl_one[bi]  = s_cl_one[bi] + s_cl_one[ai - 1];
            s_cl_mult[bi] = s_cl_mult[bi] + s_cl_mult[ai - 1];

            // s_cl_helper is binary and zero for an internal node and 1 for a
            // root node of a tree corresponding to a block
            // s_cl_blocking contains the number of nodes in each sub-tree at each
            // iteration, the data has to be kept to compute the total number of
            // eigenvalues per block that, in turn, is needed to efficiently
            // write out data in the second step
            if ((s_cl_helper[ai - 1] != 1) || (s_cl_helper[bi] != 1)) {
                // check how many childs are non terminated
                if (s_cl_helper[ai - 1] == 1) {
                    // mark as terminated
                    s_cl_helper[bi] = 1;
                }
                else if (s_cl_helper[bi] == 1) {
                    // mark as terminated
                    s_cl_helper[ai - 1] = 1;
                }
                else // both childs are non-terminated
                {
                    unsigned int temp = s_cl_blocking[bi] + s_cl_blocking[ai - 1];

                    if (temp > MAX_THREADS_BLOCK) {
                        // the two child trees have to form separate blocks, terminate trees
                        s_cl_helper[ai - 1] = 1;
                        s_cl_helper[bi]     = 1;
                    }
                    else {
                        // build up tree by joining subtrees
                        s_cl_blocking[bi]     = temp;
                        s_cl_blocking[ai - 1] = 0;
                    }
                }
            } // end s_cl_helper update
        }

        offset <<= 1;
    }

    // traverse down tree, this only for stream compaction, not for block
    // construction
    for (int d = 2; d < num_threads_compaction; d <<= 1) {
        offset >>= 1;
        cg::sync(cta);

        //
        if (tid < (d - 1)) {
            unsigned int ai = offset * (tid + 1) - 1;
            unsigned int bi = ai + (offset >> 1);

            s_cl_one[bi]  = s_cl_one[bi] + s_cl_one[ai];
            s_cl_mult[bi] = s_cl_mult[bi] + s_cl_mult[ai];
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
//! Store all non-empty intervals resulting from the subdivision of the interval
//! currently processed by the thread
////////////////////////////////////////////////////////////////////////////////
__device__ void storeNonEmptyIntervalsLarge(unsigned int         addr,
                                            const unsigned int   num_threads_active,
                                            float               *s_left,
                                            float               *s_right,
                                            unsigned short      *s_left_count,
                                            unsigned short      *s_right_count,
                                            float                left,
                                            float                mid,
                                            float                right,
                                            const unsigned short left_count,
                                            const unsigned short mid_count,
                                            const unsigned short right_count,
                                            float                epsilon,
                                            unsigned int        &compact_second_chunk,
                                            unsigned short      *s_compaction_list,
                                            unsigned int        &is_active_second)
{
    // check if both child intervals are valid
    if ((left_count != mid_count) && (mid_count != right_count)) {
        storeInterval(addr, s_left, s_right, s_left_count, s_right_count, left, mid, left_count, mid_count, epsilon);

        is_active_second               = 1;
        s_compaction_list[threadIdx.x] = 1;
        atomicExch(&compact_second_chunk, 1);
    }
    else {
        // only one non-empty child interval

        // mark that no second child
        is_active_second               = 0;
        s_compaction_list[threadIdx.x] = 0;

        // store the one valid child interval
        if (left_count != mid_count) {
            storeInterval(
                addr, s_left, s_right, s_left_count, s_right_count, left, mid, left_count, mid_count, epsilon);
        }
        else {
            storeInterval(
                addr, s_left, s_right, s_left_count, s_right_count, mid, right, mid_count, right_count, epsilon);
        }
    }
}

#endif // #ifndef _BISECT_KERNEL_LARGE_H_
