Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

populating arrays

47 views
Skip to first unread message

George

unread,
Dec 12, 2008, 10:43:01 PM12/12/08
to


I've been reading a lot of fortran reference lately and have a grab bag of
questions that are all over the map.

If I want to populate this array as in pg 255 of the Adams text, how do I
do it with the least effort? Now, it's an obvious shortcoming:

implicit none

integer a(3,3)

a = [0, 3, 6, 1, 4, 7, 2, 5, 8]

print *, a

endprogram

!g95 -Wall setback1.f03 -o e.exe

C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe

In file setback1.f03:5

a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
1
Error: Incompatible ranks in assignment at (1) (2/1)

C:\MinGW\source>

Some might think the obvious answer is to use reshape. In light of today's
reading, I can certainly believe there to be a slicker way of doing things.

Let me ask another question. Pg. 252

FORALL (I=1:10, J=1:10, A(I)=0.0 .AND. B(J)>0.0
...
end forall

If this is correct triplet notation, when is the third a scalar-logical
expression as opposed to a stride, which could easily equal 2?

Oh, heck, let's go for a hat trick.

q3) With allocatable, you can, at runtime, create a 42 x 26 matrix as
opposed to a 5 x 8. Can you, portably at runtime, if heads, create a 5 x 8
x 3, if tails, a 7 x 4?

--
George

To those of you who received honours, awards and distinctions, I say well
done. And to the C students, I say you, too, can be president of the United
States.
George W. Bush

Picture of the Day http://apod.nasa.gov/apod/

Jason Blevins

unread,
Dec 13, 2008, 2:02:50 PM12/13/08
to
On Dec 12, 2008, at 10:43 PM, geo...@example.invalid wrote:
> If I want to populate this array as in pg 255 of the Adams text, how do I
> do it with the least effort? Now, it's an obvious shortcoming:
>
> implicit none
>
> integer a(3,3)
>
> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
>
> print *, a
>
> endprogram
>
> !g95 -Wall setback1.f03 -o e.exe
>
> C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe
>
> In file setback1.f03:5
>
> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
> 1
> Error: Incompatible ranks in assignment at (1) (2/1)
>
> C:\MinGW\source>
>
> Some might think the obvious answer is to use reshape. In light of today's
> reading, I can certainly believe there to be a slicker way of doing things.

According to a message [1] on the GFortran list, array constructors
should always return rank-1 arrays. Given that, it looks like reshape
may be the best solution. (I hope I'm wrong!)

[1]: http://gcc.gnu.org/ml/fortran/2005-03/msg00149.html

> q3) With allocatable, you can, at runtime, create a 42 x 26 matrix as
> opposed to a 5 x 8. Can you, portably at runtime, if heads, create a 5 x 8
> x 3, if tails, a 7 x 4?

Is this what you mean?

program allocate_coin_toss
implicit none
real, dimension(:,:), allocatable :: a
real :: u

call random_number(u)

if (u < 0.5) then
allocate(a(5,8))
else
allocate(a(7,4))
end if

print *, shape(a)

deallocate(a)
end program allocate_coin_toss

It's not much of a coin toss without re-seeding each time though :)

--
Jason Blevins
Ph.D. Candidate, Department of Economics, Duke University
http://jblevins.org/

George

unread,
Dec 13, 2008, 3:29:30 PM12/13/08
to
On Sat, 13 Dec 2008 14:02:50 -0500, Jason Blevins wrote:

