fortran中成为性能瓶颈的很多时候是OpenMP,而在OpenMP中多重循环最容易成为性能瓶颈，在不考虑更深层次的优化的时候，如何快速的将多重循环移植到GPU上呢？ 使用CUDA FORTRAN快速的实现Fortran的循环移植遵循三步：
第一步，设计CUDA线程结构
线程结构的设计一般采用映射的办法，将各层循环分别映射到CUDA的线程维度上。由于block的维度为3维，grid维度为3维，因此最多可以表示成9维对应9重循环（实际中一般遇不到这么深的循环），但是block的维度比较小，万一不够呢？所以我们这里只讨论最多三维的情况（维度过多可以考虑循环合并，这里不展开）。再说block的设计，block在CUDA中意味着block内的通信和同步，这就需要对程序背后的算法的深入理解，不理解算法的话就不好掌控block，这篇文章讨论的是快速移植，那些需要根据算法设计block的就不在讨论之类了。那block该多大呢？我推荐128或者256，这个尺寸的block是最能发挥线程级别并行的，设计的过大，一个SMX上的最大线程数固定，不容易接近最大的线程数，设计的过小，一个SMX上的block数也是有限制的，也容易使得SMX上的线程不够多。128或者256也是经验数字，老一辈高人传授的，说这个好
。
有了block的大小，接下来怎么设计线程的维度呢？比如说，我有二层循环，那我就分别把第一层、第二层分别对应CUDA线程的x维度、y维度。拿第一层和x维度举例，假如我的第一层循环的次数为4096，而我的block的维度为16*16=256，也就是说我的block在x维度上的大小为16，所以我需要的设计grid在x维度上的大小为4096/16，假如这个时候不是4096，而是什么4095呢，不能整除的话那就取大于该除法结果的最大整数了。同样的可以设计出y维度和z维度。那总的grid的大小就有两种算法了，第一用算出来的它的x维度*它的y维度，第二就是直接把各层循环的乘积/block的大小，这两个得到的结果其实是一样的。
第二步，分配线程任务
线程设计好了，那哪个线程干哪条活呢？在多重循环中，我们一般用控制循环的变量i,j,k在循环体内控制任务的执行，在CUDA里面，我们依然保留这三个变量用来给线程分配任务。比如说，i，j，k分别对于第一，二，三重循环，在kernel里面，i要等于线程在x维度的绝对维度，即i = (blockIdx%x-1)*blockDim%x + threadIdx%x,j要等于线程在y维度上的绝对维度，即j = (blockIdx%y-1)*blockDim%y + threadIdx%y,同样的k可以等线程在z维度上的绝对维度，这样在其他部分代码不变（甚至连后面的i,j,k都不用变）的情况下，完成了任务的重新分配。
第三步，控制线程边界
循环的步数总不可能刚刚除的“干净”，那些超了边的线程就让他歇着不干活，在上面的grid设计中有时候不能整除，所以我们要“修边”，就是用一个if语句表示该线程是否干活，只有满足条件的（即那些i,j,k比对于的循环次数小的）才去干活。
综合上面三步，基本上就完成了一个“粗糙”的移植了。至于效果，本人亲自试验过，还是不错的，只要数据的并行度高的话，基本上和精确优化差的不是太多。
第一步，设计CUDA线程结构
线程结构的设计一般采用映射的办法，将各层循环分别映射到CUDA的线程维度上。由于block的维度为3维，grid维度为3维，因此最多可以表示成9维对应9重循环（实际中一般遇不到这么深的循环），但是block的维度比较小，万一不够呢？所以我们这里只讨论最多三维的情况（维度过多可以考虑循环合并，这里不展开）。再说block的设计，block在CUDA中意味着block内的通信和同步，这就需要对程序背后的算法的深入理解，不理解算法的话就不好掌控block，这篇文章讨论的是快速移植，那些需要根据算法设计block的就不在讨论之类了。那block该多大呢？我推荐128或者256，这个尺寸的block是最能发挥线程级别并行的，设计的过大，一个SMX上的最大线程数固定，不容易接近最大的线程数，设计的过小，一个SMX上的block数也是有限制的，也容易使得SMX上的线程不够多。128或者256也是经验数字，老一辈高人传授的，说这个好
。
有了block的大小，接下来怎么设计线程的维度呢？比如说，我有二层循环，那我就分别把第一层、第二层分别对应CUDA线程的x维度、y维度。拿第一层和x维度举例，假如我的第一层循环的次数为4096，而我的block的维度为16*16=256，也就是说我的block在x维度上的大小为16，所以我需要的设计grid在x维度上的大小为4096/16，假如这个时候不是4096，而是什么4095呢，不能整除的话那就取大于该除法结果的最大整数了。同样的可以设计出y维度和z维度。那总的grid的大小就有两种算法了，第一用算出来的它的x维度*它的y维度，第二就是直接把各层循环的乘积/block的大小，这两个得到的结果其实是一样的。
第二步，分配线程任务
线程设计好了，那哪个线程干哪条活呢？在多重循环中，我们一般用控制循环的变量i,j,k在循环体内控制任务的执行，在CUDA里面，我们依然保留这三个变量用来给线程分配任务。比如说，i，j，k分别对于第一，二，三重循环，在kernel里面，i要等于线程在x维度的绝对维度，即i = (blockIdx%x-1)*blockDim%x + threadIdx%x,j要等于线程在y维度上的绝对维度，即j = (blockIdx%y-1)*blockDim%y + threadIdx%y,同样的k可以等线程在z维度上的绝对维度，这样在其他部分代码不变（甚至连后面的i,j,k都不用变）的情况下，完成了任务的重新分配。
第三步，控制线程边界
循环的步数总不可能刚刚除的“干净”，那些超了边的线程就让他歇着不干活，在上面的grid设计中有时候不能整除，所以我们要“修边”，就是用一个if语句表示该线程是否干活，只有满足条件的（即那些i,j,k比对于的循环次数小的）才去干活。
综合上面三步，基本上就完成了一个“粗糙”的移植了。至于效果，本人亲自试验过，还是不错的，只要数据的并行度高的话，基本上和精确优化差的不是太多。
原文发布时间为：2017-4-8 18:44:06
原文由：LGZ 发布，版权归属于原作者
本文来自云栖社区合作伙伴NVIDIA，了解相关信息可以关注NVIDIA官方网站
