How to make a random function in fortran to generate the same random distribution into array?

Question

I think my code below it's not exactly give me the same random distribution.

subroutine trig_random_value()
    implicit none
    integer :: t, z, y, x
    real(real64) :: theta, r
    real(real64), parameter :: PI=4.D0*DATAN(1.D0)

    integer, dimension(12) :: date_time
    integer, dimension(12) :: seed

    call date_and_time(values=date_time)
    call random_seed
    seed = date_time(6) * date_time(7) + date_time(8)
    call random_seed(put = seed)

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                theta = rand()*2*PI
                r = 0.1*rand()
                l1(1, z, y, x) = r*cos(theta)
                l2(1, z, y, x) = r*sin(theta)

                theta = rand()*2*PI
                r = 0.1*rand()
                l1(2, z, y, x) = r*cos(theta)
                l2(2, z, y, x) = r*sin(theta)

            end do
        end do
    end do

    return
end subroutine trig_random_value

According to my code, I try to random value to l1(1,:,:,:), l1(2,:,:,:), l2(1,:,:,:) and l2(2,:,:,:) where l(t, x, y, z) is (3+1)-dimension array

Why do i use trigonometry function for my random function? because i want a circular randomization. If i plot distribution of l1(1,:,:,:) vs l2(1,:,:,:) or l1(2,:,:,:) vs l2(2,:,:,:) i will get circle shape distribution with radius 0.1

So, and why i tell you that this's not exactly give me a same distribution? because i was tried to measure a variance of them and i got

 variance_l1_t1 =   1.6670507752921395E-003
 variance_l1_t2 =   3.3313151655785292E-003
 variance_l2_t1 =   4.9965623815717321E-003
 variance_l2_t2 =   6.6641054728288360E-003

notice that (variance_l1_t2 - variance_l1_t1) = (variance_l2_t1 - variance_l1_t2) = (variance_l2_t2 - variance_l2_t1) = 0.00166

That's quite a weird result. In actually i should get almost the same variance value of l1(1,:,:,:), l1(2,:,:,:), l2(1,:,:,:) and l2(2,:,:,:) if this function if good random function. may be i did something wrong.

How to solve this problem?

Additional information from request:

real(real64) function find_variance(l)
    implicit none
    real(real64), dimension(z_size, y_size, x_size), intent(in) :: l
    integer :: z, y, x
    real(real64) :: l_avg = 0
    real(real64) :: sum_val = 0

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                l_avg = l_avg + l(z, y, x)

            end do
        end do
    end do

    l_avg = l_avg/(z_size*y_size*x_size)

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                sum_val = sum_val + (l(z , y, x) - l_avg)**2

            end do
        end do
    end do

    find_variance = sum_val/(z_size*y_size*x_size)

    return
end function find_variance

@roygvib summation of `(l(z , y, x) - l_avg)**2` divided by `N` (in this case `N` is `x_size*y_size*z_size`) — fronthem, Jun 02 '15 at 06:21
Edit: It should be really useful to put your full code for calculating variance_l1_t1, variance_l1_t2, etc in your Question, because if they have some problem, it is a different matter from random numbers. — roygvib, Jun 02 '15 at 06:43
Your 0.1 is default (single) real. Combining `real64` with `D0` and `DATAN` is strange, but should not cause an immediate problem. Btw I would not use tag Fortran 90 if you use features from Fortran 2008. — Vladimir F Героям слава, Jun 02 '15 at 07:00
I always think that this syntax of Fortran initialization is really problematic... — roygvib, Jun 02 '15 at 07:10
Any solution please purpose i just want a correct output. You can copy code above and try to check variance. — fronthem, Jun 02 '15 at 07:21
@roygvib Agreed, there were even proposals to change it, but it would brake too much code. It is too late. We have to live with the implied save. — Vladimir F Героям слава, Jun 02 '15 at 08:26
@VladimirF Actually I also got to know the historical background from http://stackoverflow.com/questions/3352741/fortran-assignment-on-declaration-and-save-attribute-gotcha. I hope compilers at least issue a warning when an initialization is done without attaching SAVE (or a SAVE statement in the scope). — roygvib, Jun 02 '15 at 08:36

score 3 · Accepted Answer · edited May 23 '17 at 10:29

3

In modern Fortran, an initialization of variables such as

real(real64) :: sum_val = 0

means that sum_val is a variable with the SAVE attribute (which is similar to static in C), which is initialized only once when the program starts. It is equivalent to

real(real64), save :: sum_val = 0

The value of the SAVEed variable is kept during the entire run and it will not be initialized to 0 again. To fix this, simply replace

real(real64) :: sum_val       !! this is a usual local variable
sum_val = 0                   !! or sum_val = real( 0, real64 )

then I guess it should be fine. Please see this page for more details. IMO this is one of the very confusing features of Fortran...

edited May 23 '17 at 10:29

Community

1
1

answered Jun 02 '15 at 07:34

roygvib

7,218
2
19
36

God... I've known this before but it's my negligence. very thanks to you @royvib – fronthem Jun 02 '15 at 07:44
`sum_val` isn't an automatic object, just a plain local one. – francescalus Jun 02 '15 at 08:40
@francescalus I checked this page http://docs.oracle.com/cd/E19957-01/805-4939/6j4m0vn7h/index.html and it also says it is not good to call `sum_val` automatic..(because F77 makes it static by default). I changed the word accordingly. Thank you :) – roygvib Jun 02 '15 at 09:07

How to make a random function in fortran to generate the same random distribution into array?

1 Answers1