> On Dec 12, 2008, at 10:43 PM, geo...@example.invalid wrote:
>> If I want to populate this array as in pg 255 of the Adams text, how do I
>> do it with the least effort? Now, it's an obvious shortcoming:
>>
>> implicit none
>>
>> integer a(3,3)
>>
>> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
>>
>> print *, a
>>
>> endprogram
>>
>> !g95 -Wall setback1.f03 -o e.exe
>>
>> C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe
>>
>> In file setback1.f03:5
>>
>> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
>> 1
>> Error: Incompatible ranks in assignment at (1) (2/1)
>>
>> C:\MinGW\source>
>>
>> Some might think the obvious answer is to use reshape. In light of today's
>> reading, I can certainly believe there to be a slicker way of doing things.
>
> According to a message [1] on the GFortran list, array constructors
> should always return rank-1 arrays. Given that, it looks like reshape
> may be the best solution. (I hope I'm wrong!)
>
> [1]: http://gcc.gnu.org/ml/fortran/2005-03/msg00149.html

Thanks, Jason, that link was helpful. Actually, the messages following up
were more revealing.

implicit none

integer a(3,3)
integer b(3,3)

a = reshape ([0, 3, 6, 1, 4, 7, 2, 5, 8], [3,3])
b = reshape ([0, 3, 6, 1, 4, 7, 2, 5, 8], shape(B))

print *, a
print *, shape(a)

print *, b
print *, shape(b)

a = matmul(b, b)
print *, a

endprogram

!g95 -Wall setback1.f03 -o e.exe

C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe

C:\MinGW\source>e
0 3 6 1 4 7 2 5 8
3 3
0 3 6 1 4 7 2 5 8
3 3
15 42 69 18 54 90 21 66 111

C:\MinGW\source>

I wonder why the print statement doesn't print them out like you'd expect?


>
>> q3) With allocatable, you can, at runtime, create a 42 x 26 matrix as
>> opposed to a 5 x 8. Can you, portably at runtime, if heads, create a 5 x 8
>> x 3, if tails, a 7 x 4?
>
> Is this what you mean?
>
> program allocate_coin_toss
> implicit none
> real, dimension(:,:), allocatable :: a
> real :: u
>
> call random_number(u)
>
> if (u < 0.5) then
> allocate(a(5,8))
> else
> allocate(a(7,4))
> end if
>
> print *, shape(a)
>
> deallocate(a)
> end program allocate_coin_toss
>
> It's not much of a coin toss without re-seeding each time though :)

I thought I would address the seeding question with the following:

program allocate_coin_toss
implicit none
real, dimension(:,:), allocatable :: a
real :: u

call init_random_seed()

call random_number(u)
print *, " u is ", u



if (u < 0.5) then
allocate(a(5,8))
else
allocate(a(7,4))
end if

print *, shape(a)

deallocate(a)

contains
SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
print *, "n=", n
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)
print *, "clock=", clock

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE

end program allocate_coin_toss

! g95 -Wall duke1.f03 -o h.exe

C:\MinGW\source>g95 -Wall duke1.f03 -o h.exe

C:\MinGW\source>h
n= 4
clock= 0
seed= 0 37 74 111
u is 0.8157187
7 4

C:\MinGW\source>h
n= 4
clock= 0
seed= 0 37 74 111
u is 0.8157187
7 4

C:\MinGW\source>

Since the clock is always zero, g95 isn't giving me the behavior I want.
Does this work for anyone else?

Q3 was a slightly differing question as it addresses the feasibility of
producing arrays of differing rank dynamically.

Go Blue Devils.
--
George

We can't allow the world's worst leaders to blackmail, threaten, hold
freedom-loving nations hostage with the world's worst weapons.

Jason Blevins

unread,
Dec 13, 2008, 4:38:24 PM12/13/08
to
On Dec 13, 2008, at 3:29 PM, geo...@example.invalid wrote:
> implicit none
>
> integer a(3,3)
> integer b(3,3)
>
> a = reshape ([0, 3, 6, 1, 4, 7, 2, 5, 8], [3,3])
> b = reshape ([0, 3, 6, 1, 4, 7, 2, 5, 8], shape(B))
>
> print *, a
> print *, shape(a)
>
> print *, b
> print *, shape(b)
>
> a = matmul(b, b)
> print *, a
>
> endprogram
>
> !g95 -Wall setback1.f03 -o e.exe
>
> C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe
>
> C:\MinGW\source>e
> 0 3 6 1 4 7 2 5 8
> 3 3
> 0 3 6 1 4 7 2 5 8
> 3 3
> 15 42 69 18 54 90 21 66 111
>
> C:\MinGW\source>
>
> I wonder why the print statement doesn't print them out like you'd expect?

As far as I know, print and write don't perform any special handling
based on the rank of the arrays you give them. I don't know of better
ways than simply looping over the rows using either an explicit or
implied do loop:

implicit none

integer a(3,3), i



a = reshape ([0, 3, 6, 1, 4, 7, 2, 5, 8], [3,3])

