Re: Is this a compiler error?



First order of business is to change the 'pure' in the functions x, y, z, since this would give a compiler error or it should give a compiler error. I've decided for completeness I should include the random number generator, but, if you want to comment about it let's start a new thread, it isn't really relevant to this thread.

System:



module base_kinds
implicit none

! Default kinds
integer, parameter:: ip = kind(1)
integer, parameter:: lp = kind(.FALSE.)
integer, parameter:: sp = kind(1.0)
integer, parameter:: dp = kind(1.0D0)

end module

program main
use base_kinds
use my_functions
implicit none

real(kind=sp):: r, t
integer:: s(4)

s = 20
t = 0

call initialize(s)
r = 2*x(t) + y(t)
write(*,*) t, r

call initialize(s)
r = y(t) + 2*x(t)
write(*,*) t, r

end program

module my_functions
use base_kinds
use random
implicit none

type(a_random):: p

contains

subroutine initialize(s)
integer, intent(in):: s(4)
call random_s(p, s)
end subroutine

function x(t)
real(kind=sp):: x
real(kind=sp), intent(in):: t
x = t
call random_sp(p, x)
end function

function y(t)
real(kind=sp):: y
real(kind=sp), intent(in):: t
y = t
call random_sp(p, y)
end function

function z(t)
real(kind=sp):: z
real(kind=sp), intent(in):: t
z = t
call random_sp(p, z)
end function

end module

MODULE random
use base_kinds
implicit none
PRIVATE

! The model parameters
! Computes a random sequence using a linear congruential generator
! d(n+1) = m*d(n)+d0 mod 2**bits
! For the sequence d(n) to reproduce all the whole numbers in the model, 2**bits,
! mod(d0, 2) = 1
! and mod(m, 4) = 1
! The trick is to embed the whole numbers in the default integer and perform all
! the calculations without overflowing the default integer. To maintain a large bit
! size and, therefore, a large cycle, the calculation is done on an array of integers
! each holding a 'digit' of the large integer.
! To prevent overflow it is necessary to ignore the sign bit and the addition will need
! up to two carry bits. So m must be chosen to not exceed model_bits - 3 - bits.
! m = 2**p + 1
! p = model_bits - 3 - bits.
! So if each 'digit' is 24 bits embedded in a 32 bit integer then
! p = 32 - 3 -24 = 5
! m = 2**5 + 1 = 33
! So the model is an array of 4 'digits', with each digit being a 24 bit whole number
! embeddied in a 32 bit integer. The complete whole number is 4*24 = 96 bits
! and has a cycle of 2**96.
integer, parameter:: model_bits = 31 ! ignoring the sign bit
integer, parameter:: bits = 24 ! each digit is 24 bits
integer, parameter:: up = model_bits - 2 - bits ! saving two carry bits
integer, parameter:: d0(4) = (/Z'DDDDDD', Z'BBBBBB', Z'999999', Z'777777'/)
integer, parameter:: mask = 2**bits - 1
real(kind=sp), parameter:: modelmax_sp = 2.0_sp**model_bits
real(kind=dp), parameter:: modelmax_dp = 2.0_dp**model_bits

type a_random
private
integer:: d(4)
end type

PUBLIC a_random, random_s, random_i, random_sp, random_dp

CONTAINS

pure subroutine random_s(p,s)
type(a_random), intent(inout):: p
integer, intent(in):: s(4)
integer:: i
! Seeds a 96 bit integer given as 4-24 bit integers, stored in 4-32 bit integers.
do i = 1, 4
p%d(i) = iand(s(i), mask)
enddo
end subroutine

pure subroutine random_inc(p)
type(a_random), intent(inout):: p
integer:: i, ovf
! 96 bit computation of (65*d+d0)(mod 96 bit)
! This operation applied sequentially has a period of 2**96
! and cycles through all 96 bit whole numbers.
ovf = 0
do i = 1, 4
p%d(i) = ishft(p%d(i), up) + p%d(i) + d0(i) + ovf
ovf = ishft(p%d(i), -bits)
p%d(i) = iand(p%d(i), mask)
enddo
end subroutine

