MKL/VSL，vRngGaussian：生成高斯分布随机数

转载地址：http://blog.sina.com.cn/s/blog_73641de30102uwho.html

以产生Gauss分布随机数为例，需要的步骤有：
1. 定义type(vsl_stream_state) :: stream。
按照手册，stream的定义为“Random stream (or stream) is an abstract source of pseudo- and quasi-random sequences of uniform distribution. Users have no direct access to these sequences and operate with stream state descriptors only. A stream state descriptor, which holds state descriptive information for particular BRNG, is a necessary parameter in each routine of a distribution generator. Only the distribution generator routines operate with random streams directly. ”按照字面理解，是一种“流”，随机数“流”，uniform分布的随机数流，VSL里所有分布的随机数都是根据uniform分布随机数通过变换产生的，相当于这个“流”是产生随机数的基础。不过这里的stream并不是完整地记录下整个流（否则会非常大），只是记录下描述信息，根据这个描述信息就能完整的推断出这个随机数流（毕竟计算机里都是伪随机，不是真随机嘛）。（p.s. 不懂pseudo- & quasi- random number的区别，不过貌似准随机数要比伪随机数在对分布满足的性质上更好一点）

1. 产生stream（即给stream赋值）
   status=vslNewStream( stream, brng, seed）
   brng 是整数型，定义BRNG（Basic Random Number Generator）的方法，即产生初始uniform 随机数的方法，包括有伪随机数 和 准随机数 的不同生成方法。其取值详见手册：10. Statistical Functions -> Random Number Generators -> Basic Generators -> BRNG Parameter Definition
   Seed ，整数型，提供的种子，相当于指定这个stream的初始值。Stream方法不变，Seed也不变，则每次调用生成的随机数流都是确定的。可以通过获取系统时间获取这个seed：call system_clock(count=seed)

2. 利用确定的stream，通过变换，产生所需要的高斯分布随机数
   status = vsrnggaussian( method, stream, n, r, a, sigma )
   如 status = vsrnggaussian( VSL_METHOD_SGAUSSIAN_BOXMULLER2, stream, 2000, r, 2.0, 1.0）
   method为变换的方法，n为个数， r为接收随机数的数组， a为均值， sigma为标准差
   前四项（method, stream, n, r）为所有分布的RNG函数都有的参数，后面2个参数是高斯分布所特有的。
   VSL的函数名也有规范，v代表vsl, s代表single（同理，d代表double, i 代表integer，即返回随机数的类型；s\d对应的是连续分布随机数：如高斯分布，i对应的是离散分布随机数：如二次分布）, 后面跟rng，再跟分布名称，即 vrng

3. status = vsldeletestream( stream )： 删除stream，虽然我自己用这个函数老是报错。。

4. 其他要加的：
   include “mkl_vsl.f90”
   use mkl_vsl
   use mkl_vsl_type
   编译选项看参见MKL的link line advisor：https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor

===========================

给出一个test（注意，utility, vsl_external， inc_common.h 为我个人的程序，有些有些函数、宏为自定义，这里程序仅作示范）

Fortran语言: 高亮代码由发芽网提供
program vsl_test_gauss
use utility
use mkl_vsl_type
use mkl_vsl
use vsl_external
implicit none
character(len=*), parameter :: PROCEDURE_NAME=”vsl_test_gauss”
integer,parameter :: n=2000
real :: r(n)
type(vsl_stream_state) :: stream
integer :: seed, i
integer :: brng=VSL_BRNG_MCG31
integer :: method=VSL_METHOD_SGAUSSIAN_BOXMULLER2
real :: a = 5.0, sigma = 2.0

call system_clock(count = seed)
print*,seed
VSL_CHECK(vslnewstream(stream, brng, seed))
VSL_CHECK(vsrnggaussian( method, stream, n, r, a, sigma ))

print*, get_average(r)
print*, sqrt(get_cov(r,r))
open(11,file=”gauss_random.txt”)
do i=1, n
write(11, *) r(i)
end do
close(11)

end program vsl_test_gauss