Mistakes in Fortran 90 Programs That Might Surprise You

FSI · 发表于 2005-11-25 22:15

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

Over the years we have made lots of interesting and fun mistakes in Fortran 90 that we would like to share with you. We welcome your contributions and experiences so that we can share your pain.

<H3>Topics</H3>These "gotchas" are nasty because they will not fail on some machines, while failing on others (given various combinations of compilers and machine platforms).
<UL>
<LI><a>Danger with Optional Arguments</A>
<LI><a>Danger with intent(out)</A>
<LI><a>A suprise with non-advancing I/O</A>
<LI><a>Suprise with locally initialized variables</A>
<LI><a>Danger of calling Fortran 90 style routines</A>
<LI><a>Danger with interfaces to Fortran 77 subroutines</A>
<LI><a>A suprise with generic functions</A>
<LI><a>Big Danger with Undefined Pointers</A>
<LI><a>Subtle danger with overloading (=) to assign pointers</A>
<LI><a>Danger with pointers to pointers</A> </LI></UL>
<H3>Danger with Optional Arguments </H3>In this example an optional argument is used to determine if a header is printed. <PRE> subroutine print_char(this,header)
 character(len=*), intent (in) :: this
 logical, optional, intent (in) :: header
 ! THIS IS THE WRONG WAY
 if (present(header) .and. header) then
 print *, 'This is the header '
 endif
 print *, this
 end subroutine print_char
</PRE>
<PRE> subroutine print_char(this,header)
 character(len=*), intent (in) :: this
 logical, optional, intent (in) :: header
 ! THIS IS THE RIGHT WAY
 if (present(header)) then
 if (header) print *, 'This is the header '
 endif
 print *, this
 end subroutine print_char
</PRE>
Explanation The first method is not safe because the compiler is allowed to evaluate the header argument before the present function is evaluated. If the header argument is not in fact present an out of bounds memory reference could occur, which could cause a failure.
<H3>Danger with intent(out) </H3>In this example we assign components of a derived type with intent(out). <PRE> program intent_gotcha
 type mytype
 integer :: x
 real :: y
 end type mytype

 type (mytype) :: a
 a%x = 1 ; a%y = 2.
 call assign(a)
 ! a%y COULD BE UNDEFINED HERE
 print *, a

 contains

 subroutine assign(this)
 type (mytype), intent (out) :: this
 ! THIS IS THE WRONG WAY
 this%x = 2
 end subroutine assign
</PRE>
<PRE> subroutine assign(this)
 type (mytype), intent (out) :: this
 ! THIS IS THE RIGHT WAY
 this%x = 2 ; this%y = 2.
 end subroutine assign
 end program intent_gotcha
</PRE>
Explanation The problem is that when intent(out) is used with a derived type, any component not assigned in a procedure could become undefined on exit. For example, even though a%y was defined on entry to this routine, it could become undefined on exit because it was never assigned within the routine. The lesson is that all components of a derived type should be assigned within a procedure, when intent(out) is used. Intent(out) behaves like the result variable in a function: all components must be assigned.
As an alternative, use intent(inout).
<H3>A suprise with non-advancing I/O </H3>Many people think that the new non-advancing I/O in Fortran 90 is the same as stream I/O in other languages. It is not. <PRE> do i = 1, 128
 write (unit=6,fmt='(a)',advance='no') 'X'
 end do
</PRE>
We expect this program to print 128 X's in a row. However, unexpected behavior may occur if the record length for unit 6 is less than 128.
One can inquire the record length in the follow way:
<PRE> open(unit=6)
 inquire(unit=6, recl=i)
 print *, 'recl =', i
</PRE>
Explanation All Fortran I/O is still record based. Non-advancing I/O allows partial reads and writes within a record. For many compilers the default record length is very large (e.g., 2147483647) giving the appearance of stream I/O. This is not true for all compilers however.
On some compilers it is possible to set the record length as follows: <PRE>open(unit=6, recl = 2147483646)
</PRE>
On other compilers unit 6 is preconnected and the record length cannot be changed. (Thanks to Jon Richards of the USGS for this tip.)
Note that unit 6 and unit * are not necessarily the same. Although they both may point to the default output device, with non-advancing I/O, each could keep track of the current location in its own record separately. Therefore we advise choosing one default unit and sticking with it.
<H3>Suprise with locally initialized variables </H3>One must be careful when initializing a locally declared variable. <PRE> real function kinetic_energy(v)
 real, dimension(:), intent(in) :: v
 integer i
 ! THIS IS THE WRONG WAY
 real :: ke = 0.0
 do i = 1, size(v)
 ke = ke + v(i)**2
 enddo
 kinetic_energy = .5*ke
 end function kinetic_energy
</PRE>
<PRE> real function kinetic_energy(v)
 real, dimension(:), intent(in) :: v
 integer i
 ! THIS IS THE RIGHT WAY
 real :: ke
 ke = 0.
 do i = 1, size(v)
 ke = ke + v(i)**2
 enddo
 kinetic_energy = .5*ke
 end function kinetic_energy
