Interface block inside a subroutine.
Daniel Carrera
I've written a Runge-Kutta integrator as a function and it works well.
Now I want to make it into a subroutine and I'm having trouble. For some
reason, it looks like the interface block is not welcome inside the
subroutine. This is my code:
pure subroutine rk4_vec(t, Y, dY, h)
real(pref), intent(inout) :: t, Y(:)
real(pref), intent(in) :: dY(size(Y)), h
real(pref), dimension(size(Y)) :: k1, k2, k3, k4
interface
pure function dY(t0, y0)
use types
real(pref), intent(in) :: t0, y0(size(Y))
real(pref) :: dY(size(Y))
end function
end interface
k1 = dY(t, Y)
k2 = dY(t + h/2, Y + k1*h/2)
k3 = dY(t + h/2, Y + k2*h/2)
k4 = dY(t + h , Y + k3*h)
Y = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6
t = t + h
end subroutine
The compiler errors are:
========================== // ======================
modules/ode.f90:74.32:
real(pref) :: dY
1
Error: Symbol 'dy' at (1) already has basic type of REAL
modules/ode.f90:93.36:
pure function dY(t0, y0)
1
Error: 'dy' at (1) has attributes specified outside its INTERFACE body
modules/ode.f90:94.25:
use types
1
Error: Unexpected USE statement in INTERFACE block at (1)
modules/ode.f90:95.21:
real(pref), intent(in) :: t0, y0(size(Y))
1
Error: Parameter 'pref' at (1) has not been declared or is a variable,
which does not reduce to a constant expression
modules/ode.f90:96.21:
real(pref) :: dY(size(Y))
1
Error: Parameter 'pref' at (1) has not been declared or is a variable,
which does not reduce to a constant expression
modules/ode.f90:97.15:
end function
========================== // ======================
Ok, so it doesn't like the type declaration of dY inside the interface
block, or for that matter, almost anything that I put inside the
interface block. But this same interface block works fine if I the outer
procedure is a function instead of a subroutine. What gives?
Cheers,
Daniel.
steve
> pure subroutine rk4_vec(t, Y, dY, h)
use types ! Presumably this is available
> real(pref), intent(inout) :: t, Y(:)
> real(pref), intent(in) :: dY(size(Y)), h
dY(size(Y)) ....
> real(pref), dimension(size(Y)) :: k1, k2, k3, k4
>
> interface
> pure function dY(t0, y0)
... dY(t0,y0) ?
--
steve
Richard Maine
> I've written a Runge-Kutta integrator as a function and it works well.
> Now I want to make it into a subroutine and I'm having trouble. For some
> reason, it looks like the interface block is not welcome inside the
> subroutine. This is my code:
>
> pure subroutine rk4_vec(t, Y, dY, h)
> real(pref), intent(inout) :: t, Y(:)
> real(pref), intent(in) :: dY(size(Y)), h
> real(pref), dimension(size(Y)) :: k1, k2, k3, k4
>
> interface
> pure function dY(t0, y0)
> use types
> real(pref), intent(in) :: t0, y0(size(Y))
> real(pref) :: dY(size(Y))
> end function
> end interface
The interface block is fine - well, except for a detail I'll note below.
What is not so fine is that you have two separate declarations of DY -
one in the interface block and one outside. You can't do that. Drop the
one outside. The interface block already declares everything relevant
about the function; you don't also redo some of the declarations outside
of the interface.
I'm also slightly confused by the intent(in) for the first declaration
of dY. Intent only applies to data arguments and pointers. dY is a
function instead of a data argument. Intent makes no sense (and is not
alowed) for a dumy function. But since you are going to delete that
declaration of dY completely, the intent(in) should become moot.
The one thing that isn't so fine in the interface block is that the
interface body in it has references to size(Y). Interface bodies do
*NOT* inherit by host association (unless you use the f2003 IMPORT
statement, which is, in my view, a hack around a design flaw in f90). So
nothing about Y (including its size) is known in the interface body.
What you probably want is for y0 to be assumed shape (i.e. dimensioned
(:)); it will get its shape from the actual arguments to the DY
function. Yes, that indirectly gets it from Y, but via argument
association, which is ok, instead of via host association, which is not.
Then you can dimension dY as size(y0) instead of size(Y); that's inside
the interface body - not the declaration outside the interface body,
which you need to delete.
Also, as Steve noted, you probably need a "use types" at the top of
subroutine rk4_vec (in addition to the one in the interface block; don't
delete that one).
Several of the other error messages are probably cascades from some of
the above.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
Daniel Carrera
> The interface block is fine - well, except for a detail I'll note below.
> What is not so fine is that you have two separate declarations of DY -
> one in the interface block and one outside. You can't do that. Drop the
> one outside. The interface block already declares everything relevant
> about the function; you don't also redo some of the declarations outside
> of the interface.
Thanks! Now it seems to work perfectly.
> I'm also slightly confused by the intent(in) for the first declaration
> of dY. Intent only applies to data arguments and pointers. dY is a
> function instead of a data argument. Intent makes no sense (and is not
> alowed) for a dumy function. But since you are going to delete that
> declaration of dY completely, the intent(in) should become moot.
Ok. I thought "intent(in)" was just something you use to help make sure
that you never overwrite a variable.
> The one thing that isn't so fine in the interface block is that the
> interface body in it has references to size(Y)... So
> nothing about Y (including its size) is known in the interface body.
> What you probably want is for y0 to be assumed shape (i.e. dimensioned
> (:)); it will get its shape from the actual arguments to the DY
> function.
I think I understand. I certainly understand the problem: That nothing
about Y is known inside the interface body.
As for the solution, I think I mostly understand. I just don't feel I
fully understand assumed shape arrays. In particular, when I actually
write the DY function, I won't use assumed shape arrays. Instead I'll
write something like:
pure function dr(t, r)
real(pref), intent(in) :: t, r(2)
real(pref) :: dr(2)
dr = (/ -sin(t) , cos(t) /)
end function
But I think that I mostly understand why that's ok. The interface...
interface
pure function dY(t0, y0)
use types
real(pref), intent(in) :: t0, y0(:)
real(pref) :: dY(size(y0))
end function
end interface
... only says that DY and Y0 have the same shape.
> Also, as Steve noted, you probably need a "use types" at the top of
> subroutine rk4_vec (in addition to the one in the interface block; don't
> delete that one).
I currently have it at the top of the module. Is that enough? It seems
to work...
> Several of the other error messages are probably cascades from some of
> the above.
Indeed, they are. Removing the first declaration of DY made all the
errors go away.
Thanks for the help.
Cheers,
Daniel.
glen herrmannsfeldt
(snip)
> Ok. I thought "intent(in)" was just something you use to help make sure
> that you never overwrite a variable.
Starting I believe with Fortran 2003 there are procedure pointer
variables, which could be passed to a called procedure.
Otherwise, from Fortran II through Fortran 95, you can pass a
function or subroutine as an actual argument, but there is no
variable.
(snip)
> As for the solution, I think I mostly understand. I just don't feel I
> fully understand assumed shape arrays. In particular, when I actually
> write the DY function, I won't use assumed shape arrays. Instead I'll
> write something like:
> pure function dr(t, r)
> real(pref), intent(in) :: t, r(2)
> real(pref) :: dr(2)
> dr = (/ -sin(t) , cos(t) /)
> end function
> But I think that I mostly understand why that's ok. The interface...
The interface should agree with the actual function. If you
don't use assumed shape, then the interface shouldn't either.
Change to either (2) or (*) as appropriate.
-- glen
Richard Maine
> Ok. I thought "intent(in)" was just something you use to help make sure
> that you never overwrite a variable.
It is (and it can also help optimizers). But your dy is not a variable;
that's exactly the problem. It is a dummy procedure, which is not a
variable. Overwriting a dummy procedure can't happen anyway; there isn't
any syntax for such a thing in the first place. (A dummy procedure
pointer would be a different matter again).
> As for the solution, I think I mostly understand. I just don't feel I
> fully understand assumed shape arrays. In particular, when I actually
> write the DY function, I won't use assumed shape arrays. Instead I'll
> write something like:
>
> pure function dr(t, r)
> real(pref), intent(in) :: t, r(2)
> real(pref) :: dr(2)
>
> dr = (/ -sin(t) , cos(t) /)
> end function
>
>
> But I think that I mostly understand why that's ok. The interface...
>
> interface
> pure function dY(t0, y0)
> use types
> real(pref), intent(in) :: t0, y0(:)
> real(pref) :: dY(size(y0))
> end function
> end interface
>
>
> ... only says that DY and Y0 have the same shape.
Uh Oh. NO!, NO! NO! That is *NOT* the only thing it says. Very
importantly, it also says that y0 is assumed shape. If the actual
procedure does not have y0 as assumed shape, all kinds of evil is likely
to result.
An interface body *MUST* describe the same interface as the procedure
actually has. A nice thing about module procedures is that you don't
have to worry about getting the interface body right. Alas, that doesn't
help you with dummy procedures. (The f2003 procedure statement can help,
but that's an f2003 feature).
The standard has a boring list of what things constitute the interface,
but for most purposes, you can just forget the list and assume that
everything must match. (The main exception is that dummy argument names
don't actually have to match, but I think one really ought to make them
match as well unless there is some good reason otherwise.) Explicit
shape versus assumed shape *DEFINITELY* has to match. It changes how the
underlying call works; namely for assumed shape the compiler has to pass
the shape information.
You might have gotten it to compile, but having it actually work
correctly is a different matter. And if it happens to work on a
particular compiler, that doesn't guarantee it for others. Having the
interface body disagree with the actual procedure violates the standard.
I suppose the main question is why you don't use assumed shape in the
actual procedure. That seems the most obvious fix.
> > Also, as Steve noted, you probably need a "use types" at the top of
> > subroutine rk4_vec (in addition to the one in the interface block; don't
> > delete that one).
>
> I currently have it at the top of the module. Is that enough? It seems
> to work...
Yes. You just didn't show that. Context matters.
Daniel Carrera
>> But I think that I mostly understand why that's ok. The interface...
>>
>> interface
>> pure function dY(t0, y0)
>> use types
>> real(pref), intent(in) :: t0, y0(:)
>> real(pref) :: dY(size(y0))
>> end function
>> end interface
>>
>>
>> ... only says that DY and Y0 have the same shape.
>
> Uh Oh. NO!, NO! NO! That is *NOT* the only thing it says. Very
> importantly, it also says that y0 is assumed shape. If the actual
> procedure does not have y0 as assumed shape, all kinds of evil is likely
> to result.
Ok. I'm glad I asked.
I want to write a general-purpose ODE integrator that can help me with
different problems. For example:
1) My simplest test is to draw a circle:
!
! This ODE gives me a circle.
!
pure function dr(t, r)
real(pref), intent(in) :: t, r(2)
real(pref) :: dr(2)
dr = (/ -sin(t) , cos(t) /)
end function
2) I want to model a test particle inside the gravity well of the Sun
(2-dimensional problem):
pure function d_vec(t, vec)
real(pref), intent(in) :: t, vec(4)
real(pref) :: V_x, V_y, r, d_vec(4), g, x, y
x = vec(1) ; V_x = vec(3) ! X position and velocity.
y = vec(2) ; V_y = vec(4) ! Y position and velocity.
r = sqrt( x*x + y*y )
!
! I have a function to compute the gravitational field.
!
g = gravity(r)
d_vec = (/ V_x, V_y, - g*x/r, - g*y/r /)
end function
3) Later I'll want to do an N-body problem in 3 dimensions.
Each of these requires arrays of different length. I would like to write
the ODE integrator once and re-use it for all these problems. In your
email you suggested that the "dr" function could use assumed shape
arrays as well. I am not sure how I would implement the above functions
using assumed-shape arrays.
> An interface body *MUST* describe the same interface as the procedure
> actually has. A nice thing about module procedures is that you don't
> have to worry about getting the interface body right. Alas, that doesn't
> help you with dummy procedures...
I don't think I fully understand the difference between a module
procedure and a "dummy" procedure.
> The standard has a boring list of what things constitute the interface,
> but for most purposes, you can just forget the list and assume that
> everything must match.
Ok. That's a simple guideline.
As an example, could you show me how you would rewrite this function so
it uses assumed-shape arrays instead of "r(2)"?
pure function dr(t, r)
real(pref), intent(in) :: t, r(2)
real(pref) :: dr(2)
dr = (/ -sin(t) , cos(t) /)
end function
It seems wrong to just write this:
pure function dr(t, r)
real(pref), intent(in) :: t, r(:)
real(pref) :: dr(size(r))
dr = (/ -sin(t) , cos(t) /)
end function
The function only makes sense if size(r) == 2.
> You might have gotten it to compile, but having it actually work
> correctly is a different matter. And if it happens to work on a
> particular compiler, that doesn't guarantee it for others. Having the
> interface body disagree with the actual procedure violates the standard.
I definitely want this to work correctly across compilers. This is for
my masters research and it must be robust.
>> I currently have it at the top of the module. Is that enough? It seems
>> to work...
>
> Yes. You just didn't show that. Context matters.
Yeah, sorry. for completeness, below is a full copy of the ODE module
including comments with a sample program. Thanks for the help.
!**********************************************************************
!
! ode.f90 -- A module with a general-purpose Runge-Kutta solver.
!
!
! NOTE: In the future I hope to add more solvers. For example:
!
! * A Bulirsch–Stoer solver.
! * A High-order Runge-Kutta solver.
! * Adaptive step sizes.
!
! PUBLIC:
!
! ode_rk4 -- Performs one step of Runge-Kutta 4.
!
!
! SAMPLE USAGE:
!
! program circle
! use types
! use ode
! implicit none
!
! real(pref) :: h = 1.0, t = 0, r(2) = (/ 1.0, 0.0 /)
! integer :: i
!
! write (*,*) r
! do i = 1,1000
! call ode_rk4(t, r, dr, h)
! write (*,*) r
! end do
!
! contains
!
! !***************************************
! ! Let "r" be a vector.
! !
! ! | x | d r | - cos(t) |
! ! r = | | ---- = | |
! ! | y | dt | sin(t) |
! !
! !=======================================
! pure function dr(t, r)
! real(pref), intent(in) :: t, r(2)
! real(pref) :: dr(2)
!
! dr = (/ -sin(t) , cos(t) /)
! end function
! end program circle
!
!======================================================================
module ode
use types
implicit none
!
! Access control.
!
private
public ode_rk4
!*****************************************
!!
!!! P U B L I C M E T H O D S
!!
!=========================================
!
! Let Y be a vector or matrix.
!
! Given dY/dt, perform one step of Runge-Kutta.
!
interface ode_rk4
module procedure rk4_val, rk4_vec, rk4_mat
end interface
contains
!*****************************************
!!
!!! R U N G E - K U T T A 4
!!
!=========================================
pure subroutine rk4_val(t, Y, dY, h)
real(pref), intent(inout) :: t, Y
real(pref), intent(in) :: h
real(pref) :: k1, k2, k3, k4
interface
pure function dY(t0, y0)
use types
real(pref), intent(in) :: t0, y0
real(pref) :: dY
end function
end interface
k1 = dY(t, Y)
k2 = dY(t + h/2, Y + k1*h/2)
k3 = dY(t + h/2, Y + k2*h/2)
k4 = dY(t + h , Y + k3*h)
Y = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6
t = t + h
end subroutine
pure subroutine rk4_vec(t, Y, dY, h)
real(pref), intent(inout) :: t, Y(:)
real(pref), intent(in) :: h
real(pref), dimension(size(Y)) :: k1, k2, k3, k4
interface
pure function dY(t0, y0)
use types
real(pref), intent(in) :: t0, y0(:)
real(pref) :: dY(size(y0))
end function
end interface
k1 = dY(t, Y)
k2 = dY(t + h/2, Y + k1*h/2)
k3 = dY(t + h/2, Y + k2*h/2)
k4 = dY(t + h , Y + k3*h)
Y = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6
t = t + h
end subroutine
pure subroutine rk4_mat(t, Y, dY, h)
real(pref), intent(inout) :: t, Y(:,:)
real(pref), intent(in) :: h
real(pref), dimension(size(Y,1),size(Y,2)) :: k1, k2, k3, k4
interface
pure function dY(t0, y0)
use types
real(pref), intent(in) :: t0, y0(:,:)
real(pref) :: dY(size(y0,1),size(y0,2))
end function
end interface
k1 = dY(t, Y)
k2 = dY(t + h/2, Y + k1*h/2)
k3 = dY(t + h/2, Y + k2*h/2)
k4 = dY(t + h , Y + k3*h)
Y = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6
t = t + h
end subroutine
end module
James Van Buskirk
news:ipdqa0$l7b$1...@dont-email.me...
> pure subroutine rk4_vec(t, Y, dY, h)
> real(pref), intent(inout) :: t, Y(:)
> real(pref), intent(in) :: h
> real(pref), dimension(size(Y)) :: k1, k2, k3, k4
> interface
> pure function dY(t0, y0)
> use types
> real(pref), intent(in) :: t0, y0(:)
> real(pref) :: dY(size(y0))
> end function
> end interface
Have you tried changing the above to:
interface
pure function dY(t0, y0)
import
implicit none
real(pref), intent(in) :: t0, y0(size(Y))
real(pref) :: dY(size(y0))
end function
end interface
That way you should be able to use your derivative functions
the way you wrote them in the first place.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
Michel Olagnon
>
>> Ok. I thought "intent(in)" was just something you use to help make sure
>> that you never overwrite a variable.
>
> It is (and it can also help optimizers).
My understanding was that it is to help, but that it is still the
programmer's duty to make sure; some cases could legally not be detected
by the compiler. Has that changed ? Are compilers now in position and
supposed to detect any overwriting of a variable with intent (in) ?
Daniel Carrera
> Have you tried changing the above to:
>
> interface
> pure function dY(t0, y0)
> import
> implicit none
> real(pref), intent(in) :: t0, y0(size(Y))
> real(pref) :: dY(size(y0))
> end function
> end interface
>
> That way you should be able to use your derivative functions
> the way you wrote them in the first place.
I thought Richard said that I wasn't allowed to use Y in the interface.
That said, doing what you suggested seems to work. It compiles fine, and
it even seems to give the right answer. In fact, I think you just fixed
a bug in my program. Yesterday I started getting the wrong answer when I
switched from functions to subroutines.
If someone can confirm that I am allowed to use Y the way you just did
in the interface, I'll be very happy and will move on with my project.
Daniel.
Tobias Burnus
> On 04/29/2011 10:01 AM, James Van Buskirk wrote:
> > interface
> > pure function dY(t0, y0)
> > import
>
> I thought Richard said that I wasn't allowed to use Y in the interface.
Well, note the word "IMPORT" in the code snippet above, which could
also be written as "import :: y".
That's what Richard correctly wrote:
| The one thing that isn't so fine in the interface block is that the
| interface body in it has references to size(Y). Interface bodies do
| *NOT* inherit by host association (unless you use the f2003 IMPORT
| statement, which is, in my view, a hack around a design flaw in
f90).
Tobias
Daniel Carrera
> On Apr 29, 12:42 pm, Daniel Carrera<dan...@gmail.com> wrote:
>> On 04/29/2011 10:01 AM, James Van Buskirk wrote:
>>> interface
>>> pure function dY(t0, y0)
>>> import
>>> real(pref), intent(in) :: t0, y0(size(Y))
>>
>> I thought Richard said that I wasn't allowed to use Y in the interface.
>
> Well, note the word "IMPORT" in the code snippet above, which could
> also be written as "import :: y".
Ooops. Thanks, I missed that entirely. Strangely, "import" works, but
"import :: y" gives me an error:
pure subroutine rk4_vec(t, Y, dY, h)
real(pref), intent(inout) :: t, Y(:)
real(pref), intent(in) :: h
real(pref), dimension(size(Y)) :: k1, k2, k3, k4
interface
pure function dY(t0, y0)
use types
import :: Y
real(pref), intent(in) :: t0, y0(size(Y))
real(pref) :: dY(size(y0))
end function
end interface
k1 = dY(t, Y)
k2 = dY(t + h/2, Y + k1*h/2)
k3 = dY(t + h/2, Y + k2*h/2)
k4 = dY(t + h , Y + k3*h)
Y = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6
t = t + h
end subroutine
The compile error is:
modules/ode.f90:100.27:
import :: Y
1
Error: Cannot IMPORT 'y' from host scoping unit at (1) - does not exist.
Daniel.
fj
No : many ways are possible to modify an INTENT(in) argument. But in that
case, your code is not "standard conforming" and could work with a compiler
but not with another one.
The FORTRAN norm provides the list of cases (rather limited) which must be
detected by a compiler. In practice, good compilers detect many additional
wrong cases not imposed by the norm.
But the list of possible non standard cases is almost illimited. Just a
question of imagination.
--
François Jacq
Richard Maine
> I don't think I fully understand the difference between a module
> procedure and a "dummy" procedure.
That one at least is easy. A module procedure is a procedure that is in
a module. A dummy procedure is a procedure that is a dummy argument.
The one thing that I imagine might confuse you (or others) on this is
the difference between actual and dummy arguments. It applies to any
dummy argument - not just procedures. A dummy argument is distinct from
the actual argument. The dummy and actual arguments get "associated",
but they still count as separate things. In particular, you can say many
things about dummy arguments without even knowing what the actual
arguments are. Heck, you can write and compile a procedure that doesn't
yet have any calls to it, so there aren't any actual arguments anywhere
yet.
A dummy procedure might end up being associated with a module procedure
so that a reference to the dummy procedure ends up actually invoking the
module procedure. But that doesn't mean the dummy procedure is a module
procedure; it is just a dummy procedure.
I'm going to use particular procedure names here because it gets too
confusing talking about things like the procedure that has a procedure
as a dummy argument.
Where this particularly shows up is in specifying the interface. For a
module procedure, just like anything else from a module, you USE the
module and that USE brings in the information necessary for compilation.
For a dummy procedure such as your dY, when compiling your rk4_vec
(which has dY as a dummy) you don't know what the actual procedure(s)
for dY will be. Well, perhaps you might know, but the compiler won't. As
noted above, the calling code that specifies the actual arguments might
not yet even be written. The compiler needs the information now, when
compiling rk4_vec; it affects how the invocation of dY is implemented in
the actual machine code.
So you have to tell the compiler what the interface of dY is. In
f90/f95, you do that with an interface body (or you can have an implicit
interface, but that's not allowed for lots of "interesting" cases,
including yours). The f2003 procedure statement provides another option,
but in any case, since the dummy procedure dY exists only inside of
rk4_vec, you can't just automatically get its interface from somewhere
else.
> As an example, could you show me how you would rewrite this function so
> it uses assumed-shape arrays instead of "r(2)"?
>
> pure function dr(t, r)
> real(pref), intent(in) :: t, r(2)
> real(pref) :: dr(2)
>
> dr = (/ -sin(t) , cos(t) /)
> end function
>
>
> It seems wrong to just write this:
>
> pure function dr(t, r)
> real(pref), intent(in) :: t, r(:)
> real(pref) :: dr(size(r))
>
> dr = (/ -sin(t) , cos(t) /)
> end function
>
>
> The function only makes sense if size(r) == 2.
Well, you have it pretty much correct. I'd write it more like
!-- Presumably pref comes in from host association.
pure function dr(t,r)
real(pref), intent(in) :: t, r(:)
real(pref) :: dr(size(r))
if (size(r) /= 2) then
!-- Oops. Things are amiss. Do what you will.
!-- The PURE atribute is going to get in your way. I'd do
call error_halt('Function dr invoked with wrong size array')
!-- Except that my subroutine error_halt is not pure.
end if
dr = [ -sin(t), cos(t) ]
return
end function dr
Simillar for the other sizes. I almost never write pure procedures
because all but the most trivial of procedures might want error
diagnostics. A recent thread here touched on that. If you don't feel the
need for the error diagnostics, just drop the IF block above completely.
It isn't as though you'd be in any better shape (ouch) with that kind of
invocation error with your explicit-shape version. Compilers typically
will not check that the size matches for explicit-shape arguments. Heck,
some cases of mismatch of size with explicit-shape are actually standard
conforming, so the compiler can't very well diagnose them. They might
not be what you had in mind and might not do anything sensible for your
application, but they are standard conforming.
While I understand in principle the objection to writing code that can
be called with any size arrays even though it is supposed to work only
for a particular size, you are actually safer to write it like my
example above than like your originakl explicit-shape one. My example
will actually diagnose the error (well, once you drop the PURE). Yours
will just compile and run anyway and do something wrong.
I think that explicit-shape dummy arguments are often a bad idea. Not
always, but often. They don't mean what it looks like they mean. In
particular, declaring the dummy argument to have some particular shape
des not mean that the actual argument has to have that shape. Heck, the
actual argument doesn't even have to have the same rank. Instead it
means that the dummy argument will get associated with either part or
all of the actual argument via the complicated mess that is sequence
association; this is regularly a source of confusion.
Dick Hendrickson
> Michel Olagnon wrote:
>
>> On 04/29/2011 12:58 AM, Richard Maine wrote:
>>> Daniel Carrera<dan...@gmail.com> wrote:
>>>
>>>> Ok. I thought "intent(in)" was just something you use to help make sure
>>>> that you never overwrite a variable.
>>>
>>> It is (and it can also help optimizers).
>>
>>
>> My understanding was that it is to help, but that it is still the
>> programmer's duty to make sure; some cases could legally not be detected
>> by the compiler. Has that changed ? Are compilers now in position and
>> supposed to detect any overwriting of a variable with intent (in) ?
>
> No : many ways are possible to modify an INTENT(in) argument. But in that
> case, your code is not "standard conforming" and could work with a compiler
> but not with another one.
>
> The FORTRAN norm provides the list of cases (rather limited) which must be
> detected by a compiler. In practice, good compilers detect many additional
> wrong cases not imposed by the norm.
I think the list of forbidden uses is complete, at least it was intended
to be a complete list. It gives all of the syntax uses that are illegal
in the subprogram with the intent(in) statement. The "only" exception
is passing an intent(in) argument on to a subroutine that either doesn't
have an explicit interface and/or doesn't specify intent for the
argument. You can always beat the compiler by passing things to a
routine it can't inspect. Other than that, I think INTENT(IN) is
bullet-proof.
Dick Hendrickson
fj
No this is not true. Try the following simple counter example :
program t37
real :: a=0.
call test(a)
write(*,*) a
contains
subroutine test(a)
real,intent(in),target :: a
real,pointer :: b
b=>a
b=37
write(*,*) a
end subroutine
end program
With as result with two different compilers :
[lcoul@localhost test]$ ifort t37.f90
[lcoul@localhost test]$ ./a.out
37.00000
[lcoul@localhost test]$ gfortran t37.f90
[lcoul@localhost test]$ ./a.out
37.000000
--
François Jacq
Phillip Helbig---undress to reply
<franco...@irsn.fr> writes:
> No this is not true. Try the following simple counter example :
Here's a counter-example. I'm not saying which compiler is doing the
right thing. (Note: mine is rather ancient, several years old.)
$ type test.f90
program t37
real :: a=0.
call test(a)
write(*,*) a
contains
subroutine test(a)
real,intent(in),target :: a
real,pointer :: b
b=>a
b=37
write(*,*) a
end subroutine
end program
$ fortran test
$ link test
$ run test
37.00000
0.0000000E+00
$ fortran/version
Compaq Fortran V7.5-2630-48C8L
Phillip Helbig---undress to reply
> write(*,*) a
>
> contains
> write(*,*) a
> end subroutine
>
> end program
>
> With as result with two different compilers :
>
> [lcoul@localhost test]$ ifort t37.f90
> [lcoul@localhost test]$ ./a.out
> 37.00000
> [lcoul@localhost test]$ gfortran t37.f90
> [lcoul@localhost test]$ ./a.out
> 37.000000
Since there is a WRITE statement in both the subroutine and the main
program, it would be strange if the above is the total output. Indeed,
most interesting is whether the two WRITE statements produce different
output.
fj
In fact, in my example, I tested alternatively the write statement at the
two positions with exactly the same result using ifort or gfortran. But when
posting, unfortunately, I left the two write statements uncommented.
Anyway, as the code is not standard conforming, a compiler may generate an
executable program writing anything, either 0., 37. or any other value.
The posted version of the test program gives the following results :
[lcoul@localhost test]$ ifort t37.f90
[lcoul@localhost test]$ ./a.out
37.00000
37.00000
[lcoul@localhost test]$ gfortran t37.f90
[lcoul@localhost test]$ ./a.out
37.000000
37.000000
--
François Jacq
Dick Hendrickson
Under INTENT it says
C545 (R517) A nonpointer object with the INTENT (IN) attribute shall not
appear in a variable definition context (16.5.7).
Clearly, "a" is a nonpointer with INTENT(IN) , so
16.5.7 Variable definition context
...
(3) A data-pointer-object or proc-pointer-object in a
pointer-assignment-stmt,
Makes the b=>a non-standard.
I'd be willing to argue that C545 is a constraint and that 16.5.7 merely
provides a definition for "variable definition context." If that's
correct, then the compiler must diagnose the b=>a line. I suppose you
could argue that 16.5.7 is text, not constraint, and the violations
don't to be diagnosed. I don't think that's a good argument since all
of the bullet items in 16.5.7 are syntax things diagnosable at compile
time and many words used in constraints (like INTENT(IN)) are defined
elsewhere in text.
Dick Hendrickson
Tobias Burnus
>> real,intent(in),target :: a
>> real,pointer :: b
>> b=>a
>> b=37
>>
>
> Under INTENT it says
>
> C545 (R517) A nonpointer object with the INTENT (IN) attribute shall not
> appear in a variable definition context (16.5.7).
>
> Clearly, "a" is a nonpointer with INTENT(IN) , so
>
> 16.5.7 Variable definition context
> ...
> (3) A data-pointer-object or proc-pointer-object in a
> pointer-assignment-stmt,
>
> Makes the b=>a non-standard.
Why? "a" is not a "data-pointer-object" but a "data-target"; it's "b"
which is the "data-pointer-object" - and as it is not "intent(in)", its
pointer association status may be changed. Cf.
R733 pointer-assignment-stmt
is data-pointer-object [ (bounds-spec-list) ] => data-target
or ...
Tobias
fj
> I'm not convinced. I think your compilers are broken! (;))
>
> Under INTENT it says
>
> C545 (R517) A nonpointer object with the INTENT (IN) attribute shall not
> appear in a variable definition context (16.5.7).
>
> Clearly, "a" is a nonpointer with INTENT(IN) , so
>
> 16.5.7 Variable definition context
> ...
> (3) A data-pointer-object or proc-pointer-object in a
> pointer-assignment-stmt,
>
> Makes the b=>a non-standard.
No, again I disagree. This statement is correct because the variable "a" is
not the subject of the variable definition context. The statement defines
"b" and not "a". To authorize that, it is just necessary to add the keyword
TARGET when declaring the dummy argument "a".
What is not correct is the statement "b=37.". This one makes the code
invalid but the compiler is not obliged to detect that this statement
indirectly involves a modification of the argument "a" declared
"intent(in)".
No : the compilers are not broken but they are unable to detect many wrong
cases, even very simple ones. And they are many more complicated cases
really impossible to detect.
--
François Jacq
Dick Hendrickson
Dick Hendrickson