add an additional argument for the function
Fatemeh
I have a problem with a subroutine in my program. ( a subroutine from
numerical recipes.)
the total shape of the subroutine is such as :
SUBROUTINE trapzd(func,a,b,s,n)
!This routine computes the nth stage of refinement of an
extended trapezoidal rule. func is
!input as the name of the function to be integrated between
limits a and b, also input. When
!called with n=1, the routine returns as s the crudest
estimate of b
!a f(x)dx. Subsequent
!calls with n=2,3,... (in that sequential order) will improve
the accuracy of s by adding 2n-2
!additional interior points. s should not be modified between
sequential calls
IMPLICIT NONE
REAL(SP), INTENT(IN) :: a,b
REAL(SP), INTENT(INOUT) :: s
INTEGER(I4B), INTENT(IN) :: n
INTERFACE
FUNCTION func(x)
USE nrtype
REAL(SP), DIMENSION(:), INTENT(IN) :: x
REAL(SP), DIMENSION(size(x)) :: func
END FUNCTION func
END INTERFACE
REAL(SP) :: del,fsum
INTEGER(I4B) :: it
if (n == 1) then
s=0.5_sp*(b-a)*sum(func( (/ a,b /) ))
**********************************
else
it=2**(n-2)
del=(b-a)/it !This is the spacing of the points to be added.
fsum=sum(func(arth(a+0.5_sp*del,del,it)))
**********************************
s=0.5_sp*(s+del*fsum) !This replaces s by its refined value.
end if
END SUBROUTINE trapzd
I comment the part
FUNCTION func(x)
USE nrtype
REAL(SP), DIMENSION(:), INTENT(IN) :: x
REAL(SP), DIMENSION(size(x)) :: func
END FUNCTION func
END INTERFACE
and also make a change in the first line such as this :
SUBROUTINE trapzd(func,a,b,s,n) ----> SUBROUTINE trapzd
(a,b,s,n,yy)
and add
REAL(SP), INTENT(IN) :: yy
and after that I define a function after the subroutine which is:
real function func(x,yy)
use nrtype
implicit none
real(SP),intent (in)::x ,yy
REAL(SP) ::bessj0_s,bessj1_s
func=(bessj0_s(x)*bessj1_s(x))/(x*(1+exp(x*yy)))
end function func
my problem is that also I should make some changes in lines with
************** mark
I do such as this:
s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
fsum=sum(func((arth(a+0.5_sp*del,del,it),yy)))
it means that I add an additional argument for the function.
but the compiler makes some errors :
fortcom: Error: aaa.f90, line 161: A constant or named constant is
required in this context.
s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
------------------------------------^
fortcom: Error: aaa.f90, line 161: An INTEGER or REAL data type is
required in this context.
s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
------------------------------------^
fortcom: Error: aaa.f90, line 161: An array-valued argument is
required in this context. [SUM]
s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
--------------------------^
fortcom: Error: aaa.f90, line 165: A constant or named constant is
required in this context. [ARTH]
fsum=sum(func((arth(a+0.5_sp*del,del,it),yy)))
----------------------^
fortcom: Error: aaa.f90, line 165: An INTEGER or REAL data type is
required in this context. [ARTH]
fsum=sum(func((arth(a+0.5_sp*del,del,it),yy)))
----------------------^
fortcom: Error: aaa.f90, line 165: An array-valued argument is
required in this context. [SUM]
fsum=sum(func((arth(a+0.5_sp*del,del,it),yy)))
----------------^
fortcom: Error: aaa.f90, line 216: Syntax error, found '.' when
expecting one of: ( * , <IDENTIFIER> <CHAR_CON_KIND_PARAM>
<CHAR_NAM_KIND_PARAM> <CHARACTER_CONSTANT> ...
call nrerror(.rtsec: exceed maximum iterations.)
---------------------^
Can anyone guide me what should I do ?
because I have to define a function with 2 argument so I need to enter
" yy " in my subroutine as an input.
I would appreciate your help.
Best,
Fatemeh
Arjen Markus
> Dear all ;
>
> I have a problem with a subroutine in my program. ( a subroutine from
> numerical recipes.)
> the total shape of the subroutine is such as :
>
> SUBROUTINE trapzd(func,a,b,s,n)
> extended trapezoidal rule. func is
The problem you sketch is quite a tricky one. There are no really
satisfactory
solutions, but if you do not use multiple threads, you can do
something like:
module functions
use nrtype
implicit none
real :: yy
contains
subroutine set_yy( new_yy )
real :: new_yy
yy = new_yy
end subroutine
real function func(x)
real(SP), intent(in) :: x
func=(bessj0_s(x)*bessj1_s(x))/(x*(1+exp(x*yy)))
end function func
end module functions
And use the trapzd routine like:
use functions
...
call set_yy( yy )
call trapzd( func, a, ... )
The problem is that you can not expand the function's interface at
will.
In Fortran 2003 there are more possibilities, but they still require
cooperation of the trapzd routine.
An alternative could be to customise the trapzd, so that it takes an
extra argument and the function it expects takes that argument as
well.
Regards,
Arjen
e p chandler
Thanks for quickly outlining a possible solution for the OP.
Several problems remain.
1. Func takes a vector as an argument.
2. Function arth is undefined. I suppose it returns a vector
containing values in an arithmetic progression.
3. IMO "Numerical Recipes" is a bad choice. The addition of Fortran 90
to the mix only serves to confuse the student more. In some respects,
the Fortran 77 version would have been a better choice, but not much.
Please excuse my comments on style and logistics.
1. It helps to know what operating system, compiler and compiler
version the OP is using. One strategy for solving these types of
problems is to just compile program units separately. It is difficult
to know which options to suggest without knowing which compiler, etc.
For example, I would look at the original subroutine trapzd and
compile and then stop with
g95 -c foo.f90
2. We still have not seen any main program listed. Often the calling
program and its declarations are vital to understanding a problem.
3. This discussion has been scattered over a number of different
postings in this newsgroup. It is far better to start a topic as one
posting and see it progress.
4. If the OP has Fortran 90+, then he should use free form to full
advantage. There is no reason to indent each line of source code. Save
indentation for demonstrating block structure.
5. I can not tell for sure, but the 8 leading spaces on each line lead
me to wonder if the OP has used TAB characters in his source code.
Please, please do not do this!
6. ******* is a Fortran 77 style comment. Replace the first * with !
if you must use this style of comment.
7. The OP needs to find a good Fortran 95+ textbook and RTFM (read the
Fortran Manual).
8. I also second all comments suggesting that the OP is trying to do
too much at one time. Build a program from individual pieces that
work. Test each piece first. Then refine your program. This is a
strategy for success.
----- e
James Van Buskirk
news:9a4b1822-baef-407f...@l39g2000yqn.googlegroups.com...
> I do such as this:
> s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
> fsum=sum(func((arth(a+0.5_sp*del,del,it),yy)))
> it means that I add an additional argument for the function.
> but the compiler makes some errors :
> fortcom: Error: aaa.f90, line 161: A constant or named constant is
> required in this context.
> s=0.5_sp*(b-a)*sum(func(( (/ a,b /),yy) ))
> ------------------------------------^
Your program seems to have hit a blind spot: the syntax for a COMPLEX
literal is:
(a,b)
where a and b are both numeric literals. a and b can't be variables
or even named constants (well, named constants are allowed in f2003,
see N1601.pdf, section 4.4.3).
What's happening to you is that you're putting extra parentheses
around a pair of arguments. If you put extra parentheses around a
single argument it has the effect of creating an expression out of
the argument and passing (probably a reference to) the expression.
Fine if that's what you want to do or if the single argument is
already an expression and not a variable. But if the argument is
a variable that gets changed within the procedure, then it's wrong
and likely won't result in the variable being changed on return or
can cause a crash within the called procedure.
With a pair or arguments, the compiler thinks you are trying to
write out a COMPLEX literal, but since you aren't following the
rules for COMPLEX literals it's giving you all sorts of errors.
I don't know what ifort would do if you enclosed three or more
arguments with extra parentheses. Maybe you can do that and report
the results here. To fix your problem though, remove the extra
parentheses:
s=0.5_sp*(b-a)*sum(func( (/ a,b /),yy) )
fsum=sum(func(arth(a+0.5_sp*del,del,it),yy))
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
James Van Buskirk
news:a03d3f3c-3830-4da0...@v19g2000yqn.googlegroups.com...
> The problem is that you can not expand the function's interface at
> will.
> In Fortran 2003 there are more possibilities, but they still require
> cooperation of the trapzd routine.
In Fortran 2008 you can pass internal procedures as actual arguments
and this would solve the problem without modification of trapzd in a
thread-safe way. The OP is using ifort which has supported this since
it became the new dvf. I have posted some examples of this usage
before:
http://groups.google.com/group/comp.lang.fortran/msg/a86c0c018bf16607?hl=en
> An alternative could be to customise the trapzd, so that it takes an
> extra argument and the function it expects takes that argument as
> well.
This is the approach the OP has taken, but was frustrated by syntax
errors.
e p chandler
[snip]
Even if the OP does manage to use this technique or if he choses to
add the extra argument YY to the argument lists of trapzd AND func, he
will still run into a numerical brick wall at high speed.
From a different (one of many) thread in this newsgroup it appears
that he wants to integrate from zero to infinity.
But his integrand is
func=(bessj0_s(x)*bessj1_s(x))/(x*(1+exp(x*yy)))
which is singular at x == 0 and ill behaved thereabout.
I was able, passing yy as an extra argument, to write a program that
does compile and run using g95. That required writing my own
replacements for nrtypes and for a module which contains the function
arth. In addition, the bessel functions used are those supplied with
NR. I replaced them with the ones supplied as extensions with g95. [I
hope but do not know if they are equivalent.] Also, I have to guess
what the program that calls trapzd looks like.
Here is my program:
--- start text ---
! my own modules to replace those of NR
MODULE nrtype
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: I4B = KIND(1)
END MODULE nrtype
MODULE nrutil
USE nrtype
IMPLICIT NONE
CONTAINS
FUNCTION arth(start,delta,intervals)
REAL(SP), INTENT(IN) :: start, delta
INTEGER(I4B), INTENT(IN) :: intervals
REAL(SP), DIMENSION(intervals) :: arth
INTEGER(I4B) :: i
arth = (/ (start + (i-1) * delta, i=1,intervals) /)
END FUNCTION arth
END MODULE nrutil
SUBROUTINE trapzd(func,a,b,s,n,yy)
USE nrtype
USE nrutil
IMPLICIT NONE
REAL(SP), INTENT(IN) :: a,b,yy
REAL(SP), INTENT(INOUT) :: s
INTEGER(I4B), INTENT(IN) :: n
INTERFACE
FUNCTION func(x,yy)
USE nrtype
REAL(SP), DIMENSION(:), INTENT(IN) :: x
REAL(SP), DIMENSION(size(x)) :: func
REAL(SP), INTENT(IN) :: yy
END FUNCTION func
END INTERFACE
REAL(SP) :: del,fsum
INTEGER(I4B) :: it
if (n == 1) then
s=0.5_sp*(b-a)*sum(func( (/ a,b /), yy ))
else
it=2**(n-2)
del=(b-a)/it
fsum=sum(func(arth(a+0.5_sp*del,del,it), yy ))
s=0.5_sp*(s+del*fsum)
end if
END SUBROUTINE trapzd
FUNCTION func(x,yy)
USE nrtype
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: x
REAL(SP), DIMENSION(size(x)) :: func
REAL(SP), INTENT(IN) :: yy
! original func=(bessj0_s(x)*bessj1_s(x))/(x*(1+exp(x*yy)))
! avoid x = 0, where the function is singular
! change bessel functions to extensions supplied with g95
func=(besj0(x)*besj1(x)) / ( x * (1+exp(x*yy) ) )
END FUNCTION func
PROGRAM nr90
USE nrtype
IMPLICIT NONE
REAL(SP) :: a,b,s0,s,eps,yy
INTEGER(I4B) :: n
REAL(SP) :: func
EXTERNAL func
INTEGER(I4B), PARAMETER :: NMAX = 10
eps = 1e-4
a = .05 ! singularity at 0
! is this valid ??? who knows ???
b = 10.0
yy = 1
n = 1
call trapzd(func,a,b,s,n,yy)
do
n = n + 1
if (n > NMAX) then
print *,'failed to converge'
print *,'s=',s
print *,'at n=',NMAX
exit ! do
end if
s0=s
call trapzd(func,a,b,s,n,yy)
print *,n,s,s0
if (abs(s-s0) < eps*abs(s0)) then
print *,'yy=',yy,'s=',s,'n=',n
exit ! do
end if
end do
END PROGRAM nr90
--- end text ---
Here is the output:
2 0.60612494 1.2115268
3 0.30071536 0.60612494
4 0.22001435 0.30071536
5 0.21307474 0.22001435
6 0.19193985 0.21307474
7 0.15673803 0.19193985
8 0.12085123 0.15673803
9 0.091628835 0.12085123
10 0.070696205 0.091628835
failed to converge
s= 0.070696205
at n= 10
--- end text ---
Notes to the OP.
1. Func must be an external procedure to be passed as an argument. So
it can't be put in a module. Thus it needs an interface body written
out and that interface body must be repeated in the function itself.
2. Func needs to be decared as REAL(SP) and EXTERNAL in the calling
program. External because it is an external function passed to a
procedure as an argument.
I also took a sneak peek at how arth was defined in NR. IIRC it
calculates the entries in its output vector by successively adding DEL
to its start value. It could also be that my function arth is
mismatched with what trapzd expects. ICK!
So there you go. S*** with icing, Fortran 90 style.
--- e
Richard Maine
> 1. Func must be an external procedure to be passed as an argument. So
> it can't be put in a module. Thus it needs an interface body written
> out and that interface body must be repeated in the function itself.
I haven't paid a lot of attention to most of this thread, but the above
caught my eye. That's just not true. There is nothing at all wrong with
passing a module procedure as an actual argument. In fact, I recommend
making pretty much all procedures module procedures. I'm not sure what
would lead you to the above incorrect conclusion. My guesses would be
1. You might be confusing module procedures with internal ones. Although
they both follow a CONTAINS statement, they are not at all the same
thing. Module procedures don't have the issues that have prompted
disallowance of internal procedures as actual arguments.
2. You might be mislead by the use of the EXTERNAL statement to declare
dummy procedures. That particular bit of syntax is horribly misleading
for historical reasons. For a dummy procedure, it just declares that the
dummy is a procedure. It does *NOT* have anything to do with declaring
that the actual argument has to be an external procedure - it doesn't
have to be.
If it isn't one of those confusions, I'm at loss as to what it might be.
But in any case, passing a module procedure as an actual argument is
perfectly fine, normal, and even recommended.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
James Van Buskirk
news:cbaad509-9ba3-4894...@z9g2000yqi.googlegroups.com...
> Even if the OP does manage to use this technique or if he choses to
> add the extra argument YY to the argument lists of trapzd AND func, he
> will still run into a numerical brick wall at high speed.
> From a different (one of many) thread in this newsgroup it appears
> that he wants to integrate from zero to infinity.
> But his integrand is
> func=(bessj0_s(x)*bessj1_s(x))/(x*(1+exp(x*yy)))
> which is singular at x == 0 and ill behaved thereabout.
The singularity at x = 0 is removable since
lim_{x -> 0) BESSEL_J0(x)/x = 0.5, so if this singularity is removed
the integrand will be well-behaved.
Also I think trapzd is intended to be called by qromb, which is a
much more reasonable numerical quadrature method.
James Van Buskirk
news:gpi7tc$7kj$1...@news.motzarella.org...
> lim_{x -> 0) BESSEL_J0(x)/x = 0.5, so if this singularity is removed
> the integrand will be well-behaved.
<gack> I meant BESSEL_J1(x)/x, of course!
glen herrmannsfeldt
(snip)
> What's happening to you is that you're putting extra parentheses
> around a pair of arguments. If you put extra parentheses around a
> single argument it has the effect of creating an expression out of
> the argument and passing (probably a reference to) the expression.
PL/I does that, but as far as I know Fortran doesn't, or at
least doesn't guarantee it. Parenthesis around a variable
may or may not pass a temporary copy, and you aren't allowed
to modify the copy if it does.
-- glen
Richard Maine
Fortran doesn't talk about things at the level of implementation of
making copies. But you *ARE* allowed to pass the same variable as a
separate argument, which then is allowed to be modified. Yes, the
standard requires that. Compilers have been known to get it wrong, but
that's called a bug. The "obvious" method of implementation is to make a
copy, perhaps unless your optimizer can figure out that it isn't needed.
It does pretty much have to act like a copy.
e p chandler
OK. You are correct about both 1 and 2. The last time I looked at this
miserable trapezoidal integration code from NR (the Fortran 77
version), I *did* put my integrand inside a module. Yes, I did need to
declare its type.
I was caught by CONTAINS and EXTERNAL as you explain. As you say,
experience comes from bad judgment. [smile]
With a bit of work I was able to get a version of the program with
FUNC in a module to run, but I had to retain the interface block from
the original program to do so, even though FUNC is now a module
procedure.
So I'm still confused about what TRAPZD "inherits" when a module is
USED in its parent program.
---- e
e p chandler
> "James Van Buskirk" <not_va...@comcast.net> wrote in messagenews:gpi7tc$7kj$1...@news.motzarella.org...
>
> > lim_{x -> 0) BESSEL_J0(x)/x = 0.5, so if this singularity is removed
> > the integrand will be well-behaved.
>
> <gack> I meant BESSEL_J1(x)/x, of course!
This is generally over my head, BUT
If I want to deal with the singularity by evaluating X,
say IF (X < some_epsilon) then
func = something
else
func = something else
HOWEVER, func has a vector argument here which the program is trying
to evaluate with a single statement.
1. I should probably RTFM for the solution. func = WHERE ...
ELSEWHERE .... ?
2. My real point is that by being "cute" and playing with fancy
features of Fortran 90, the authors of NR have made the program more
complicated, the interface more complicated, and the issue of dealing
with singularities more complicated, in addition to all of this
confusing the student (both myself AND the OP)!
--- e
Richard Maine
> So I'm still confused about what TRAPZD "inherits" when a module is
> USED in its parent program.
The same kinds of things it would inherit if were used as an internal
procedure of a program or subprogram - most identifiers in the host.
Names of types, variables, procedures etc. That part has nothing to do
with internal versus module procedures.
If anything, it is simpler for module procedures because there can never
be simultaneous multiple instances of things in the module. That's a
complication of internal procedures as actual arguments. If it is an
internal procedure of a recursive procedure, you have to worry about
which instance of the recursive procedure it is internal to.
James Van Buskirk
news:1iwm5rn.ejyvh5hrvji2N%nos...@see.signature...
>e p chandler <ep...@juno.com> wrote:
>> So I'm still confused about what TRAPZD "inherits" when a module is
>> USED in its parent program.
> The same kinds of things it would inherit if were used as an internal
> procedure of a program or subprogram - most identifiers in the host.
> Names of types, variables, procedures etc. That part has nothing to do
> with internal versus module procedures.
> If anything, it is simpler for module procedures because there can never
> be simultaneous multiple instances of things in the module. That's a
> complication of internal procedures as actual arguments. If it is an
> internal procedure of a recursive procedure, you have to worry about
> which instance of the recursive procedure it is internal to.
One thing that wasn't pointed out in this subthread is that with
f2003 you don't need to write out interface blocks if there is a
module procedure or abstract interface available to clone the dummy
procedure from:
C:\gfortran\clf\proc_iface>type proc_iface.f90
module funcs1
implicit none
contains
function fun(a,b,n)
real a, b
integer n
real fun(n)
integer i
if(n > 0) fun = (/a,(a+(b-a)*i/(n-1),i=1,n-1)/)
end function fun
end module funcs1
module funcs2
implicit none
contains
subroutine sub(f,a,b,n)
use funcs1
procedure(fun) f
real a, b
integer n
write(*,*) f(a,b,n)
end subroutine sub
end module funcs2
program prog
use funcs1
use funcs2
implicit none
real a, b
integer n
a = 7
b = 11
n = 9
call sub(fun,a,b,n)
end program prog
C:\gfortran\clf\proc_iface>gfortran proc_iface.f90 -oproc_iface
C:\gfortran\clf\proc_iface>proc_iface
7.0000000 7.5000000 8.0000000 8.5000000 9.0000000
9.5000000 10.000000 10.500000 11.000000
e p chandler
> "Richard Maine" <nos...@see.signature> wrote in message
>
> news:1iwm5rn.ejyvh5hrvji2N%nos...@see.signature...
>
Very interesting! That's exactly the sub-problem I was interested in -
how to avoid writing an interface block for a procedure passed as an
argument.
Thanks to both you and Richard for clarifying things.
-- elliot
Ron Shepard
<3d2d22e0-b6bf-48ac...@v15g2000yqn.googlegroups.com>,
e p chandler <ep...@juno.com> wrote:
> With a bit of work I was able to get a version of the program with
> FUNC in a module to run, but I had to retain the interface block from
> the original program to do so, even though FUNC is now a module
> procedure.
Most of fortran seems to be designed so that multiple definitions of
explicit interfaces can be avoided. In many situations it isn't
even possible to have multiple definitions because the language
prevents it, which is nice because it prevents the programmer from
shooting himself in the foot. But one exception to this is the
interface block for a dummy procedure. In this case, it seems that
the programmer MUST specify it redundantly. Furthermore, some
compilers do not check to ensure that the interface of the actual
procedure matches that of the dummy procedure. Maybe in some cases
this is not possible, but in the simple cases it is. Some compilers
do check, which is a nice implementation feature.
$.02 -Ron Shepard
Richard Maine
> Most of fortran seems to be designed so that multiple definitions of
> explicit interfaces can be avoided. In many situations it isn't
> even possible to have multiple definitions because the language
> prevents it, which is nice because it prevents the programmer from
> shooting himself in the foot. But one exception to this is the
> interface block for a dummy procedure. In this case, it seems that
> the programmer MUST specify it redundantly.
See James' comment about the new f2003 feature for that.
glen herrmannsfeldt
> If anything, it is simpler for module procedures because there can never
> be simultaneous multiple instances of things in the module. That's a
> complication of internal procedures as actual arguments. If it is an
> internal procedure of a recursive procedure, you have to worry about
> which instance of the recursive procedure it is internal to.
So they could have had the rule that only internal procedures of
recursive procedures not be allowed as actual arguments.
(Not that Fortran doesn't already have enough special case
rules, but it would seem that it could have been done that way.)
-- glen
glen herrmannsfeldt
> Fortran doesn't talk about things at the level of implementation of
> making copies. But you *ARE* allowed to pass the same variable as a
> separate argument, which then is allowed to be modified. Yes, the
> standard requires that. Compilers have been known to get it wrong, but
> that's called a bug. The "obvious" method of implementation is to make a
> copy, perhaps unless your optimizer can figure out that it isn't needed.
> It does pretty much have to act like a copy.
To be specific, the case of
call sub((x))
is (x) an expression or a variable?
In PL/I it is an expression, and a temporary variable is passed
as it would for any expression. Also, PL/I allows the called
routine to modify the temporary, and guarantees not to
modify x.
It is my understanding that in that case, or even in
call sub(x+0)
Fortran allows, but does not require, a copy. That the
called routine should not modify the dummy argument, as it
may be an unmodifiable copy, or x.
-- glen
e p chandler
argument
On Mar 15, 12:46 am, e p chandler <e...@juno.com> wrote:
> Here is the output:
>
> 2 0.60612494 1.2115268
> 3 0.30071536 0.60612494
> 4 0.22001435 0.30071536
> 5 0.21307474 0.22001435
> 6 0.19193985 0.21307474
> 7 0.15673803 0.19193985
> 8 0.12085123 0.15673803
> 9 0.091628835 0.12085123
> 10 0.070696205 0.091628835
> failed to converge
> s= 0.070696205
> at n= 10
The program runs properly when compiled with gfortran but not g95, but
it took me a long time to find out why.
Here is my integrand function
func=(besj0(x)*besj1(x)) / ( x * (1+exp(x*yy) ) )
But X is not a scalar. X is a vector, thanks to the design of NR for
F90.
Under gfortran, its besj0 and besj1 are both elemental, but under g95,
they are not.
Here's a short demonstration program:
C:\Users\epc\temp>type z.f90
real x(2)
x(1)=1
x(2)=2
print *,x,besj0(x(1)),besj0(x(2))
print *,x,besj0(x)
end
C:\Users\epc\temp>g95 z.f90
C:\Users\epc\temp>a
1. 2. 0.7651977 0.22389078
1. 2. 0.7651977
C:\Users\epc\temp>gfortran z.f90
C:\Users\epc\temp>a
1.0000000 2.0000000 0.76519769 0.22389078
1.0000000 2.0000000 0.76519769 0.22389078
In fact, under g95, if an array is passed to besjN, elements 2 and on
are not defined.
GRrrrrrrrrrrrrrr. I won't tell anyone how long it took me to
systematicly whittle down my test program to find the "bug".
[So long IT career!]
---- e
Richard Maine
> Richard Maine <nos...@see.signature> wrote:
>
> > If anything, it is simpler for module procedures because there can never
> > be simultaneous multiple instances of things in the module. That's a
> > complication of internal procedures as actual arguments. If it is an
> > internal procedure of a recursive procedure, you have to worry about
> > which instance of the recursive procedure it is internal to.
>
> So they could have had the rule that only internal procedures of
> recursive procedures not be allowed as actual arguments.
They could have. But they didn't.
Richard Maine
> call sub((x))
>
> is (x) an expression or a variable?
That's trivially easy. Of course it is an expression. There is no
variable syntax like that.
> It is my understanding that in that case, or even in
>
> call sub(x+0)
>
> Fortran allows, but does not require, a copy. That the
> called routine should not modify the dummy argument, as it
> may be an unmodifiable copy, or x.
Fortran doesn't talk about copying at all, so that's pretty much
irrelevant. And that example doesn't have any of the features that make
the question interesting anyway. See my prior comments about the
interactions with *OTHER* arguments. There have actually been extensive
discussions of it in this newsgroup in the past. Much more interesting
examples are things like
call sub(x,(x))
in which case the first argument *MAY* be modified, but such
modifications don't affect the second one.
Ron Shepard
glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
> To be specific, the case of
>
> call sub((x))
>
> is (x) an expression or a variable?
It is an expression, but this isn't one of the interesting cases.
Consider a call like
call sub( (x), (x), x )
In this case, the first two arguments are expressions and the last
one is a variable. The compiler is allowed to evaluate the
expression only once, store that value in a temporary location, and
then pass that same location address as the first two arguments.
Unless there is an explicit interface that places additional
requirements on the arguments, the subroutine can modify its third
argument, but it cannot modify the first two.
$.02 -Ron Shepard
glen herrmannsfeldt
(snip)
> call sub( (x), (x), x )
> In this case, the first two arguments are expressions and the last
> one is a variable. The compiler is allowed to evaluate the
> expression only once, store that value in a temporary location, and
> then pass that same location address as the first two arguments.
> Unless there is an explicit interface that places additional
> requirements on the arguments, the subroutine can modify its third
> argument, but it cannot modify the first two.
OK, but in those cases the temporary is required, yes?
In the:
call sub((x))
case, as I understand it, no temporary is required and the
called routine is not allowed to modify the dummy argument.
-- glen
Richard Maine
> Ron Shepard <ron-s...@nospam.comcast.net> wrote:
> (snip)
>
> > call sub( (x), (x), x )
> OK, but in those cases the temporary is required, yes?
No. And I'm reasonably confident you know the standard well enough to
know why that's the answer, but...
The standard doesn't talk about things like temporaries being required
or not. (Well, the VALUE attribute might say something related to that,
but it's a bit different in that way).
I realize that you like to think of things in terms of the most likely
implementations. And that can be fine to a certain extent. But that is
*NOT* what the standard requires. So as soon as you start asking whether
something is required at that level, the answer is pretty much always
"no". It might be the most likely implementation, but no, it is not
required and that is not the same thing. It is one of these things where
you don't even need to hear the rest of the question. If the question
starts out something like "Does the standard require a temporary copy
when..." then the answer is "no".
Even if you couldn't think of any other way to meet the standard, "I
can't think of any other way" does not constitute a requirement of the
standard.
One can certainly think of other ways here. And I don't really care
about the likely objections that these wouldn't apply to all cases or
might be unlikely for an actual compiler to do. That's quite irrelevant
to whether this is a requirement of the standard.
One "obvious" way is for the compiler's optimizer to notice that, even
though the standard allows the programer to modify the 3rd dummy
argument, the code doesn't actually do that and therefore can be
optimized.
glen herrmannsfeldt
> glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
>> Ron Shepard <ron-s...@nospam.comcast.net> wrote:
>> (snip)
>
>> > call sub( (x), (x), x )
>> OK, but in those cases the temporary is required, yes?
> No. And I'm reasonably confident you know the standard well enough to
> know why that's the answer, but...
> The standard doesn't talk about things like temporaries being required
> or not. (Well, the VALUE attribute might say something related to that,
> but it's a bit different in that way).
OK, then I will make it more specific. If the called routine
modifies the third dummy argument does the standard not
allow the first dummy argument to change?
As you say, the compiler might notice that you don't do that,
but for a separate compilation it couldn't notice that.
-- glen
nm...@cam.ac.uk
>modifies the third dummy argument does the standard not
>allow the first dummy argument to change?
>
>As you say, the compiler might notice that you don't do that,
>but for a separate compilation it couldn't notice that.
It is undefined behaviour if it does. With some exceptions.
Regards,
Nick Maclaren.
James Van Buskirk
news:gpld30$oeu$1...@soup.linux.pwf.cam.ac.uk...
I'm not sure what the first sentence above is supposed to mean.
Does 'if it does' mean if sub defines its first dummy argument?
If so, then the statement is incorrect because the first dummy
argument is not definable. That isn't undefined behavior, that's
a nonstandard program.
If 'if it does' means that if sub defines its third dummy argument,
the statement is again incorrect because the first actual arguments
are given by expressions so that by the time sub can see them they
are distinct entities from the variable third dummy. Defining
that third dummy in that case is fine and behavior is also well-
defined: the first and second dummies must go unaffected.
N1601.pdf, section 6:
"R601 variable is designator"
"R602 variable-name is name"
"R603 designator is object-name
or array-element
or array-section
or structure-component
or substring"
Section 7.1.1.1:
"R701 primary is constant
or designator
...
or (expr)"
So x has to be a variable unless it's a named constant or a
procedure name, and (x) is an expression, not a variable.
Modifying a variable doesn't do anything to the value of an
expression that has already been evaluated.
I don't see where there is anything complicated in any of this.
Ron Shepard
glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
> >> > call sub( (x), (x), x )
>
> >> OK, but in those cases the temporary is required, yes?
>
> > No. And I'm reasonably confident you know the standard well enough to
> > know why that's the answer, but...
>
> > The standard doesn't talk about things like temporaries being required
> > or not. (Well, the VALUE attribute might say something related to that,
> > but it's a bit different in that way).
>
> OK, then I will make it more specific. If the called routine
> modifies the third dummy argument does the standard not
> allow the first dummy argument to change?
If I understand your question correctly, then I think the standard
requires that the first (and second) dummy argument does not change
when the third dummy argument value is changed.
If you use a LOC() function, or something similar, then the
addresses of the first two dummy arguments could be the same, but
the last one would be different. The standard does not require this
behavior, of course, but if you want to understand what is
happening, that is one way to figure it out.
On the other hand, if the third dummy argument is declared as
INTENT(IN), then the compiler could use the same address for all
three arguments. It is not required to do so, but I think it could.
It is this common address thing that I think is the root cause for
the prohibition in fortran against modifying dummy arguments that
are associated with expressions. In a sequence of subprogram calls
call sub1( 1, 1, 1, x )
call sub2( 1, 1, 1, y )
that common integer constant "1" might end up having the same
address for the corresponding dummy arguments. It it were allowed,
then changing one of those values would then incorrectly change the
others, causing all kinds of difficult to find problems in the code.
When expressions were allowed in argument lists, they inherited
these same prohibitions against being changed through argument
association for, I think, the same general reasons.
$.02 -Ron Shepard