</PRE>
Explanation A local variable that is initialized when declared has an implicit save attribute. ke is initialized only the first time the function is called. On subsequent calls the old value of ke is retained. This is a real suprise to C programmers.
To avoid confusion it is best to add the save attribute to such locally initialized variables explicitly, even though this is redundant.
<H3>Danger of calling Fortran 90 style routines </H3><PRE> program main
 real, dimension(5) :: x

 x = 0.
! THIS IS WRONG
 call incb(x)
 print *, x

 end program main

 subroutine incb(a)
! this is a fortran90 style subroutine
 real, dimension(:) :: a
 a = a + 1.
 end subroutine incb
</PRE>
Explanation The subroutine incb uses a Fortran 90 style assumed shape array (containing dimension(:)). Such routines must either be in a module, or have an explicit interface wherever they are used. In this example, neither one was true.
One correct way to call such procedures is to use an explicit interface as follows: <PRE> program main
 real, dimension(5) :: x

! THIS IS THE RIGHT WAY
 interface
 subroutine incb(a)
 real, dimension(:) :: a
 end subroutine incb
 end interface

 x = 0.
 call incb(x)
 print *, x

 end program main

 subroutine incb(a)
! this is a fortran90 style subroutine
 real, dimension(:) :: a
 a = a + 1.
 end subroutine incb
</PRE>
If the routine is in a module interfaces are generated automatically and do not need to be explicitly written. <PRE>! THIS IS ANOTHER RIGHT WAY
 module inc
 contains
 subroutine incb(a)
! this is a fortran90 style subroutine
 real, dimension(:) :: a
 a = a + 1.
 end subroutine incb
 end module inc

 program main
 use inc
 real, dimension(5) :: x

 x = 0.
 call incb(x)
 print *, x

 end program main
</PRE>If interfaces are used, the interface MUST match the actual function. 
<H3>Danger with interfaces to Fortran 77 subroutines </H3><PRE> program main
 real, dimension(5) :: x

! interface to Fortran 77 style routine
 interface
 subroutine inca(a,n)
 integer :: n
! THIS IS THE WRONG WAY
 real, dimension(:) :: a
! THIS IS THE RIGHT WAY
 real, dimension(n) :: a
 end subroutine inca
 end interface

 x = 0.
 call inca(x,5)
 print *, x

 end program main

 subroutine inca(a,n)
! this is a fortran77 style subroutine
 dimension a(n)
 do 10 j = 1, n
 a(j) = a(j) + 1.
10 continue
 return
 end
</PRE>
Explanation The interface declaration must always match the actual subroutine declaration. In this case, the interface statement refers to a Fortran 90 style assumed shape array. The actual subroutine refers to a Fortran 77 explicit shape array. The lesson here is: Interfaces to Fortran 77 style routines must only use Fortran 77 style constructs.
In this example, it is permitted to leave out the interface altogether since routines without interfaces are treated as Fortran77 style routines by default. However, if the interface is left out, the compiler will no longer check whether the arguments of calling procedures agree with the arguments listed in the interface.

FSI · 发表于 2005-11-25 22:15

<H3>A Suprise with Generic Functions (Function Overloading) </H3>Fortran 90 allows the same function name to be used for different actual functions, so long as the arguments to the functions differ. One would expect that the functions first_sub and second_sub below would be different, because in first_sub, the first argument is a real and the second is an integer, while in second_sub the arguments are reversed.
<PRE> subroutine first_sub(a,i)
 real :: a
 integer :: i
 ...
 end subroutine first_sub
!
 subroutine second_sub(i,a)
 integer :: i
 real :: a
 ...
 end subroutine second_sub
</PRE>
So that one could define a generic function first_or_second below: <PRE> interface first_or_second
 module procedure first, second
 end interface
</PRE>
This is NOT so.
Explanation The reason is that Fortran 90 allows procedures to be called by name (keyword) arguments. The following <PRE> real :: b
 integer :: n
 call first_or_second(i=n,a=b)
</PRE>
does not work because when called by keyword, first_sub and second_sub are indistinguishable, <PRE> call first_sub(i=n,a=b)
 call second_sub(i=n,a=b)
</PRE>
and therefore a generic function cannot be defined. A generic function must be able to distinguish its arguments by type AND by name.
The solution is to not use the same dummy argument name in both procedures. For example, the following would work: <PRE> subroutine second_sub(i,aa)
 integer :: i
 real :: aa
 ...
 end subroutine second_sub
</PRE>

<H3>Dangers with Pointers </H3>Fortran 90 has 3 ways to implement dynamic memory: Automatic arrays, allocatable arrays, and pointers.
Automatic arrays are automatically created on entry and deleted on exit from a procedure, and they are safest and easiest to use. Allocatable arrays require the user to manually create and delete them, and should only be used if automatic creation and deletion is not the desired behavior.
Pointers are the most error prone and should only be used when allocatable arrays are not possible, e.g., when one desires an array to be a component of a derived type.
<H3>Big Danger with Undefined Pointers </H3>Many people think that the status of a pointer which has never been associated is .not. associated. This is false.
In this example we are allocating a local_table on first entry that is to be reused on subsequent entries. <PRE> subroutine local_pointer(this)
 real, dimension(:) :: this
 real, dimension(:), save, pointer :: local_table
 ! THIS IS THE WRONG WAY
 if (.not. associated(local_table)) then
 allocate(local_table(size(this)))
 endif
 local_table = ...
 ...
 end subroutine local_pointer