do i = 1, size(a,1)
print *, a(i,:)
end do

print '(3i4)', ( a(i,:), i = 1, 3 )

endprogram

It works for me with GFortran:

% gfortran -o time time.f90
% ./time
n= 8
clock= 843121476
seed= 843121476 843121513 843121550 843121587 843121624 843121661 843121698 843121735
u is 0.77183896
7 4
% ./time
n= 8
clock= 843122101
seed= 843122101 843122138 843122175 843122212 843122249 843122286 843122323 843122360
u is 0.40753603
5 8

> Q3 was a slightly differing question as it addresses the feasibility of
> producing arrays of differing rank dynamically.

Ahh, now I see. I think I was confused by the line-wrapping. I
didn't see the trailing "x 3" part in the first array! I don't know
of any way to allocate arrays of different ranks. I think the rank
has to be known at compile time.

Best

Dick Hendrickson

unread,
Dec 15, 2008, 2:10:49 PM12/15/08
to
George wrote:
>
>
> I've been reading a lot of fortran reference lately and have a grab bag of
> questions that are all over the map.
>
> If I want to populate this array as in pg 255 of the Adams text, how do I
> do it with the least effort? Now, it's an obvious shortcoming:
>
> implicit none
>
> integer a(3,3)
>
> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
>
> print *, a
>
> endprogram
>
> !g95 -Wall setback1.f03 -o e.exe
>
> C:\MinGW\source>g95 -Wall setback1.f03 -o e.exe
>
> In file setback1.f03:5
>
> a = [0, 3, 6, 1, 4, 7, 2, 5, 8]
> 1
> Error: Incompatible ranks in assignment at (1) (2/1)

I think this is trivial with a defined assignment. You merely
need to provide assignment subroutines for each type and rank
that "a" and the right hand side could have.

>
> C:\MinGW\source>
>
> Some might think the obvious answer is to use reshape. In light of today's
> reading, I can certainly believe there to be a slicker way of doing things.
>
> Let me ask another question. Pg. 252
>
> FORALL (I=1:10, J=1:10, A(I)=0.0 .AND. B(J)>0.0
> ...
> end forall
>
> If this is correct triplet notation, when is the third a scalar-logical
> expression as opposed to a stride, which could easily equal 2?

Looks almost correct to me. You need == instead of = in the logical
expression and you're missing the closing ).

The third thing is a scalar-logical expression because it just is ;).
A FORALL header has a sequence of integer names, an = sign, and 2 (or 3)
integer expressions separated by colon(s). There can be an arbitrary
long list of these, separated by commas. Eventually, the parser comes
to a closing expression or a scalar logical expression (or a syntax
error). In your example, it can't be a "stride" because it's separated
from the final value by a colon and it's not an integer expression.

The logical expression can't (not even the hard way) be "2". It can
only be TRUE or FALSE.

Dick Hendrickson

>
> Oh, heck, let's go for a hat trick.
>
> q3) With allocatable, you can, at runtime, create a 42 x 26 matrix as
> opposed to a 5 x 8. Can you, portably at runtime, if heads, create a 5 x 8
> x 3, if tails, a 7 x 4?
>

You could probably also do it with something like

a = reshape (value_list, shape = merge( [5,8], [7,4], heads ))
where a is a 2d allocatable array. The magic deallocate/allocate
should give a either a [5,8] or [7,4] shape, depending on
whether or not heads is true.

Dick Hendrickson

George

unread,
Dec 17, 2008, 9:05:23 PM12/17/08
to

I broke my bug-reporting cherry today, as I heard from elliot that he gets
zeros from the systemclock call on g95 as well. Elsethread, I've been
trying to work out what the taxonomy of a bug is, and I think this might
fall into a category where the standard allows an implementation to return
zero on a systemclock call, but it completely undermines its use to be
identical from one call to the next here.

Richard has talked about this before, and I've gotten very much the
impression that standardizing behavior for randomseed calls was a moment
that fortran missed.

BTW, Andy's got a slick captcha screen. The part I like best is that I
didn't have to ingest LSD to read it correctly:

http://i429.photobucket.com/albums/qq15/george196884/fortran38.jpg
--
George

Everywhere that freedom stirs, let tyrants fear.

George

