Fortran - design question

Author

design question

Bart Vandewoestyne

2007-03-15, 7:11 pm

I currently have the following subroutine for generating random
integers:

! Returns a random integer from the closed interval 0..max_int-1
!
subroutine random_integer(max_int, res)

integer(kind=i4b), intent(in) :: max_int
integer(kind=i4b), intent(out) :: res

real(kind=qp) :: x

call random_number(x)

res = floor(max_int*x)

end subroutine random_integer

Now it would be very convenient for me if i was able to call this
subroutine with res being a rank-1 array instead of a scalar... or even
a higher rank array...

Making the subroutine elemental here does not work because of the call to the
random_number intrinsic...

So my question is: am I forced to write a subroutine for each possible rank
that i need and then write an generic procedure interface for these
subroutines?

Or does anybody see a much better and cleaner solution?

Thanks,
Bart

--
"Share what you know. Learn what you don't."

Rich Townsend

2007-03-15, 7:11 pm

Bart Vandewoestyne wrote:
> I currently have the following subroutine for generating random
> integers:
>
> ! Returns a random integer from the closed interval 0..max_int-1
> !
> subroutine random_integer(max_int, res)
>
> integer(kind=i4b), intent(in) :: max_int
> integer(kind=i4b), intent(out) :: res
>
> real(kind=qp) :: x
>
> call random_number(x)
>
> res = floor(max_int*x)
>
> end subroutine random_integer
>
>
> Now it would be very convenient for me if i was able to call this
> subroutine with res being a rank-1 array instead of a scalar... or even
> a higher rank array...
>
> Making the subroutine elemental here does not work because of the call to the
> random_number intrinsic...
>
> So my question is: am I forced to write a subroutine for each possible rank
> that i need and then write an generic procedure interface for these
> subroutines?

I would think so; as you point out, you can never have RANDOM_NUMBER inside an
elemental routine, or indeed in the callgraph of an elemental routine.

However, it's worth recalling that RANDOM_NUMBER itself can take array
arguments. So you could just inline the above code by hand.

cheers

Rich

>
> Or does anybody see a much better and cleaner solution?
>
> Thanks,
> Bart
>

Michael Metcalf

2007-03-15, 7:11 pm

"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:S_dKh.67424$h63.346172@phobos.telenet-ops.be...
>
> Or does anybody see a much better and cleaner solution?
>
This came up on another thread recently. For rank-1 this seemed the fastest:

function scatter(nn)
integer :: nn, scatter(nn), ii, k, itemp
real :: numbers(nn)