</PRE>
<PRE> subroutine local_pointer(this)
 real, dimension(:) :: this
 real, dimension(:), save, pointer :: local_table
 ! THIS IS THE RIGHT WAY
 logical, save :: first_entry = .true.
 if (first_entry) then
 nullify(local_table) ; first_entry = .false.
 end if
 if (.not. associated(local_table)) then
 allocate(local_table(size(this)))
 endif
 local_table = ...
 ...
 end subroutine local_pointer
</PRE>
Explanation When a pointer is declared its status is undefined, and cannot be safely queried with the associated intrinsic. A second variable is introduced to nullify the pointer on first entry so that its status can be safely tested. This is not a problem in Fortran 95 which allows one to nullify a pointer on declaration.
Note that the save attribute for local_table is necessary to guarantee that the array and the pointer status are preserved on subsequent entries. We recommend that the save attribute should always be used when pointers and allocatable arrays are allocated in procedures.
<H3>Subtle danger with overloading (=) to assign pointers </H3>One must be careful with overloading the assignment operator.
In this module we have created a private type which contains a pointer and a public procedure to assign that pointer. <PRE> module assign_pointer_class
 type mytype
 private
 real, pointer :: pr
 end type mytype
 interface assignment (=)
 module procedure assign_pointer
 end interface
 contains
 subroutine assign_pointer(this, a)
 type (mytype), intent(out) :: this
 real, target, intent(in) :: a
 this%pr => a
 end subroutine assign_pointer
 end module assign_pointer_class
</PRE>
In this main program we intend to assign the pointer component x%pr to the variable a, x%pr =>a. We cannot do so directly because the components of mytype are private. One must use a public procedure to do so. Furthermore, to simplify the syntax one might be tempted to use an overloaded assignment operator (=). <PRE> program main
 use assign_pointer_class
 type (mytype) :: x
 real :: a = 0
 ! THIS IS THE WRONG WAY
 x = a
 end program main
</PRE>
Don't give into this temptation! The only safe way to accomplish this is to call the procedure directly. <PRE> program main
 use assign_pointer_class
 type (mytype) :: x
 ! THIS IS THE RIGHT WAY
 real, target :: a = 0
 call assign_pointer(x,a)
 end program main
</PRE>
Explanation The Fortran 90 standard says that the right hand side of an assignment operator is an expression that may potentially only persist for the duration of the call. In other words, x%pr could inadvertently point to a temporary copy of the variable a.
Thanks to Henry Zongaro of IBM for pointing this out. (We never would have figured this one out on our own.)
Also, James Giles found a subtle point regarding this example. We did not include "target" in the declaration of the real variable "a" (this has been corrected above). In James' words:
"Notice that for this to really work, the actual argument, 'a', must be declared with the target attribute. You correctly declare the dummy argument in the assign_pointer routine with the target attribute, but the actual argument must also have that attribute (otherwise it's illegal for any pointer to be associated with it). Just a minor point..."

<H3>Danger with pointers to pointers </H3>When creating a hierarchy of pointers to pointers, each level of pointers must be allocated before being used. <PRE> program main

 type mytype
 real, dimension(:), pointer :: p
 end type mytype

 type (mytype), pointer :: x


! BOTH OF THESE ARE THE WRONG WAY
! AND THE COMPILER WON'T CATCH IT
! nullify(x%p)
! allocate(x%p(5))


! ONE SHOULD ALWAYS IMMEDIATELY NULLIFY THE PARENT POINTER
! OR ALLOCATE IT
 nullify(x) ! or allocate(x)
 ...
! THEN LATER NULLIFY OR ALLOCATE THE CHILD POINTER
 call child_construct(x,5)
 if (associated(x%p)) print *, x%p

 contains

 subroutine child_construct(this,len)
! child constructor for pointer within mytype
! if len is present, then allocate it, otherwise nullify it.
! mytype is assumed to be already nullified or allocated
 type (mytype), pointer :: this
 integer, optional, intent(in) :: len
 if (.not.associated(x)) allocate(x)
 if (present(len)) then
 allocate(x%p(len))
 x%p = 0.
 else
 nullify(x%p)
 endif
 end subroutine child_construct

 end program main
</PRE>
Explanation This example creates a pointer to a pointer to an array of reals where the first pointer has not been allocated. For safety one should always either allocate or nullify the parent pointer immediately after its declaration. The child pointer cannot be allocated before the parent. Since the child pointer may be allocated elsewhere in the code, it is convenient to use constructor routines for this purpose.
Each child constructor can safely allocate or nullify its pointers only when it can be sure that its parent's pointers have been allocated or nullified.

账号		自动登录	找回密码
密码			我要加入

[其他相关] Mistakes in Fortran 90 Programs That Might Surprise You

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。