unread,
Dec 18, 2008, 3:04:00 AM12/18/08
to
On Wed, 17 Dec 2008 19:05:23 -0700, George wrote:

> I broke my bug-reporting cherry today, as I heard from elliot that he gets
> zeros from the systemclock call on g95 as well. Elsethread, I've been
> trying to work out what the taxonomy of a bug is, and I think this might
> fall into a category where the standard allows an implementation to return
> zero on a systemclock call, but it completely undermines its use to be
> identical from one call to the next here.
>
> Richard has talked about this before, and I've gotten very much the
> impression that standardizing behavior for randomseed calls was a moment
> that fortran missed.

Above, I conflated things that have nothing to do with each other.

Andy responded with this:

It isn't a bug. SYSTEM_CLOCK() is intended as an interval timer. To
that end, I've rigged it so that it always returns zero the first time
that it is called. This gives you the maximum amount of time before you
have to worry about the ticks rolling over. In g95, the first call starts
the clock ticking.

If you want to seed the random number generator with the time, call
DATE_AND_TIME() with the value= parameter, and use the
minutes/seconds/milliseconds values. 60 gives you almost 6 bits, 1000
gives you almost 10. So that's 6+6+10=22, call it 20 bits, or a million
possible values. If you need more, hours/day/month/year should be plenty.

! end response

Looks like I'll need to work this up.
--
George

You teach a child to read, and he or her will be able to pass a literacy
test.

George

unread,
Dec 19, 2008, 12:10:36 AM12/19/08
to
On Thu, 18 Dec 2008 01:04:00 -0700, George wrote:

> Andy responded with this:
>
> It isn't a bug. SYSTEM_CLOCK() is intended as an interval timer. To
> that end, I've rigged it so that it always returns zero the first time
> that it is called. This gives you the maximum amount of time before you
> have to worry about the ticks rolling over. In g95, the first call starts
> the clock ticking.
>
> If you want to seed the random number generator with the time, call
> DATE_AND_TIME() with the value= parameter, and use the
> minutes/seconds/milliseconds values. 60 gives you almost 6 bits, 1000
> gives you almost 10. So that's 6+6+10=22, call it 20 bits, or a million
> possible values. If you need more, hours/day/month/year should be plenty.

I'm a little out of my depth here. I don't understand Andy's last
sentence, and I get partial results.

program allocate_coin_toss
!implicit none


real, dimension(:,:), allocatable :: a
real :: u

call init_g95_seed()



call random_number(u)
print *, " u is ", u

if (u < 0.5) then
allocate(a(5,8))
else
allocate(a(7,4))
end if

print *, shape(a)

deallocate(a)

contains
SUBROUTINE init_g95_seed()
INTEGER :: i, n


INTEGER, DIMENSION(:), ALLOCATABLE :: seed

INTEGER date_time(8)
CHARACTER(LEN=10) big_ben(3)

CALL DATE_AND_TIME(big_ben(1), big_ben(2), big_ben(3), date_time)
PRINT *,'date_time array values:'
PRINT *,'year=',date_time(1)
PRINT *,'month_of_year=',date_time(2)
PRINT *,'day_of_month=',date_time(3)
PRINT *,'time difference in minutes=',date_time(4)
PRINT *,'hour of day=',date_time(5)
PRINT *,'minutes of hour=',date_time(6)
PRINT *,'seconds of minute=',date_time(7)
PRINT *,'milliseconds of second=',date_time(8)
PRINT *, 'DATE=',big_ben(1)
PRINT *, 'TIME=',big_ben(2)
PRINT *, 'ZONE=',big_ben(3)

CALL RANDOM_SEED(size = n)
print *, "n=", n


ALLOCATE(seed(n))

do i = 1, n
seed(i) = date_time(9-i)
end do

print *, seed


! seed = date_time(8) + 37 * (/ (i - 1, i = 1, n) /)


CALL RANDOM_SEED(PUT = seed)
print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE

end program allocate_coin_toss

! g95 -Wall duke2.f03 -o h.exe

C:\MinGW\source>g95 -Wall duke2.f03 -o h.exe

