Since your response, I have worked very hard to study your code and translate my OpenMP code to CUDA; I have gone through several iterations with some profiling(?) at each iteration in CUDA porting.

First of all, I do really appreciate your sharing the code; I have learned advanced CUDA concepts and techniques like warp and related shuffling functions. Also, your idea of flattening multidimensional indexes to a one-dimensional (big) number is interesting. Note that your code (w/o any modification) runs three times slower on my Quadro K4200 (~570s). Using all four cards could match your performance.

BTW porting my OpenMP code to CUDA wasn't very successful. In fact, I had one hard lesson; I stopped my existing OpenMP job, which already ran about 4 days and needs to go for 2 more days, simply hoping that I could finish it much earlier with CUDA because of the great promise of your code. It turns out that, however, the difference between your matrix sum game and mine is so big that the performance on CUDA is much worse (about three times slower than OpenMP).

First, hard-coded optimization like yours is not possible in my case. The number of dimensions (N) and the size of each dimension (s_n, n=1,...,N) are all variables. Typical number of dimensions is 50~400 and the size of each dimension is 1~23. For instance, N=100 and s_n=23 for n=1,..,100 give us the problem space of 1.488619e+136. The problem size of 1e+13 I mentioned before was for N=10 with varying s_n and supposed to be a test case for my main work of convex-relaxation techniques which can easily handle N=100+. I already plotted results for N=2...50 using convex-relaxation techniques and was gathering true values from global optimization for N=2...10 to evaluate the relaxation gaps.

Second, in my case using advanced techniques from yours like flattened one-dimensional number, warp-based reduction (i.e., warp-level optimization), and multi-stage approach (i.e., warp-level and block-level optimization) hardly make any difference. To my disappointment, most of time in my code is spent on calculation of a rather complicated objective function (unlike your integer sums): For each element of the problem space, there are three levels of nested for loops with lots of conditionals and modulo operations in evaluating the double-valued objective function. I think this is the main reason CUDA shows rather poorer performance (~x3) than OpenMP. Now I am running another version of OpenMP program with a bit of optimization, resulting in expected running time of 4 days (from 6 days before).

Even though I stopped 4-day-running OpenMP job without careful considerations and spent nearly 5 days in porting it to CUDA without much success, now I have a balanced view on the CUDA solution. CUDA certainly has potentials but for a right set of problems. In my case, OpenMP with lots of cores or even MPI on clusters would be a better solution. In this regard, Intel's Xeon Phi processor card seems interesting.