/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's 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.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "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 Intel Corporation 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.
//
//M*/

#ifndef __OPENCV_GPU_UTILITY_HPP__
#define __OPENCV_GPU_UTILITY_HPP__

#include "internal_shared.hpp"
#include "saturate_cast.hpp"

#ifndef __CUDA_ARCH__
	#define __CUDA_ARCH__ 0
#endif

#define OPENCV_GPU_LOG_WARP_SIZE	    (5)
#define OPENCV_GPU_WARP_SIZE	        (1 << OPENCV_GPU_LOG_WARP_SIZE)
#define OPENCV_GPU_LOG_MEM_BANKS        ((__CUDA_ARCH__ >= 200) ? 5 : 4) // 32 banks on fermi, 16 on tesla
#define OPENCV_GPU_MEM_BANKS            (1 << OPENCV_GPU_LOG_MEM_BANKS)

#if defined(_WIN64) || defined(__LP64__)		
    // 64-bit register modifier for inlined asm
    #define OPENCV_GPU_ASM_PTR "l"
#else	
    // 32-bit register modifier for inlined asm
    #define OPENCV_GPU_ASM_PTR "r"
#endif

namespace cv {  namespace gpu { namespace device
{
    template <typename T> void __host__ __device__ __forceinline__ swap(T& a, T& b) 
    {
        const T temp = a;
        a = b;
        b = temp;
    }

    // warp-synchronous 32 elements reduction
    template <typename T, typename Op> __device__ __forceinline__ void warpReduce32(volatile T* data, T& partial_reduction, int tid, const Op& op)
    {
        data[tid] = partial_reduction;

        if (tid < 16)
        {
            data[tid] = partial_reduction = op(partial_reduction, data[tid + 16]);
            data[tid] = partial_reduction = op(partial_reduction, data[tid + 8 ]);
            data[tid] = partial_reduction = op(partial_reduction, data[tid + 4 ]);
            data[tid] = partial_reduction = op(partial_reduction, data[tid + 2 ]);
            data[tid] = partial_reduction = op(partial_reduction, data[tid + 1 ]);
        }
    }

    // warp-synchronous 16 elements reduction
    template <typename T, typename Op> __device__ __forceinline__ void warpReduce16(volatile T* data, T& partial_reduction, int tid, const Op& op)
    {
        data[tid] = partial_reduction;

        if (tid < 8)
        {
            data[tid] = partial_reduction = op(partial_reduction, (T)data[tid + 8 ]);
            data[tid] = partial_reduction = op(partial_reduction, (T)data[tid + 4 ]);
            data[tid] = partial_reduction = op(partial_reduction, (T)data[tid + 2 ]);
            data[tid] = partial_reduction = op(partial_reduction, (T)data[tid + 1 ]);
        }
    }

    // warp-synchronous reduction
    template <int n, typename T, typename Op> __device__ __forceinline__ void warpReduce(volatile T* data, T& partial_reduction, int tid, const Op& op)
    {
        if (tid < n)
            data[tid] = partial_reduction;

        if (n > 16)
        {
            if (tid < n - 16) 
                data[tid] = partial_reduction = op(partial_reduction, data[tid + 16]);
            if (tid < 8)
            {
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  8]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  4]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  2]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  1]);
            }
        }
        else if (n > 8)
        {
            if (tid < n - 8) 
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  8]);
            if (tid < 4)
            {
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  4]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  2]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  1]);
            }
        }
        else if (n > 4)
        {
            if (tid < n - 4) 
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  4]);
            if (tid < 2)
            {
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  2]);
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  1]);
            }
        }   
        else if (n > 2)
        {
            if (tid < n - 2) 
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  2]);
            if (tid < 2)
            {
                data[tid] = partial_reduction = op(partial_reduction, data[tid +  1]);
            }
        }      
    }

    // solve 2x2 linear system Ax=b
    template <typename T> __device__ __forceinline__ bool solve2x2(const T A[2][2], const T b[2], T x[2])
    {
        T det = A[0][0] * A[1][1] - A[1][0] * A[0][1];

        if (det != 0)
        {
            double invdet = 1.0 / det;

            x[0] = saturate_cast<T>(invdet * (b[0] * A[1][1] - b[1] * A[0][1]));

            x[1] = saturate_cast<T>(invdet * (A[0][0] * b[1] - A[1][0] * b[0]));

            return true;
        }

        return false;
    }

    // solve 3x3 linear system Ax=b
    template <typename T> __device__ __forceinline__ bool solve3x3(const T A[3][3], const T b[3], T x[3])
    {
        T det = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
              - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0])
              + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);

        if (det != 0)
        {
            double invdet = 1.0 / det;

            x[0] = saturate_cast<T>(invdet * 
                (b[0]    * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) -
                 A[0][1] * (b[1]    * A[2][2] - A[1][2] * b[2]   ) +
                 A[0][2] * (b[1]    * A[2][1] - A[1][1] * b[2]   )));

            x[1] = saturate_cast<T>(invdet * 
                (A[0][0] * (b[1]    * A[2][2] - A[1][2] * b[2]   ) -
                 b[0]    * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) +
                 A[0][2] * (A[1][0] * b[2]    - b[1]    * A[2][0])));

            x[2] = saturate_cast<T>(invdet * 
                (A[0][0] * (A[1][1] * b[2]    - b[1]    * A[2][1]) -
                 A[0][1] * (A[1][0] * b[2]    - b[1]    * A[2][0]) +
                 b[0]    * (A[1][0] * A[2][1] - A[1][1] * A[2][0])));

            return true;
        }

        return false;
    }
}}}

#endif // __OPENCV_GPU_UTILITY_HPP__