C:\MinGW\source>h
date_time array values:
year= 2008
month_of_year= 12
day_of_month= 18
time difference in minutes= -420
hour of day= 22
minutes of hour= 1
seconds of minute= 32
milliseconds of second= 625
DATE=20081218
TIME=220132.625
ZONE=-0700
n= 4
625 32 1 22
seed= 625 32 1 22
u is 0.8378754
7 4

C:\MinGW\source>h
date_time array values:
year= 2008
month_of_year= 12
day_of_month= 18
time difference in minutes= -420
hour of day= 22
minutes of hour= 1
seconds of minute= 41
milliseconds of second= 78
DATE=20081218
TIME=220141.078
ZONE=-0700
n= 4
78 41 1 22
seed= 78 41 1 22
u is 0.812684
7 4

C:\MinGW\source>h
date_time array values:
year= 2008
month_of_year= 12
day_of_month= 18
time difference in minutes= -420
hour of day= 22
minutes of hour= 2
seconds of minute= 7
milliseconds of second= 968
DATE=20081218
TIME=220207.968
ZONE=-0700
n= 4
968 7 2 22
seed= 968 7 2 22
u is 0.83060235
7 4

C:\MinGW\source>

I thought I would vary the first input as much as I could, yet I only get
numbers near .82 . Fishing for tips.
--
George

America is a Nation with a mission - and that mission comes from our most
basic beliefs. We have no desire to dominate, no ambitions of empire. Our
aim is a democratic peace - a peace founded upon the dignity and rights of
every man and woman.

George

unread,
Dec 19, 2008, 12:18:26 AM12/19/08
to
On Mon, 15 Dec 2008 19:10:49 GMT, Dick Hendrickson wrote:

> George wrote:

> I think this is trivial with a defined assignment. You merely
> need to provide assignment subroutines for each type and rank
> that "a" and the right hand side could have.

I'll take this up later.

>> FORALL (I=1:10, J=1:10, A(I)=0.0 .AND. B(J)>0.0
>> ...
>> end forall
>>
>> If this is correct triplet notation, when is the third a scalar-logical
>> expression as opposed to a stride, which could easily equal 2?
>
> Looks almost correct to me. You need == instead of = in the logical
> expression and you're missing the closing ).

Actually it's /= .


>
> The third thing is a scalar-logical expression because it just is ;).
> A FORALL header has a sequence of integer names, an = sign, and 2 (or 3)
> integer expressions separated by colon(s). There can be an arbitrary
> long list of these, separated by commas. Eventually, the parser comes
> to a closing expression or a scalar logical expression (or a syntax
> error). In your example, it can't be a "stride" because it's separated
> from the final value by a colon and it's not an integer expression.
>
> The logical expression can't (not even the hard way) be "2". It can
> only be TRUE or FALSE.

ok


>> q3) With allocatable, you can, at runtime, create a 42 x 26 matrix as
>> opposed to a 5 x 8. Can you, portably at runtime, if heads, create a 5 x 8
>> x 3, if tails, a 7 x 4?
>>
>
> You could probably also do it with something like
>
> a = reshape (value_list, shape = merge( [5,8], [7,4], heads ))
> where a is a 2d allocatable array. The magic deallocate/allocate
> should give a either a [5,8] or [7,4] shape, depending on
> whether or not heads is true.
>
> Dick Hendrickson

Thanks for your reply, Dick, I'll take this up after I hit page 713 of your
text. I've got more fortran to learn along the way than do others; I'm
only beginning chapter 9.
--
George

There's no bigger task than protecting the homeland of our country.

user1

unread,
Dec 19, 2008, 1:03:23 PM12/19/08
to

Why not just use a minor variation of what you had originally, with SYSTEM_CLOCK
replaced with DATE_AND_TIME ??

SUBROUTINE init_random_seed()
INTEGER :: i, n, clock, ival(8)


INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))
CALL DATE_AND_TIME(values=ival)
clock = 3600*ival(5) + 60*ival(6) + ival(7)
clock = 1000*clock + ival(8)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE

George

unread,
Dec 20, 2008, 11:08:55 PM12/20/08
to

Thanks for your post, user1. I don't think we're quite there yet, but this
is a big step.


implicit none

real :: u, sum1
integer :: i, j

sum1 = 0.0
do i = 1, 20
call init_random_seed()

call random_number(u)
print *, " u is ", u

sum1 = sum1 + u

