/*
 * Copyright 1993-2007 NVIDIA Corporation.  All rights reserved.
 *
 * NOTICE TO USER:
 *
 * This source code is subject to NVIDIA ownership rights under U.S. and
 * international Copyright laws.  Users and possessors of this source code
 * are hereby granted a nonexclusive, royalty-free license to use this code
 * in individual and commercial software.
 *
 * NVIDIA MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THIS SOURCE
 * CODE FOR ANY PURPOSE.  IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR
 * IMPLIED WARRANTY OF ANY KIND.  NVIDIA DISCLAIMS ALL WARRANTIES WITH
 * REGARD TO THIS SOURCE CODE, INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY, NONINFRINGEMENT, AND FITNESS FOR A PARTICULAR PURPOSE.
 * IN NO EVENT SHALL NVIDIA BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL,
 * OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS,  WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
 * OR OTHER TORTIOUS ACTION,  ARISING OUT OF OR IN CONNECTION WITH THE USE
 * OR PERFORMANCE OF THIS SOURCE CODE.
 *
 * U.S. Government End Users.   This source code is a "commercial item" as
 * that term is defined at  48 C.F.R. 2.101 (OCT 1995), consisting  of
 * "commercial computer  software"  and "commercial computer software
 * documentation" as such terms are  used in 48 C.F.R. 12.212 (SEPT 1995)
 * and is provided to the U.S. Government only as a commercial end item.
 * Consistent with 48 C.F.R.12.212 and 48 C.F.R. 227.7202-1 through
 * 227.7202-4 (JUNE 1995), all U.S. Government End Users acquire the
 * source code with only those rights set forth herein.
 *
 * Any use of this source code in individual and commercial software must
 * include, in the user documentation and internal comments to the code,
 * the above Disclaimer and U.S. Government End Users Notice.
 */

/* Perform second step of bisection algorithm for large matrices for
 * intervals that contained after the first step more than one eigenvalue
 */

#ifndef _BISECT_KERNEL_LARGE_MULTI_H_
#define _BISECT_KERNEL_LARGE_MULTI_H_

// includes, project
#include "config.h"
#include "util.h"

// additional kernel
#include "bisect_util.cu"

