safegaurd Box-Muller normal random number generation against u=0.0, a…

…lso improves performance ~10x

src/nf/nf_conv2d_layer_submodule.f90

@@ -1,7 +1,7 @@
 submodule(nf_conv2d_layer) nf_conv2d_layer_submodule
 
   use nf_activation, only: activation_function
-  use nf_random, only: randn
+  use nf_random, only: randn_vec
 
   implicit none
 
@@ -40,9 +40,8 @@ module subroutine init(self, input_shape)
                            self % kernel_size, self % kernel_size))
 
     ! Initialize the kernel with random values with a normal distribution.
-    self % kernel = randn(self % filters, self % channels, &
-                          self % kernel_size, self % kernel_size) &
-                  / self % kernel_size**2 !TODO kernel_width * kernel_height
+    call randn_vec(self % kernel, self % filters * self % channels * self % kernel_size**2)
+    self % kernel = self % kernel / self % kernel_size**2
 
     allocate(self % biases(self % filters))
     self % biases = 0

src/nf/nf_dense_layer_submodule.f90

@@ -2,7 +2,7 @@
 
   use nf_activation, only: activation_function
   use nf_base_layer, only: base_layer
-  use nf_random, only: randn
+  use nf_random, only: randn_vec
 
   implicit none
 
@@ -114,8 +114,8 @@ module subroutine init(self, input_shape)
     ! Weights are a 2-d array of shape previous layer size
     ! times this layer size.
     allocate(self % weights(self % input_size, self % output_size))
-    self % weights = randn(self % input_size, self % output_size) &
-                   / self % input_size
+    call randn_vec(self % weights, self % input_size * self % output_size)
+    self % weights = self % weights / self % input_size
 
     ! Broadcast weights to all other images, if any.
     call co_broadcast(self % weights, 1)

src/nf/nf_random.f90

@@ -6,7 +6,7 @@ module nf_random
   implicit none
 
   private
-  public :: randn
+  public :: randn, randn_vec
 
   interface randn
 
@@ -36,4 +36,26 @@ end function randn_4d
 
   end interface randn
 
+  real, parameter :: pi = 4 * atan(1.d0)
+
+  contains
+
+  impure subroutine randn_vec(x, n)
+    real, intent(out) :: x(*)
+    integer, intent(in), value :: n
+    real :: u((n+1)/2), v((n+1)/2), r((n+1)/2)
+    integer :: i, nu
+    call random_number(u)
+    nu = size(u)
+    do i=1,nu
+      do while (u(i) == 0.0)
+        call random_number(u(i))
+      end do
+    end do
+    call random_number(v)
+    r = sqrt(-2.0*log(u))
+    x(1:nu) = r*sin(2.0*pi*v)
+    if (n > (nu+1)) x(nu+1:n) = r(1:(n-nu))*cos(2.0*pi*v(1:(n-nu)))
+  end subroutine randn_vec
+
 end module nf_random

src/nf/nf_random_submodule.f90

@@ -1,8 +1,6 @@
 submodule(nf_random) nf_random_submodule
   implicit none
 
-  real, parameter :: pi = 4 * atan(1.d0)
-
 contains
 
   module function randn_1d(i) result(x)