! I need a pause here of 30 milliseconds
do j = 1, 1000000
call random_number(u)
end do

end do

print *, "sum over 20 trials is ", sum1

contains


SUBROUTINE init_random_seed()
INTEGER :: i, n, clock, ival(8)
INTEGER, DIMENSION(:), ALLOCATABLE :: seed
CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))
CALL DATE_AND_TIME(values=ival)
clock = 3600*ival(5) + 60*ival(6) + ival(7)
clock = 1000*clock + ival(8)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)

print *, seed


CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE

endprogram


! g95 -Wall duke3.f03 -o h.exe


C:\MinGW\source> g95 -Wall duke3.f03 -o h.exe

C:\MinGW\source>h
75261828 75261865 75261902 75261939
u is 0.4071415
75261843 75261880 75261917 75261954
u is 0.40404147
75261875 75261912 75261949 75261986
u is 0.4026565
75261906 75261943 75261980 75262017
u is 0.4047169
75261937 75261974 75262011 75262048
u is 0.4053532
75261968 75262005 75262042 75262079
u is 0.46627986
75262000 75262037 75262074 75262111
u is 0.4654197
75262015 75262052 75262089 75262126
u is 0.46570802
75262046 75262083 75262120 75262157
u is 0.46641546
75262078 75262115 75262152 75262189
u is 0.46503907
75262109 75262146 75262183 75262220
u is 0.46179205
75262140 75262177 75262214 75262251
u is 0.46218264
75262156 75262193 75262230 75262267
u is 0.46388954
75262187 75262224 75262261 75262298
u is 0.4631936
75262218 75262255 75262292 75262329
u is 0.4556914
75262250 75262287 75262324 75262361
u is 0.45635337
75262281 75262318 75262355 75262392
u is 0.45500475
75262312 75262349 75262386 75262423
u is 0.45462817
75262328 75262365 75262402 75262439
u is 0.45458508
75262359 75262396 75262433 75262470
u is 0.47835493
sum over 20 trials is 8.9584465

C:\MinGW\source>

I've concocted a program that shows lack of randomness above. Normally,
you would think to call init_any_seed once, but I do so serially here and
with measurable time differences.

I can only speculate how this runs on other machines. Please comment if
you can. When I had longer trial runs, the output for this seemed to have
a "period".

Andy sent me some source to try out. I'll see how that goes. One question
here:

Is there a portable way to waste 30 milliseconds?
--
George

Do I think faith will be an important part of being a good president? Yes,
I do.

user1

unread,
Dec 21, 2008, 9:36:15 AM12/21/08
to

Why do you want to reseed the random number generation on each loop iteration ?
I moved it outside the loop, and the values of u produce on two successive runs
are listed below.

> implicit none
>
> real :: u, sum1
> integer :: i, j
>
> sum1 = 0.0

> call init_random_seed()

>
> do i = 1, 20

> call random_number(u)
> print *, " u is ", u
> sum1 = sum1 + u
>
> ! I need a pause here of 30 milliseconds
> do j = 1, 1000000
> call random_number(u)
> end do
>
> end do
>

------------------------------------------------------------------------------------------------------------
34479890 34479927 34479964 34480001
u is 0.62775445
u is 0.08166063
u is 0.44228584
u is 0.031528175
u is 0.65559906
u is 0.24459541
u is 0.89109623
u is 0.6606597
u is 0.22105855
u is 0.13889772
u is 0.1406675
u is 0.29403526
u is 0.8427418
u is 0.28213042
u is 0.73694146
u is 0.013806403
u is 0.74734735
u is 0.40877485
u is 0.051763177
u is 0.18899906
sum over 20 trials is 7.702343

------------------------------------------------------------------------------------------------------------------


34483593 34483630 34483667 34483704
u is 0.8823633
u is 0.06700063
u is 0.82393557
u is 0.14112788
u is 0.66519684
u is 0.034879684
u is 0.5340405
u is 0.10210335
u is 0.41274947
u is 0.77753377
u is 0.18015385
u is 0.25273865
u is 0.43116927
u is 0.5845781
u is 0.90469617
u is 0.9799248
u is 0.37720376
u is 0.21182173
u is 0.53122145
u is 0.27423805
sum over 20 trials is 9.168675


user1