////////////////////////////////////////////////////////////////////////////////
//! Perform second step of bisection algorithm for large matrices for
//! intervals that after the first step contained more than one eigenvalue
//! @param  g_d  diagonal elements of symmetric, tridiagonal matrix
//! @param  g_s  superdiagonal elements of symmetric, tridiagonal matrix
//! @param  n    matrix size
//! @param  blocks_mult  start addresses of blocks of intervals that are 
//!                      processed by one block of threads, each of the 
//!                      intervals contains more than one eigenvalue
//! @param  blocks_mult_sum  total number of eigenvalues / singleton intervals
//!                          in one block of intervals
//! @param  g_left  left limits of intervals
//! @param  g_right  right limits of intervals
//! @param  g_left_count  number of eigenvalues less than left limits
//! @param  g_right_count  number of eigenvalues less than right limits
//! @param  g_lambda  final eigenvalue
//! @param  g_pos  index of eigenvalue (in ascending order)
//! @param  precision  desired precision of eigenvalues
////////////////////////////////////////////////////////////////////////////////
__global__
void
bisectKernelLarge_MultIntervals( float* g_d, float* g_s, const unsigned int n,
                                 unsigned int* blocks_mult,
                                 unsigned int* blocks_mult_sum,
                                 float* g_left, float* g_right, 
                                 unsigned int* g_left_count,
                                 unsigned int* g_right_count,
                                 float* g_lambda, unsigned int* g_pos,
                                 float precision
                                )
{
  const unsigned int tid = threadIdx.x;

  // left and right limits of interval
  __shared__  float  s_left[2 * MAX_THREADS_BLOCK];
  __shared__  float  s_right[2 * MAX_THREADS_BLOCK];

  // number of eigenvalues smaller than interval limits
  __shared__  unsigned int  s_left_count[2 * MAX_THREADS_BLOCK];
  __shared__  unsigned int  s_right_count[2 * MAX_THREADS_BLOCK];

  // helper array for chunk compaction of second chunk
  __shared__  unsigned int  s_compaction_list[2 * MAX_THREADS_BLOCK + 1];
  // compaction list helper for exclusive scan
  unsigned int* s_compaction_list_exc = s_compaction_list + 1;

  // flag if all threads are converged
  __shared__  unsigned int  all_threads_converged;
  // number of active threads
  __shared__  unsigned int  num_threads_active;
  // number of threads to employ for compaction
  __shared__  unsigned int  num_threads_compaction;
  // flag if second chunk has to be compacted
  __shared__  unsigned int compact_second_chunk;

  // parameters of block of intervals processed by this block of threads
  __shared__  unsigned int  c_block_start;
  __shared__  unsigned int  c_block_end;
  __shared__  unsigned int  c_block_offset_output;

  // midpoint of currently active interval of the thread
  float mid = 0.0f;
  // number of eigenvalues smaller than \a mid
  unsigned int  mid_count = 0;
  // current interval parameter
  float  left;
  float  right;
  unsigned int  left_count;
  unsigned int  right_count;
  // helper for compaction, keep track which threads have a second child
  unsigned int  is_active_second = 0;

  // initialize common start conditions
  if( 0 == tid) {

    c_block_start = blocks_mult[blockIdx.x];
    c_block_end = blocks_mult[blockIdx.x + 1];
    c_block_offset_output = blocks_mult_sum[blockIdx.x];

    num_threads_active = c_block_end - c_block_start;
    s_compaction_list[0] = 0;
    num_threads_compaction = ceilPow2( num_threads_active);

    all_threads_converged = 1;
    compact_second_chunk = 0;
  }

  __syncthreads();

  // read data into shared memory
  if( tid < num_threads_active) {

    s_left[tid]  = g_left[c_block_start + tid];
    s_right[tid] = g_right[c_block_start + tid];
    s_left_count[tid]  = g_left_count[c_block_start + tid];
    s_right_count[tid] = g_right_count[c_block_start + tid];
  }

  __syncthreads();

  // do until all threads converged
  while(true) {
  //for (int iter=0; iter < 0; iter++) {

    // subdivide interval if currently active and not already converged
    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);

    __syncthreads();

    // stop if all eigenvalues have been found
    if( 1 == all_threads_converged) {
      break;
    }

    // compute number of eigenvalues smaller than mid for active and not 
    // converged intervals, use all threads for loading data from gmem and 
    // s_left and s_right as scratch space to store the data load from gmem
    // in shared memory
    mid_count = computeNumSmallerEigenvalsLarge( g_d, g_s, n,
                                                 mid, tid, num_threads_active,
                                                 s_left, s_right,
                                                 (left == right) );

    __syncthreads();

    if( tid < num_threads_active) {

      // store intervals
      if( left != right) {

        storeNonEmptyIntervals( tid, num_threads_active,
                                s_left, s_right, s_left_count, s_right_count,
                                left, mid, right, 
                                left_count, mid_count, right_count,
                                precision, compact_second_chunk,
                                s_compaction_list_exc, 
                                is_active_second ); 
      }
      else {

        storeIntervalConverged( s_left, s_right, s_left_count, s_right_count,
                                left, mid, right, 
                                left_count, mid_count, right_count, 
                                s_compaction_list_exc, compact_second_chunk,
                                num_threads_active,
                                is_active_second );

      }
    }

    __syncthreads();

    // compact second chunk of intervals if any of the threads generated 
    // two child intervals
    if( 1 == compact_second_chunk) {

      createIndicesCompaction( s_compaction_list_exc, num_threads_compaction);

      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 );
    }

    __syncthreads();

    // update state variables
    if( 0 == tid) {
      num_threads_active += s_compaction_list[num_threads_active];
      num_threads_compaction = ceilPow2( num_threads_active);

      compact_second_chunk = 0;
      all_threads_converged = 1;
    }
 
    __syncthreads();

    // clear
    s_compaction_list_exc[threadIdx.x] = 0;
    s_compaction_list_exc[threadIdx.x + blockDim.x] = 0;

    __syncthreads();

  }  // end until all threads converged

  // write data back to global memory
  if( tid < num_threads_active) {
  
    unsigned int addr = c_block_offset_output + tid;

    g_lambda[addr]  = s_left[tid];
    g_pos[addr]   = s_right_count[tid];
  }
}

#endif // #ifndef _BISECT_KERNEL_LARGE_MULTI_H_