pure subroutine random_i(p, r)
type(a_random), intent(inout):: p
integer, intent(out):: r
integer:: i1, i2
! Result is a random 31 bit whole number.
! ieor is an addition without the carry.
call random_inc(p)
i1 = ieor(p%d(1),p%d(3))
i2 = ishft(ieor(p%d(2), p%d(4)), model_bits - bits)
r = ieor(i1, i2)
end subroutine

pure subroutine random_sp(p, r)
type(a_random), intent(inout):: p
real(kind=sp), intent(out):: r
integer:: i
call random_i(p, i)
r = i/modelmax_sp
end subroutine

pure subroutine random_dp(p, r)
type(a_random), intent(inout):: p
real(kind=dp), intent(out):: r
integer:: i1, i2
call random_inc(p)
i1 = ieor(ishft(p%d(4), model_bits - bits), p%d(3))
i2 = ieor(ishft(p%d(2), model_bits - bits),p%d(1))
r = (i1+i2/modelmax_dp)/modelmax_dp
end subroutine

END MODULE random

One thing I wanted from the program is that the 'error' should be eye opening. Even though I understand Richard's comment about Fortran not being associative and commutative, I must stress that this is an obvious breakdown that must be dealt with. This is all the thunder I've got, please, feign a little awe. Obviously, Fortran has to follow the mathematics at least approximately.

You can get a break with the underlying mathematics because of the finite representation of type real. In the most coarse case the problem could deal with numbers outside the range of representation, overflow and underflow. It used to be these conditions would be runtime errors, but with IEEE you can get silent infinities and silent NAN's. Integers are even worse. It is very common that integer overflow is silent (the processor can catch these, but it is usually turned off). In fact, I used to write random number generators and just ignore integer overflow. Trouble is different compilers can use different ways of handling integer overflow, so the results were compiler dependent.

Another way the finite repesentation causes trouble is that Fortran follows the math only to within the round off error. The trouble here is the round off error depends on the problem. As Francois noted the difference this program produces appears to be quite large for round off error. You would have to do a lot of work before the round off error got that big. Of course, that's possible, so this should certainly be kept in mind. Round off error can cause a complete break of the stored number and the underlying math. It's like using a 4 byte real to store the sum of 1--50,000,000 times and not coming out with 50,000,000.

Now, what I consider to be the most common break with the underlying math is this side effects problem. Judging from the posts everyone else considers this the most common. So it quickly becomes the candidate of choice when presented with this difficulty.

Note it is possible that a compiler could evaluate these two expressions

r = 2*x+y
r = y+2*x

using exactly the same code. The compiler is free to use any mathematically equivalent expression, sort of like 'Send a picture of you or somebody who looks just like you'. In which case you wouldn't see the error at all. Point is, this problem doesn't necessarily produce two different numbers, but what it produces is compiler dependent. I was tickled to find you could determine how the compiler evaluated the expression with this problem. It was evaluating the expression from left to right as written.

MikeT
.



Relevant Pages

  • Re: How can I make this module less repetitive?
    ... now I'm using an include file to store the common body between two subroutines. ... pure subroutine ck54_val ... end interface ...
    (comp.lang.fortran)
  • Re: DPROD issues
    ... a switch like that typically ... makes a compiler nonstandard in that mode. ... treatment of specific intrinsics is one ... subroutine sub1a ...
    (comp.lang.fortran)
  • Re: Jumping into block of an if construct
    ... (For that matter a clever enough compiler could replace this PUT DATA ... routine which itself executes the loops around element handling. ... Either way I think the cost of element handling will usually ... So locally based on subroutine arguments, but not on, for example, ...
    (comp.lang.fortran)
  • Re: Bus error/ segmentation fault--help?
    ... When I compile with the intel fortran compiler, ... This subroutine integrates the function y3 up one step, ... implicit double precision (a-h, o-z) ...
    (comp.lang.fortran)
  • Re: Question about name conflicts in Fortran
    ... I am hacking the g95 compiler so that it dumps the structure ... Subroutine foo ... The para after the 2nd numbered list in 16.2 covers the case where the ... One of the exceptions is "a generic name may be the ...
    (comp.lang.fortran)