unread,
Dec 21, 2008, 10:40:17 AM12/21/08
to
George wrote:

>
> Is there a portable way to waste 30 milliseconds?

I suppose you could build a time delay around repeated calls to DATE_AND_TIME.
It seems that the clock updates every 15.6 milliseconds, so you would not get
exactly 30 milliseconds. (probably 31,32)

Possible demo code follows.


implicit none
!
! demo time delay
!
integer timer
print *,timer()
call delay(30)
print *,timer()
end


integer function timer()
!
! returns milliseconds past midnight
!
implicit none
integer seconds,milliseconds
integer iv(8)
call DATE_AND_TIME(values=iv)
seconds = 3600*iv(5) + 60*iv(6) + iv(7)
milliseconds = 1000*seconds+ iv(8)
timer=milliseconds
return
end


subroutine delay (dt)
!
! delay for dt milliseconds
!
implicit none
integer dt,timer,t1,t2
t1=timer()
t2=t1
do while (t2.lt.t1+dt)
t2=timer()
enddo
return
end


Ron Shepard

unread,
Dec 21, 2008, 12:14:29 PM12/21/08
to
> Is there a portable way to waste 30 milliseconds?

I know of no standard way in fortran to pause a program for a
specified time. A nonstandard way might be something like

call system( 'sleep 10' )

to pause for 10 seconds, for example; this does not work for
fractions of a second, which is the original question. You can
write a C function that accesses the standard usleep() function,
which uses units of microseconds, or the nanosleep() function, which
uses units of nanoseconds. If your compiler supports C
interoperability, then this approach is standard conforming (through
C and POSIX). There are some I/O and message passing libraries that
support sleep functionality (usually by using usleep() or
nanosleep() under the hood), so if you don't want to write your own
interface you could link with one of these libraries.

The general advantage of this kind of approach is that the operating
system sets a timer rather than relying on the CPU spinning with
wasted work for the sleep period. Perhaps this has more to do with
esthetics than with practicalities, but I really don't like the idea
of spinning the CPU with do loops that accomplish nothing useful.

$.02 -Ron Shepard

user1

unread,
Dec 21, 2008, 7:17:48 PM12/21/08
to

There are practicalities involved. The method I posted will monopolize the CPU.
A "call delay(10000)" in the code I posted will put 99% cpu usage on the windows
taskmgr for 10 seconds. Fortran could sure use a portable sleep / usleep /
nanosleep function.


Gary Scott

unread,
Dec 21, 2008, 8:06:48 PM12/21/08
to
yeah, but on Windows, the API for timers is minimum of 1 millisecond.
If you want a faster clock, you have to add hardware. I recently
purchased a discrete card with 3 1-usec timers. You can gang them
together to get a larger signal range. Still, i can only read them
about 7 microseconds per byte. So a 16-bit timer takes 14 microseconds
on average (2Ghz emachines T4200). There was a thread about IVF's
sleepqq having to do with the imprecision of the wakeup process. Even
though the resolution is 1 millisecond, it has a fairly large tolerance
factor.

--

Gary Scott
mailto:garylscott@sbcglobal dot net

Fortran Library: http://www.fortranlib.com

Support the Original G95 Project: http://www.g95.org
-OR-
Support the GNU GFortran Project: http://gcc.gnu.org/fortran/index.html

If you want to do the impossible, don't hire an expert because he knows
it can't be done.

-- Henry Ford

e p chandler

unread,
Dec 21, 2008, 11:51:09 PM12/21/08
to
On Dec 17, 9:05 pm, George <geo...@example.invalid> wrote:
[snip]

>
> I broke my bug-reporting cherry today, as I heard from elliot that he gets
> zeros from the systemclock call on g95 as well.  Elsethread, I've been
> trying to work out what the taxonomy of a bug is, and I think this might
> fall into a category where the standard allows an implementation to return
> zero on a systemclock call, but it completely undermines its use to be
> identical from one call to the next here.  

Sorry about that George. I goofed [again]. MR&C says that system_clock
can return 0 on the first call. If so, it is an interval timeer.
I ran your test program, confirmed it returned 0 on the first call,
and let it go at that.

-- e

George

unread,
Dec 22, 2008, 3:51:18 AM12/22/08
to
On Sun, 21 Dec 2008 09:36:15 -0500, user1 wrote:

