allocate(G(N, N))

  G = 0

  ! Loop through columns of G
  do i = 1, N
     call random_number(G(i:N, i))
     G(i:N, i) = G(i:N, i) + 1
  end do

  A = matmul(G, transpose(G))

  deallocate(G)