scatter = (/ (ii, ii=0,nn-1) /)
call random_number(numbers)
do ii = 1, nn
k = int(numbers(ii)*(nn-ii+1)) + ii ! based on H. D. Knoble's
`RPERM'
itemp = scatter(ii)
scatter(ii) = scatter(k)
scatter(k) = itemp
end do

end function scatter

Regards,

Mike Metcalf

Douglas Wells

2007-03-16, 4:31 am

In article <S_dKh.67424$h63.346172@phobos.telenet-ops.be>, Bart Vandewoestyne <MyFirstName.MyLastName@telenet.be> writes:
> Now it would be very convenient for me if i was able to call this
> subroutine with res being a rank-1 array instead of a scalar... or even
> a higher rank array...
>
> Or does anybody see a much better and cleaner solution?

Well, the following version appears to meet your basic desire,
makes minimal modifications to your original procedure, seems
reasonably clean, and is accepted by at least two compilers here:

! Returns a random integer from the closed interval 0..max_int-1
!
subroutine random_integer(max_int, res)
INTEGER, PARAMETER :: i4b = KIND(1)
INTEGER, PARAMETER :: qp = KIND(1d0)

integer(kind=i4b), intent(in) :: max_int
! > integer(kind=i4b), intent(out) :: res
INTEGER(KIND=i4b), INTENT(OUT) :: res(:)

! > real(kind=qp) :: x
REAL(KIND=qp) :: x(SIZE(res))

call random_number(x)

res = floor(max_int*x)

end subroutine random_integer

It does not, of course, address your desire to handle arbitrary
higher-level ranks, but it can be straight-forwardly extended for
specific ranks, e.g.:

INTEGER(KIND=i4b), INTENT(OUT) :: res(:,:)
REAL(KIND=qp) :: x(SIZE(res,1),SIZE(res,2))

- dmw

--
.. Douglas Wells . Connection Technologies .
.. Internet: -sp9804- -at - contek.com- .

Beliavsky

2007-03-16, 4:31 am

On Mar 15, 10:54 am, Bart Vandewoestyne
<MyFirstName.MyLastN...@telenet.be> wrote:
> I currently have the following subroutine for generating random
> integers:
>
> ! Returns a random integer from the closed interval 0..max_int-1
> !
> subroutine random_integer(max_int, res)
>
> integer(kind=i4b), intent(in) :: max_int
> integer(kind=i4b), intent(out) :: res
>
> real(kind=qp) :: x
>
> call random_number(x)
>
> res = floor(max_int*x)
>
> end subroutine random_integer
>
> Now it would be very convenient for me if i was able to call this
> subroutine with res being a rank-1 array instead of a scalar... or even
> a higher rank array...
>
> Making the subroutine elemental here does not work because of the call to the
> random_number intrinsic...
>
> So my question is: am I forced to write a subroutine for each possible rank
> that i need and then write an generic procedure interface for these
> subroutines?

This is easy to do in Fortran using RESHAPE, which can be used either
inside or outside the procedures, as shown below:

module ran_mod
implicit none
private
public :: random,ran_int
interface random
module procedure random_int_vec,random_int_mat,random_int
_tensor
end interface random
contains
subroutine random_int_vec(max_int,iarr)
integer, intent(in) :: max_int
integer, intent(out) :: iarr(:)
iarr = ran_int(max_int,size(iarr))
end subroutine random_int_vec
!
subroutine random_int_mat(max_int,iarr)
integer, intent(in) :: max_int
integer, intent(out) :: iarr(:,:)
iarr = reshape(ran_int(max_int,size(iarr)),shap
e(iarr))
end subroutine random_int_mat
!
subroutine random_int_tensor(max_int,iarr)
integer, intent(in) :: max_int
integer, intent(out) :: iarr(:,:,:)
iarr = reshape(ran_int(max_int,size(iarr)),shap
e(iarr))
end subroutine random_int_tensor
!
function ran_int(max_int,n) result(ivec)
! return a vector of random integers from 1 to n, inclusive
integer, intent(in) :: max_int,n
integer :: ivec(n)
real :: x(n)
call random_number(x)
ivec = 1 + int(max_int*x)
end function ran_int
end module ran_mod

program xran
use ran_mod, only: ran_int,random
implicit none
integer, parameter :: n1 = 5, n2 = 10, n3 = 20, max_int = 3
integer ::
i,ivec(n1*n2*n3),imat(n1*n2,n3),itens(n1
,n2,n3),imethod
call random_seed()
do imethod = 1,2
if (imethod == 1) then
ivec = ran_int(max_int,size(ivec))
imat = reshape(ran_int(max_int,size(imat)) ,shape(imat))
itens = reshape(ran_int(max_int,size(itens)),sha
pe(itens))
else
call random(max_int,ivec)
call random(max_int,imat)
call random(max_int,itens)
end if
do i=1,max_int
print*,"#obs = " ,i,count(ivec==i),count(imat==i),count(i
tens==i)
end do
write (*,*)
end do
end program xran

output using g95:

#obs = 1 300 350 347
#obs = 2 342 311 314
#obs = 3 358 339 339

#obs = 1 311 359 340
#obs = 2 361 311 344
#obs = 3 328 330 316

Sponsored Links