> Why do you want to reseed the random number generation on each loop iteration ?
> I moved it outside the loop, and the values of u produce on two successive runs
> are listed below.
>


I reseeded to show nonrandom behavior.

Andy's routine with only minor modifications is looking pretty good:

implicit none

real :: u, sum1
integer :: i, j

sum1 = 0.0


do i = 1, 20

call init_seed()

call random_number(u)
print *, " u is ", u
sum1 = sum1 + u

! I need a pause here of 30 milliseconds
do j = 1, 1000000
call random_number(u)
end do

end do

print *, "sum over 20 trials is ", sum1

contains
subroutine init_seed()
integer :: n, ival(8), v(3)
integer, allocatable :: seed(:)

call date_and_time(values=ival)

v(1) = ival(8) + 2048*ival(7)
v(2) = ival(6) + 64*ival(5) ! value(4) isn't really 'random'
v(3) = ival(3) + 32*ival(2) + 32*8*ival(1)

call random_seed(size=n)

allocate(seed(n))

call random_seed() ! Give the seed an implementation-dependent kick

call random_seed(get=seed)

do i=1, n
seed(i) = seed(i) + v(mod(i-1, 3) + 1)
enddo

print *, "seed is ", seed

call random_seed(put=seed)

deallocate(seed)

end subroutine

endprogram


! g95 -Wall duke4.f03 -o h.exe


C:\MinGW\source>g95 -Wall duke4.f03 -o h.exe

C:\MinGW\source>h
seed is 2108286485 -1152138943 -1023650606 1735419315
u is 0.4804541
seed is -936698053 2122440835 -1688078280 -392786557
u is 0.20808166
seed is -942294415 -1451826177 -1743292148 -79510488
u is 0.39475608
seed is 405325042 465516843 -1140849534 -1403313516
u is 0.3664204
seed is -1783054544 -1597365025 -2017858480 1947845786
u is 0.5765057
seed is 355350964 -449964185 1754174651 -154195882
u is 0.32717717
seed is 766811532 495024727 -142102485 -127875334
u is 0.84118825
seed is 1596777905 2029074541 -1569086009 1136507580
u is 0.56891173
seed is 780423977 690355590 1254707161 -695145504
u is 0.35492295
seed is -932995370 -1795388557 -960056318 566888942
u is 0.43557996
seed is 1846858861 261613666 -629179891 -214364602
u is 0.89750624
seed is 143088918 1778048088 -859707451 -1093445463
u is 0.82791734
seed is -1867416918 54633948 1814479614 1644428741
u is 0.53665435
seed is -607328660 -2140823894 -614226858 45923750
u is 0.6718992
seed is 1001648945 1589912045 -369821113 -1009466424
u is 0.3016531
seed is -117694449 -1473505439 115419968 1923804805
u is 0.21691966
seed is -452804293 -751410959 301797425 -966800346
u is 0.07954723
seed is -343129457 1410051559 -218053683 -2050648185
u is 0.6903682
seed is 1102399073 -1553256912 -1468195856 2146152744
u is 0.8821783
seed is -110615999 543205906 -650382501 -1177194024
u is 0.17106211
sum over 20 trials is 9.829703

C:\MinGW\source>

--
George

There's no bigger task than protecting the homeland of our country.

user1

unread,
Dec 22, 2008, 11:29:10 AM12/22/08
to
George wrote:
> On Sun, 21 Dec 2008 09:36:15 -0500, user1 wrote:
>
>> Why do you want to reseed the random number generation on each loop iteration ?
>> I moved it outside the loop, and the values of u produce on two successive runs
>> are listed below.
>>
>
>
> I reseeded to show nonrandom behavior.
>
>

OK, but I don't see the point.

It seems that what you have demonstrated is that the first random number
in the sequence doesn't change very much when the generator is reseeded
using a time of day based seed that hasn't really changed very much.

George

unread,
Dec 22, 2008, 6:28:06 PM12/22/08
to

Exactly, but the way Andy does it, it looks completely random. I was
surprised how much computing fits into 30 milliseconds though.
--
George

It is clear our nation is reliant upon big foreign oil. More and more of
our imports come from overseas.

0 new messages