Partial derivative
Allamarein
REAL*8 FUNCTION F(x,y,z)
IMPLICIT NONE
REAL*8, INTENT(IN) :: x,y,z
f = log(x)+y**2+z
END FUNCTION
I would compute automatically the three partial derivatives of the
above function by means of a subroutine.
How can I do?
Paul van Delst
a subroutine that implements those equations. (I'm not being facetious, that is what I would do).
cheers,
paulv
p.s. you might want to add some input checking on the value of the x dummy argument too.
ralf.schaa
> f = log(x)+y**2+z
> I would compute automatically the three partial derivatives of the
> above function by means of a subroutine.
> How can I do?
real function partials_f(x,y)
real, intent(in) :: x,y
partials_f = [1/x, 2*y, 1 ]
end function partials_f
e p chandler
"Allamarein" wrote
Please clarify. Do you want numerical routines which approximate the partial
derivatives for ANY F(x,y,z) OR do you want the partials of THIS particular
function?
ralf.schaa
ok, i think i forgot: dimension(3), partials_f
Gib Bogle
Homework?
glen herrmannsfeldt
Some do it numerically, though the results aren't quite as
good as one might hope.
There are ways to do symbolic derivatives, including calculators
such as the HP-28C and TI-92. Also, Mathematica can do it and
also output in FortranForm[].
If the function is known at compile time, it is common to use
symbolic derivatives.
-- glen
Allamarein
>derivatives for ANY F(x,y,z) OR do you want the partials of THIS particular
>function?
No. It was just an example.
It should derive ANY F(x,y,z) by a numerical method.
I would not pre-process my function by Matlab Symbolic Toolbox or
similar.
I am looking for a genuine numerical method to find partial
derivatives for ANY F.
e p chandler
"Allamarein" <matteo.d...@gmail.com> wrote in message
news:a2567bd7-eb33-41b8...@d12g2000vbz.googlegroups.com...
First of all, you have to deal with the arguments to your function and its
returned value. Full generality would be very difficult. But let's suppose
that your function, for the sake of argument, does take default reals x, y
and z and returns a default real.
Of course you could start with some simple numerical approximation that uses
differences between tabulated or calculated values. Or you could look at it
as an interpolation (or extrapolation) problem.
But IMO, just finding a numerical method that calculates a partial
derivative is premature. The answer is related to
1) What problem are you trying to solve?
2) What general method are you trying to use?
The chances are that someone else has already written programs or developed
methods in your area of interest. For example, you can find any number of
video tutorials on the internet which discuss the diffusion of heat.
Likewise, there are well known techniques for dealing with partial
differential equations of many types and relating to many fields.
Just my two cents worth.
-- e
Arjen Markus
You could expand my automatic differentiation module at http://flibs.sf.net.
It currently only works for functions of one variable only, but that
is a
(mathematical and programmatical) detail.
Regards,
Arjen
Allamarein
> You could expand my automatic differentiation module athttp://flibs.sf.net.
> It currently only works for functions of one variable only, but that
> is a
> (mathematical and programmatical) detail.
>
> Regards,
>
> Arjen
You are right.
Effectively I could use the subroutine for functions R->R, enabling a
variable at each time.
Below you can find a draft of the supposed code.
I supposed that the function 'auto_derive' handling function R->R, as
Arjen suggested: given a function and a point, it turns out the value
of the derivative in that point.
'Pot' is the potential function R3->R, which I need its gradient.
'PD' is the subroutine that computes the gradient
I have dropped some formal details (e.g. declarations of some
variables, etc..)
In your opinion a code like this could approach well my problem?
program demo
use variables
real*8, dimension(1:3) :: grad
x = 2.0
y = 3.0
z = 1.0
!results
call PD(x,y,z,grad)
end
real*8 function pot(x,y,z)
fun = x**2+log(y)+z
end function
subroutine PD(x,y,z,yout)
yout(1) = auto_derive(deriveX(pot,x))
yout(2) = auto_derive(deriveY(pot,y))
yout(3) = auto_derive(deriveZ(pot,z))
end
real*8 function deriveX(fun,x0)
use variables, only :: y,z
deriveX = fun(x0, y, z)
end function
real*8 function deriveY(fun,y0)
use variables, only :: x,z
deriveY = fun(x, y0, z)
end function
real*8 function deriveZ(fun,z0)
use variables, only :: x,y
deriveZ = fun(x, y, z0)
end function
module variables
real*8 :: x,y,z
end module
Arjen Markus
The technique I use is based on:
- Mathematically exact expressions for the derivatives of
common functions and composition of mathematical functions
- Overloading of functions
That way the compiler handles all the difficult details.
What you could do is provide three versions of the function:
function pot_dx(x,y,z) result(fun)
real :: y, z
type(auto_deriv) :: x, fun
fun = x**2+log(y)+z
end function
and similar for the derivatives in y and z direction.
Then:
type(auto_deriv) :: grad, var
real :: gradient
var%v = x0
var%dv = 1.0
grad(1) = pot_dx(var, y0, z0)
var%v = y0
grad(2) = pot_dy(x0, var, z0)
var%v = z0
grad(3) = pot_dz(x0, y0, var)
gradient = grad%dv
Or something along those lines
Regards,
Arjen
fj
Notice that I posted, several months ago (one year?), a FORTRAN module
var_dep.f90 working in a similar way as the Arjen's tool (but still more
genaral) : grouping together a value and its partial derivatives in a
derived type variable. All basic operators and mathematical functions are
implemented betwwen vardep/vardep vardep/scalar and scalar/vardep
combinations. So a TYPE(vardep) variable may be used as easily as a scalar
value.
My post was pointing about performance issues between fixed size arrays and
allocatables (on needs to retain, in principle, an unknown number of partial
derivatives).
The module is presently in service in an actual scientific code but,
unfortunately, only the fixed size version has finally been selected (much
faster)...
--
François Jacq
Tobias Burnus
> My post was pointing about performance issues between fixed size arrays and
> allocatables (on needs to retain, in principle, an unknown number of partial
> derivatives).
> The module is presently in service in an actual scientific code but,
> unfortunately, only the fixed size version has finally been selected (much
> faster)...
In principle there is no reason why a fixed size version should be
(much) faster than an allocatable-array version.
It is not completely clear how the terms translate into standard speech
(Fixed size = explicit size with compile-time known size? Allocatable
array = the dummy argument is allocatable?)
I actually see only two (three) potential performance issues (speaking
of a common implementation, not of all possible implementations the
standard permits):
a) Passing an assumed-shape array to an explicit-size dummy
b) Using an assumed-shaped array in loops or on whole-array operations
c) TARGET and POINTER ...
For (a), the compiler needs to do a copy in/copy out unless it know that
the actual argument is contiguous, which it usually does not for
assumed-shape arrays.
For (b) the issue is that the compiler needs to assume that the array is
not contiguous, thus it needs to take the strides into account.
Depending on the code, either (a) or (b) causes the greater slow down.
However, if the argument is allocatable, it is contiguous. Thus, both
issues shouldn't occur and there shouldn't be any larger performance
difference for your program.
(For assumed-shape arrays and for pointers, Fortran 2008 offers the
CONTIGUOUS attribute.)
(c) Well, that's obvious: If the compiler does not know whether two
variables point to the same memory, a lot of optimizations are prevented
and unneeded temporary variables might be introduced.
Tobias
fj
> On 03/08/2011 07:16 PM, fj wrote:
>> My post was pointing about performance issues between fixed size arrays
>> and allocatables (on needs to retain, in principle, an unknown number of
>> partial derivatives).
>> The module is presently in service in an actual scientific code but,
>> unfortunately, only the fixed size version has finally been selected
>> (much faster)...
>
> In principle there is no reason why a fixed size version should be
> (much) faster than an allocatable-array version.
I agree with you but I observe a clear difference
>
> It is not completely clear how the terms translate into standard speech
> (Fixed size = explicit size with compile-time known size? Allocatable
> array = the dummy argument is allocatable?)
A fixed size vardep is defined as follows
INTEGER,PARAMETER :: ndmax=39
TYPE :: vardep
REAL(c_double) :: value
REAL(c_double) :: deriv(ndmax)
INTEGER(c_int) :: index(ndmax)
INTEGER(c_int) :: nderiv=-1
END TYPE
The version with allocatable is :
TYPE :: vardep
REAL(c_double) :: value
REAL(c_double),allocatable :: deriv(:)
INTEGER(c_int),allocatable :: index(:)
INTEGER(c_int) :: nderiv=-1
END TYPE
>
> I actually see only two (three) potential performance issues (speaking
> of a common implementation, not of all possible implementations the
> standard permits):
>
> a) Passing an assumed-shape array to an explicit-size dummy
No
>
> b) Using an assumed-shaped array in loops or on whole-array operations
They are operations on whole arrays but few because constructing derivatives
usually needs explicit loops.
>
> c) TARGET and POINTER ...
No
>
> For (a), the compiler needs to do a copy in/copy out unless it know that
> the actual argument is contiguous, which it usually does not for
> assumed-shape arrays.
No
>
> For (b) the issue is that the compiler needs to assume that the array is
> not contiguous, thus it needs to take the strides into account.
All the arrays are contiguous here
>
> Depending on the code, either (a) or (b) causes the greater slow down.
> However, if the argument is allocatable, it is contiguous. Thus, both
> issues shouldn't occur and there shouldn't be any larger performance
> difference for your program.
>
> (For assumed-shape arrays and for pointers, Fortran 2008 offers the
> CONTIGUOUS attribute.)
>
> (c) Well, that's obvious: If the compiler does not know whether two
> variables point to the same memory, a lot of optimizations are prevented
> and unneeded temporary variables might be introduced.
>
> Tobias
The problem comes partly from the way physicists are using the module : they
prefer taking advantage of operators + - * ** / ... instead of using
subroutines calls to update existing "vardep" variables (a more efficient
way from the CPU point of view but much less elegant from the programming
one). So many temporaries are created by the compiler with many
allocate/deallocate statements when derivatives and associated indexes are
stored in allocatable arrays.
Replacing allocatable vectors by vectors dimensioned to 39 (hopefully a
relatively low value because of the problems solved) allows to divide the
CPU by a factor between 2 and 5. Why ? I don't know exactly because
temporaries are also involved in that case. This factor has been measured
with several compilers (GCC, INTEL, NAG, SUN, g95) with about the same
result. I already discussed with you about these performance trouble...
without finding any improvement.
--
François Jacq
glen herrmannsfeldt
(snip, someone wrote)
>> In principle there is no reason why a fixed size version should be
>> (much) faster than an allocatable-array version.
> I agree with you but I observe a clear difference
(snip)
> The problem comes partly from the way physicists are using
> the module : they prefer taking advantage of operators + - *
> ** / ... instead of using subroutines calls to update existing
> "vardep" variables (a more efficient way from the CPU point of
> view but much less elegant from the programming one).
> So many temporaries are created by the compiler with many
> allocate/deallocate statements when derivatives and associated
> indexes are stored in allocatable arrays.
Yes, you want to keep the allocate/deallocate out of the nested
loops as much as possible. Usually the array size should depend
on the number of variables, and usually that doesn't change inside
the loop.
The reference to an allocatable array might be slightly slower,
but should be only slightly. It is the repeated allocation
that you want to avoid.
-- glen
Arjen Markus
>
> - Tekst uit oorspronkelijk bericht weergeven -
Interesting, I do not remember seeing it, though.
Is it the one at: http://parfum-echecs.chez-alice.fr/download/var_dep.f90
Have you considered using the parametrised type feature of Fortran
2003?
(I know it is not widely implemented yet, but it could be a way out of
the fixed arrays you now use)
Regards,
Arjen
Aris
Hello Arjen,
I fail to see how your program could be useful. If I understand well,
the user needs to provide an expression. This expression has to come
from somewhere, and I cannot imagine a situation in practice in which
the user herself cannot also provide the derivative of the expression.
The expression may have been generated by a computer algebra system,
and these usually can also provide derivatives. Furthermore, one would
like to be able to manipulate expressions for derivatives. Think of
simple examples like f(x)=log(x)/(1-x) which has "naive" derivatives
that become more and more numerically unstable. Does your program take
care of this?
I don't want to be a b*tch, I am just interested ;-)
Arjen Markus
> Arjen Markus <arjen.markus...@gmail.com> wrote:
> > On 3 mrt, 22:26, Allamarein <matteo.diplom...@gmail.com> wrote:
> >> Let's say I have this function:
>
> >> REAL*8 FUNCTION F(x,y,z)
> >> IMPLICIT NONE
> >> REAL*8, INTENT(IN) :: x,y,z
> >> f = log(x)+y**2+z
> >> END FUNCTION
>
> >> I would compute automatically the three partial derivatives of the
> >> above function by means of a subroutine.
> >> How can I do?
>
> > It currently only works for functions of one variable only, but that
> > is a
> > (mathematical and programmatical) detail.
>
> Hello Arjen,
>
> I fail to see how your program could be useful. If I understand well,
> the user needs to provide an expression. This expression has to come
> from somewhere, and I cannot imagine a situation in practice in which
> the user herself cannot also provide the derivative of the expression.
> The expression may have been generated by a computer algebra system,
> and these usually can also provide derivatives. Furthermore, one would
> like to be able to manipulate expressions for derivatives. Think of
> simple examples like f(x)=log(x)/(1-x) which has "naive" derivatives
> that become more and more numerically unstable. Does your program take
> care of this?
>
> - Tekst uit oorspronkelijk bericht weergeven -
What happens in my module (similar to Francois', btw, only simpler)
is
that the function f(x) = log(x)/(1-x) is evaluated as follows:
x is replaced by a vector (x0,Dx) (because x is the independent
variable).
log(x) is replaced by a vector (log(x0), Dx/x0)
1-x is replaced by a vector (1-x0,-Dx)
etc.
That is:
- For each function the value and the derivative are computed
- For each operation the value and the derivative are computed
- For each part of the complete expression both the value and the
derivative
are propagated
Because of the overloading of the functions, the compiler does the
hardest
part - making a parse tree and evaluating the entire expression.
So, rather than determining the derivative yourself and putting it
into
code (or have it explicitly generated), it is all done implicitly.
It is actually very much like introducing quaternions as an
alternative
to real and complex numbers: define the basic operations and functions
and let the compiler figure out how to combine everything.
Regards,
Arjen
fj
Yes but I can provide a better version if somebody is interested in.
>
> Have you considered using the parametrised type feature of Fortran
> 2003?
> (I know it is not widely implemented yet, but it could be a way out of
> the fixed arrays you now use)
No : our projects need to work with at least four FORTRAN compilers (Intel,
GCC, NAG and SUN or g95)
I don't see the interest of parametrized types because they depend of
parameters which are compilation dependent. I don't understand the
limitation of using parameters. Why not variables ?
For instance, the language ESOPE (a pre-compiler above FORTRAN-77) was able
to generate dynamic segments (kind of derived types depending on integer
variables) at the beginning of eighties ! In additions, arrays inside a
segment were extensible. I speak in the past even if this language is still
in use in many codes developed by CEA (French nuclear research agency), in
particular in mechanics (CAST3M) and neutronics (EURANOS).
>
> Regards,
>
> Arjen
--
François Jacq
fj
Sorry ! Read ERANOS insteed of EURANOS
http://www.oecd-nea.org/tools/abstract/detail/nea-1683
Richard Maine
> No : our projects need to work with at least four FORTRAN compilers (Intel,
> GCC, NAG and SUN or g95)
That certainly rules out parameterized derived types.
> I don't see the interest of parametrized types because they depend of
> parameters which are compilation dependent. I don't understand the
> limitation of using parameters. Why not variables ?
The term is probably confusing you. Type parameters, including those in
parameterized derived types have very little to do with those Fortran
parameter things that I wish had more accurately been called constants,
other than the simillar spelling of the terms. (The standard even does
call parameters named constants. I'm not sure that it ever uses
"parameter" to refer to them, but that's still the keyword and
attribute).
Type parameters do not have to be declared using parameters. Type
parameters come in two flavors. Kind type parameters need to be declared
with initialization expressions, which is to say that they must be known
at compile time. Length type parameters are those that do not have to be
known at compile time; they can be declared with variables. The only
intrinsic length type parameter is the one for character, but the f2003
parameterized derived type feature allows user-defined ones.
However, I'm not at all convinced that PDTs would do anything to help
your problem anyway, even if they were available on the compilers you
use. It seems that the mentioned performance problem relates to
repeatedly allocating and deallocating storage. Just because the term
"parameterized derived type" doesn't include the word "allocatable",
that doesn't magically make the necessity to allocate things go away.
You'll just end up having to declare your variables as allocatable and
pretty much the same allocations will happen.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
Aris
What I wanted to say is that I would not blindly trust the result
of a program that calculates a derivative algebraicly. I would
always want to check the expression for the derivative it comes
up with. And then, I can just as well use more powerfull algebraic
systems than your program to manipulate this expression. Below is
an example of use of your program where you can see the problem.
I just edited your test_diff. In the contains, I write
type(AUTODERIV) function f( x )
type(AUTODERIV), intent(in) :: x
f = log(x)/(1.0-x)
end function f
and in the body
x = derivvar(0.99999)
write(*,*) 'x = ', x
write(*,*) 'log(x)/(1-x) = ', f(x)
with the result
x = 0.99999 1.
log(x)/(1-x) = -1.000005 0.5011338
while the correct result is
correct: -1.000005 0.5000067
I don't think this problem can be solved differently than by analysing
the expression for the derivative, and manipulating it by hand (this
case probably needs a Taylor expansion around x=1).
Arjen Markus
> However, I'm not at all convinced that PDTs would do anything to help
> your problem anyway, even if they were available on the compilers you
> use. It seems that the mentioned performance problem relates to
> repeatedly allocating and deallocating storage. Just because the term
> "parameterized derived type" doesn't include the word "allocatable",
> that doesn't magically make the necessity to allocate things go away.
> You'll just end up having to declare your variables as allocatable and
> pretty much the same allocations will happen.
>
Interesting, I had not thought of that.
An alternative might be to use pointers, rather than allocatables
and then finely control when the allocations/deallocations ought to
happen. F2003's move_alloc statement might help too. (It would
complicate the code, there is no doubt about that!)
Regards,
Arjen
Arjen Markus
>
> What I wanted to say is that I would not blindly trust the result
> of a program that calculates a derivative algebraicly. I would
> always want to check the expression for the derivative it comes
> up with. And then, I can just as well use more powerfull algebraic
> systems than your program to manipulate this expression. Below is
> an example of use of your program where you can see the problem.
> I just edited your test_diff. In the contains, I write
>
> type(AUTODERIV) function f( x )
> type(AUTODERIV), intent(in) :: x
> f = log(x)/(1.0-x)
> end function f
>
> and in the body
>
> x = derivvar(0.99999)
> write(*,*) 'x = ', x
> write(*,*) 'log(x)/(1-x) = ', f(x)
>
> with the result
>
> x = 0.99999 1.
> log(x)/(1-x) = -1.000005 0.5011338
>
> while the correct result is
>
> correct: -1.000005 0.5000067
>
> I don't think this problem can be solved differently than by analysing
> the expression for the derivative, and manipulating it by hand (this
>
> - Tekst uit oorspronkelijk bericht weergeven -
Interesting. This might be a matter of catestrophic cancellations.
Even with double precision you are not entirely out of the woods.
And of course for x = 1, it will fail miserably with a division of
0 by 0.
Another example where this things go awry, but that is with the
function value itself is:
f(x) = (1.0-cos(x))/x
with x near 0
If you reformulate it as:
f(x) = 0.5 * sin(x)**2 / x
you get a much more accurate answer.
I will have to document these issues.
Regards,
Arjen
Tobias Burnus
> On 9 mrt, 20:18, nos...@see.signature (Richard Maine) wrote:
>> However, I'm not at all convinced that PDTs would do anything to help
>> your problem anyway, even if they were available on the compilers you
>> use. It seems that the mentioned performance problem relates to
>> repeatedly allocating and deallocating storage. Just because the term
>> "parameterized derived type" doesn't include the word "allocatable",
>> that doesn't magically make the necessity to allocate things go away.
>> You'll just end up having to declare your variables as allocatable and
>> pretty much the same allocations will happen.
That's actually a quite interesting issue. For "normal" dummy
arguments/variables, one has a nice machinery to avoid performance
issues: explicit shape arrays, contiguous, etc.
However, some solutions do not help any more if one deals with derived
types. With (re)allocate on assignment and other items, one might have
many implicit allocatable checks, allocations and deallocations.
> An alternative might be to use pointers, rather than allocatables
> and then finely control when the allocations/deallocations ought to
> happen.
I am not sure that this helps. You avoid the implicit allocation checks
and (re,de)allocation, but with a high cost: POINTER (and to a lesser
extend: TARGET) mean that different variables can point to the same
memory, which prevents many optimizations. For instance:
real :: a(:), b(:)
! pointer :: a, b
a = b
If "a" and "b" are pointer - or one is pointer and the other target, the
compiler has to create a temporary variable, assign "b" to it and then
assign the temporary to "a" and finally free the temporary. If neither
is a pointer (or one is a pointer and the other is not a target), no
temporary is needed.
In that sense, I wonder whether Fortran needs an equivalent to C's
"restrict" to mark pointers as restricted, combining the properties of
Fortran's allocatables (not alising, unless target) with Fortran's
pointers (no automatic memory handling).
I think for some special cases (like the one in this thread), it seems
to be useful addition. However, in the general case, I think would
always recommend to use "allocatable" where possible and use "pointer"
only when needed (e.g. for linked lists).
Still, one could consider it as proposal for the next Fortran standard
(2008+5 = 2013?).
Tobias
PS: On my wishlist for Fortran 2013(?) is the possibility to disable
host association and only import selected symbols. And I think ANDTHEN
and ORELSE would be useful (Cf. J3/05-134) but seemingly there are some
issues to implement it properly.
PPS: The next Fortran addition will be a Technical Report (TR) on
further interoperability with C. See http://www.nag.co.uk/SC22WG5/ for
the details. If you have time, you should read the document and send
comments to J3 and/or WG5.
James Van Buskirk
news:1jxuzwn.znugnt3ez7neN%nos...@see.signature...
> However, I'm not at all convinced that PDTs would do anything to help
> your problem anyway, even if they were available on the compilers you
> use. It seems that the mentioned performance problem relates to
> repeatedly allocating and deallocating storage. Just because the term
> "parameterized derived type" doesn't include the word "allocatable",
> that doesn't magically make the necessity to allocate things go away.
> You'll just end up having to declare your variables as allocatable and
> pretty much the same allocations will happen.
In a parameterized derived type it is at least in principle possible
to keep the variable sized array next to the other stuff in the
derived type in memory. You can achieve this kind of layout without
parameterized derived types, of course; I was just programming that
kind of stuff today, but it's harder than simply having a
parameterized derived type take care of the layout of data in memory
for you.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
Richard Maine
> On 9 mrt, 20:18, nos...@see.signature (Richard Maine) wrote:
>
> > However, I'm not at all convinced that PDTs would do anything to help
> > your problem anyway, even if they were available on the compilers you
> > use. It seems that the mentioned performance problem relates to
> > repeatedly allocating and deallocating storage. Just because the term
> > "parameterized derived type" doesn't include the word "allocatable",
> > that doesn't magically make the necessity to allocate things go away.
> > You'll just end up having to declare your variables as allocatable and
> > pretty much the same allocations will happen.
>
> An alternative might be to use pointers, rather than allocatables
> and then finely control when the allocations/deallocations ought to
> happen.
I would strongly recommend against that. You'll end up with code that
doesn't work instead of code that works slowly. I note fj's earlier
comment that
>>> The problem comes partly from the way physicists are using the
>>> module : they prefer taking advantage of operators + - * ** / ...
>>> instead of using subroutines calls to update existing "vardep"
>>> variables (a more efficient way from the CPU point of view but much
>>> less elegant from the programming one).
That doesn't sound like the kind of environment where people are
carefully thinking about the implications of using pointers at every
step along the way. Sounds more like one where they want to be insulated
from the issue. Making it so that the pointer thing works right in all
cases and doesn't get accidentally messed up by undisciplined use can be
a tremendous amount of work. You better start with throwing out all
generics that aren't type-bound. Probably better also make components
private.
Remember how long it took before allocatable components took to get to
where they worked reasonably reliably? They were supposed to go in f95,
but got dropped from it because they weren't ready. So that was early
90's when work on them started. A TR came out, but ended up needing a
revision because it still had bugs. Even the revised TR omitted what
turned out to be important parts that were added in f2003. That's all
just in terms of the standard specifying the darn things. Then there was
the "minor detail" of getting at least the majority of the bugs out of
compiler imnplementations. You'd end up retracing much of that history
in essentially redoing the work, presumably with less resources devoted
to it than the Fortran standard had. See you in a decade or two.
Nah. Much better to just not go there at all.
fj
> fj <franco...@irsn.fr> wrote:
>
>> No : our projects need to work with at least four FORTRAN compilers
>> (Intel, GCC, NAG and SUN or g95)
>
> That certainly rules out parameterized derived types.
>
>> I don't see the interest of parametrized types because they depend of
>> parameters which are compilation dependent. I don't understand the
>> limitation of using parameters. Why not variables ?
>
> The term is probably confusing you. Type parameters, including those in
> parameterized derived types have very little to do with those Fortran
> parameter things that I wish had more accurately been called constants,
> other than the simillar spelling of the terms. (The standard even does
> call parameters named constants. I'm not sure that it ever uses
> "parameter" to refer to them, but that's still the keyword and
> attribute).
You are right ! As usual I missed an important point : "parameter" does not
mean "constant" in that context.
>
> Type parameters do not have to be declared using parameters. Type
> parameters come in two flavors. Kind type parameters need to be declared
> with initialization expressions, which is to say that they must be known
> at compile time. Length type parameters are those that do not have to be
> known at compile time; they can be declared with variables. The only
> intrinsic length type parameter is the one for character, but the f2003
> parameterized derived type feature allows user-defined ones.
>
> However, I'm not at all convinced that PDTs would do anything to help
> your problem anyway, even if they were available on the compilers you
> use. It seems that the mentioned performance problem relates to
> repeatedly allocating and deallocating storage. Just because the term
> "parameterized derived type" doesn't include the word "allocatable",
> that doesn't magically make the necessity to allocate things go away.
> You'll just end up having to declare your variables as allocatable and
> pretty much the same allocations will happen.
I am not sure that you are right here. For instance, look at the following
piece of programming :
MODULE var_dep
TYPE vardep(n)
DOUBLE PRECISION :: value
DOUBLE PRECISION :: deriv(n)
INTEGER :: index(n)
INTEGER :: nderiv=-1
INTEGER,LEN :: n=10
END TYPE
...
END
subroutine heat_coefficient(Hconv)
use var_dep
type(vardep),intent(inout) :: Hconv(:)
type(vardep(Hconv(1)%n)) :: rho,v,lambda,mu,T,cp,Re,Pr,Nu,dH
integer :: imesh
do imesh=1,size(Hconv)
call get_fluid_property(imesh,'lambda',lambda)
call get_fluid_property(imesh,'mu',mu)
call get_fluid_property(imesh,'rho',rho)
call get_fluid_property(imesh,'velocity',v)
call get_fluid_property(imesh,'cp',cp)
call get_fluid_property(imesh,'dH',dh)
! Dittus-Boelter
Re=rho*v*dh/mu
Pr=mug*cp/lambda
Nu=0.0243d0*Re**0.8d0*Pr**0.3d0
Hconv(imesh)=lambda*Nu/dH
enddo
end subroutine
In the procedure heat_coefficient, the PDT's are allocated only at the
entrance and not in the loop. With allocatable vardep, allocations normally
occur only when needed, often once for imesh=1, so globally the same work.
In fact, the problem of the allocatable version comes from expressions like
Nu=0.0243d0*Re**0.8d0*Pr**0.3d0. Indeed 4 operators are called, working with
temporary vardeps, which implies at least four masked allocate/deallocate
statements at each loop index for this single instruction.
With PDTs, like with my fixed size version,the temporaries may be allocated
only once, outside the loop (if the compiler is clever enough).
Compared to the fixed version, the dimension of arrays,( Hconv(1)%n)) may be
adjusted to an optimal value less than the overestimation 39. Indeed, I
noticed for instance that the dimension 19 was about 30% less CPU consuming
than 39, probably because of a better localization of data improving the
cache memory management.
So the proposal of Arjen of using PDT to optimize my module, deserves to be
tested. Unfortunately, I cannot do that because I don't have any fortran
compiler accepting PDTs.
--
François Jacq
Richard Maine
> Richard Maine wrote:
> > However, I'm not at all convinced that PDTs would do anything to help
> > your problem anyway, even if they were available on the compilers you
> > use. It seems that the mentioned performance problem relates to
> > repeatedly allocating and deallocating storage. Just because the term
> > "parameterized derived type" doesn't include the word "allocatable",
> > that doesn't magically make the necessity to allocate things go away.
> > You'll just end up having to declare your variables as allocatable and
> > pretty much the same allocations will happen.
>
> I am not sure that you are right here.
...
> With PDTs, like with my fixed size version,the temporaries may be allocated
> only once, outside the loop (if the compiler is clever enough).
I don't see why do PDTs make it inherently more likely for the compiler
to be able to hoist the allocation outside of the loop. Certainly, I
agree that in principle it is possible, but then it seems like it ought
to be possible with allocatable components as well.
Insomuch as vendors mostly don't yet have PDTs implemented at all, I
wouldn't hold my breath while making assumptions about what a great job
of optimization they are likely to do. Even if it turns out that it is
easier to do the optimization with PDTs, compilers have about a decade
of head start in working on allocatable components.
> So the proposal of Arjen of using PDT to optimize my module, deserves to be
> tested. Unfortunately, I cannot do that because I don't have any fortran
> compiler accepting PDTs.
Sure it seems worth testing. I could easily be wrong. The above is all
just my guesses about optimizability issues. I won't try to claim to be
an expert on the matter.
fj
> fj <franco...@irsn.fr> wrote:
>
>> Richard Maine wrote:
>
>> > However, I'm not at all convinced that PDTs would do anything to help
>> > your problem anyway, even if they were available on the compilers you
>> > use. It seems that the mentioned performance problem relates to
>> > repeatedly allocating and deallocating storage. Just because the term
>> > "parameterized derived type" doesn't include the word "allocatable",
>> > that doesn't magically make the necessity to allocate things go away.
>> > You'll just end up having to declare your variables as allocatable and
>> > pretty much the same allocations will happen.
>>
>> I am not sure that you are right here.
> ...
>> With PDTs, like with my fixed size version,the temporaries may be
>> allocated only once, outside the loop (if the compiler is clever enough).
>
> I don't see why do PDTs make it inherently more likely for the compiler
> to be able to hoist the allocation outside of the loop. Certainly, I
> agree that in principle it is possible, but then it seems like it ought
> to be possible with allocatable components as well.
For the temporary variable itself I agree, but * NOT * for a member of that
temporary variable which has the allocatable attribute : even a very clever
compiler cannot take the risk not to deallocate such a member before calling
an operator returning this temporary variable. For instance the use of the
functions SIZE and ALLOCATED inside such operator would make the result
highly unpredictable... even if ironically, with my version of allocatable
vardep, such oversight would be handled successfully !
>
> Insomuch as vendors mostly don't yet have PDTs implemented at all, I
> wouldn't hold my breath while making assumptions about what a great job
> of optimization they are likely to do. Even if it turns out that it is
> easier to do the optimization with PDTs, compilers have about a decade
> of head start in working on allocatable components.
>
>> So the proposal of Arjen of using PDT to optimize my module, deserves to
>> be tested. Unfortunately, I cannot do that because I don't have any
>> fortran compiler accepting PDTs.
>
> Sure it seems worth testing. I could easily be wrong. The above is all
> just my guesses about optimizability issues. I won't try to claim to be
> an expert on the matter.
>
--
François Jacq