常见fortran系统错误

Mistakes in Fortran 90 Programs That Might Surprise YouOver 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.TopicsThese "gotchas" are nasty because they will not fail on some machines, while failing on others (given various combinations of compilers and machine platforms).∙Danger with Optional Arguments∙Danger with intent(out)∙ A suprise with non-advancing I/O∙Suprise with locally initialized variables∙Danger of calling Fortran 90 style routines∙Danger with interfaces to Fortran 77 subroutines∙ A suprise with generic functions∙Big Danger with Undefined Pointers∙Subtle danger with overloading (=) to assign pointers∙Danger with pointers to pointersDanger with Optional ArgumentsIn this example an optional argument is used to determine if a header is printed.subroutine print_char(this,header)character(len=*), intent (in) :: thislogical, optional, intent (in) :: header! THIS IS THE WRONG WAYif (present(header) .and. header) thenprint *, 'This is the header 'endifprint *, thissubroutine print_char(this,header)character(len=*), intent (in) :: thislogical, optional, intent (in) :: header! THIS IS THE RIGHT WAYif (present(header)) thenif (header) print *, 'This is the header 'endifprint *, thisend subroutine print_charExplanationThe 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.Danger with intent(out)In this example we assign components of a derived type with intent(out). program intent_gotchatype mytypeinteger :: xreal :: yend type mytypetype (mytype) :: aa%x = 1 ; a%y = 2.call assign(a)! a%y COULD BE UNDEFINED HEREprint *, acontainssubroutine assign(this)type (mytype), intent (out) :: this! THIS IS THE WRONG WAYthis%x = 2end subroutine assignsubroutine assign(this)type (mytype), intent (out) :: this! THIS IS THE RIGHT WAYthis%x = 2 ; this%y = 2.end subroutine assignExplanationThe 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).A suprise with non-advancing I/OMany 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.do i = 1, 128write (unit=6,fmt='(a)',advance='no') 'X'end doWe 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:open(unit=6)inquire(unit=6, recl=i)print *, 'recl =', iExplanationAll 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:open(unit=6, recl = 2147483646)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.Suprise with locally initialized variablesOne must be careful when initializing a locally declared variable.real function kinetic_energy(v)real, dimension(:), intent(in) :: vinteger i! THIS IS THE WRONG WAYreal :: ke = 0.0do i = 1, size(v)ke = ke + v(i)**2enddokinetic_energy = .5*keend function kinetic_energyreal function kinetic_energy(v)real, dimension(:), intent(in) :: vinteger i! THIS IS THE RIGHT WAYreal :: keke = 0.do i = 1, size(v)ke = ke + v(i)**2enddokinetic_energy = .5*keend function kinetic_energyExplanationA 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.Danger of calling Fortran 90 style routinesprogram mainreal, dimension(5) :: xx = 0.! THIS IS WRONGcall incb(x)print *, xend program mainsubroutine incb(a)! this is a fortran90 style subroutinereal, dimension(:) :: aa = a + 1.end subroutine incbExplanationThe 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:program mainreal, dimension(5) :: x! THIS IS THE RIGHT WAYinterfacesubroutine incb(a)real, dimension(:) :: aend subroutine incbend interfacex = 0.call incb(x)print *, xend program mainsubroutine incb(a)! this is a fortran90 style subroutinereal, dimension(:) :: aa = a + 1.end subroutine incbIf the routine is in a module interfaces are generated automatically and do not need to be explicitly written.! THIS IS ANOTHER RIGHT WAYmodule inccontainssubroutine incb(a)! this is a fortran90 style subroutinereal, dimension(:) :: aa = a + 1.end subroutine incbend module incprogram mainuse increal, dimension(5) :: xx = 0.call incb(x)print *, xend program mainIf interfaces are used, the interface MUST match the actual function. Danger with interfaces to Fortran 77 subroutinesprogram mainreal, dimension(5) :: x! interface to Fortran 77 style routineinterfacesubroutine inca(a,n)integer :: n! THIS IS THE WRONG WAYreal, dimension(:) :: a! THIS IS THE RIGHT WAYreal, dimension(n) :: aend subroutine incaend interfacex = 0.call inca(x,5)print *, xend program mainsubroutine inca(a,n)! this is a fortran77 style subroutinedimension a(n)do 10 j = 1, na(j) = a(j) + 1.10 continuereturnendExplanationThe 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.A Suprise with Generic Functions (Function Overloading)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. subroutine first_sub(a,i)real :: ainteger :: i...end subroutine first_sub!subroutine second_sub(i,a)integer :: ireal :: a...end subroutine second_subSo that one could define a generic function first_or_second below:interface first_or_secondmodule procedure first, secondend interfaceThis is NOT so.ExplanationThe reason is that Fortran 90 allows procedures to be called by name (keyword) arguments. The followingreal :: binteger :: ncall first_or_second(i=n,a=b)does not work because when called by keyword, first_sub and second_sub are indistinguishable,call first_sub(i=n,a=b)call second_sub(i=n,a=b)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:subroutine second_sub(i,aa)integer :: ireal :: aa...end subroutine second_subDangers with PointersFortran 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.Big Danger with Undefined PointersMany 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.subroutine local_pointer(this)real, dimension(:) :: thisreal, dimension(:), save, pointer :: local_table! THIS IS THE WRONG WAYif (.not. associated(local_table)) thenallocate(local_table(size(this)))endiflocal_table = ......end subroutine local_pointersubroutine local_pointer(this)real, dimension(:) :: thisreal, dimension(:), save, pointer :: local_table! THIS IS THE RIGHT WAYlogical, save :: first_entry = .true.if (first_entry) thennullify(local_table) ; first_entry = .false.end ifif (.not. associated(local_table)) thenallocate(local_table(size(this)))endiflocal_table = ......end subroutine local_pointerExplanationWhen 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.Subtle danger with overloading (=) to assign pointersOne 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.module assign_pointer_classtype mytypeprivatereal, pointer :: prend type mytypeinterface assignment (=)module procedure assign_pointerend interfacecontainssubroutine assign_pointer(this, a)type (mytype), intent(out) :: thisreal, target, intent(in) :: athis%pr => aend subroutine assign_pointerend module assign_pointer_classIn 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 (=).program mainuse assign_pointer_classtype (mytype) :: xreal :: a = 0! THIS IS THE WRONG WAYx = aend program mainDon't give into this temptation! The only safe way to accomplish this is to call the procedure directly.program mainuse assign_pointer_classtype (mytype) :: x! THIS IS THE RIGHT WAYreal, target :: a = 0call assign_pointer(x,a)end program mainExplanationThe 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..."Danger with pointers to pointersWhen creating a hierarchy of pointers to pointers, each level of pointers must be allocated before being used.program maintype mytypereal, dimension(:), pointer :: pend type mytypetype (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 ITnullify(x) ! or allocate(x)...! THEN LATER NULLIFY OR ALLOCATE THE CHILD POINTERcall child_construct(x,5)if (associated(x%p)) print *, x%pcontainssubroutine 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 allocatedtype (mytype), pointer :: thisinteger, optional, intent(in) :: lenif (.not.associated(x)) allocate(x)if (present(len)) thenallocate(x%p(len))x%p = 0.elsenullify(x%p)endifend subroutine child_constructend program mainExplanationThis 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.。