Fortran templates

610 views
Skip to first unread message

relaxmike

unread,
Apr 15, 2008, 10:16:41 AM4/15/08
to
Hi all fortran gurus !

I would like to discuss the current methods to manage
templates in fortran, following the thread

http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/f055f870774d98a0/110fa05197d7b709?lnk=gst&q=templates+fortran#110fa05197d7b709

which was very interesting, but showed no details on
how to pratically manage the source code.
This feature is very well-known by C++ developers as "templates".

One example of the problems solved by C++ templates is to
have a sorting source code which is able to manage for any
data type, including integers, reals, or even abstract data
types.

I currently know 3 ways of dealing with templates in fortran,
even if none of them is included in any fortran norm
(and none of them is detailed in a fortran book, to my knowledge) :
- pre-processing macros,
- clever use of the "include" statement,
- m4 macros.

Let's begin with the pre-processing macros.
Suppose that the file "sorting_template.f90" contains a template
sorting module, parametrized by the _QSORT_TYPE macro.
The following is an example of parametrized argument declaration :

subroutine qsort ( array, compare )
_QSORT_TYPE, dimension(:) :: array

Here "_QSORT_TYPE" may be an integer, a real or any fortran
derived type.
This name has been chosen because a fortran variable cannot
begin with an underscore so that no confusion can occur
between a preprocessing macro and a variable name.
Before using the template, one must instanciate it,
which is done with :

#define _QSORT_TYPE type(SORT_DATA)
#include "sorting_template.f90"

The "instanciation" is based, for example, on the following
SORT_DATA derived type, which may be defined in the
module "m_testsorting.f90" :

type SORT_DATA
integer key
integer index
end type SORT_DATA

Even if this derived type is quite simple, practical uses of this
method may include more complex data types.
The module "m_testsorting" has to be pre-processed before
being used. After pre-processing, the previous line has
been transformed to :

subroutine qsort ( array, compare )
type(SORT_DATA), dimension(:) :: array

With this method, generic source code templates can be
configured at will and allow to generate specific
sorting algorithms for whatever type of data.
This method can lead to more complex templates, with
no limit in the number of macros, that is, no limit
in the number of abstract data types in the template.

The main drawback is that the debugging process is not
possible interactively. This is because the source code
is generated at compile-time and the "template" does
not correspond to the post-processed one anymore.
One solution is to use manually the pre-processor to
generate the pre-processed template and to debug on that
later file.

Two other methods exist to manage to do generic programming
in fortran : the "include" statement and m4 macros.

The "include" fortran statement can be used so that the
target source code contains the definitions of an
abstract data type used in the template.
This method is used by Arjen Markus in the Flibs
library, for example in the linked list abstract
data structure :

http://flibs.sourceforge.net/linked_list.html

Here we suppose that the file "sorting_template.f90"
contains one sorting subroutine, which arguments
are defined like this :

subroutine qsort_array( array, compare )
type(SORT_DATA), dimension(:) :: array

The trick is to use the include fortran statement
to put that template source code into a context in
which the abstract data type is defined.
The following source code defines a module, the derived
type SORT_DATA and then include the template source code.

module m_testsorting
type SORT_DATA
integer key
integer index
end type SORT_DATA
contains
include "sorting_template.f90"
end module m_testsorting

The module m_testsorting is now providing the derived type
SORT_DATA and the methods to manage it, which are included
in the template. With a little more work, it is even possible
to make so that the final derived type has a name which
corresponds more to the one implemented by the abstraction,
as shown by Arjen Markus on the linkedlist example.

module MYDATA_MODULE
type MYDATA
character(len=20) :: string
end type MYDATA
end module
module MYDATA_LISTS
use MYDATA_MODULE, LIST_DATA => MYDATA
include "linkedlist.f90"
end module MYDATA_LISTS

The main advantage of the "include" method is that
it uses only fortran statements so that the debugging
process is possible interactively. Moreover, it is very
elegant.

Another method is to use m4 macros to generate source code
at compile-time. Gnu m4 is an implementation of the traditional
Unix macro processor. This method is used by Toby White the
in the Fox library :

http://uszla.me.uk/space/software/FoX/

The template source code is defined in files which have the .m4.
extension and the corresponding source code is defined in the
.F90 file. The source code may be difficult to maintain, but the
method is extremely powerful, thanks to the m4 features.
For example, one can use a m4 "for" loop to generate
several subroutine based on one template subroutine.
But, as for the pre-processing method, the interactive
debugging of the .m4 templates is not possible :
instead, the processed .f90 generated files are debugged easily
and directly.

Are there other well-known methods ?
Is there a definitive drawback for one of these methods so
that another way should be chosen?

All comments will be appreciated.

Best regards,

Michaël Baudin

Arjen Markus

unread,
Apr 15, 2008, 10:30:30 AM4/15/08
to
On 15 apr, 16:16, relaxmike <michael.bau...@gmail.com> wrote:
> Hi all fortran gurus !
>
> I would like to discuss the current methods to manage
> templates in fortran, following the thread
>
> http://groups.google.com/group/comp.lang.fortran/browse_thread/thread...

(One remark: part of the code you refer to was written
by Joe Krahn)

In Fortran 2003 you can use the class facilities and unlimited
polymorphic variables to solve these problems, but material
showing how to do that is limited so far.

Regards,

Arjen

Bil Kleb

unread,
Apr 15, 2008, 12:33:46 PM4/15/08
to
relaxmike wrote:
> Hi all fortran gurus !

Hi.

> Are there other well-known methods ?

Apparently not well-known, but there is:

http://www.macanics.net/forpedo/

Regards,
--
Bil Kleb
http://fun3d.larc.nasa.gov

Terence

unread,
Apr 15, 2008, 7:46:12 PM4/15/08
to
I wrote a commercial general sort/merge program (still in general
public use) that just reads in a text control file with :-
a) the input file name
b) the second input file name (for merges) or blank
c) the output file name
d) record length in bytes
e) then pretty much as many fields as you want, coded as
type,bytestart,bytelength,A/D [,and repeat]

You can also request a tag sort, where the file byte offeset addresses
of the sorted records are left in the output file (for cases where you
don't have twice the original data size available).

Of course it can't deal with arbitary fields beyond the expected
strings, integers, reals and T/F, nor does it sort variable length
records (but could easily changed to do so, nor does it cater for non-
MS byte order and FP formatting), but I see no need for the formal
route suggested.

I always prefer simple and practical.

Damian

unread,
Apr 16, 2008, 1:25:33 AM4/16/08
to
On Apr 15, 7:16 am, relaxmike <michael.bau...@gmail.com> wrote:
> Hi all fortran gurus !
>
> I would like to discuss the current methods to manage
> templates in fortran, following the thread
>
> http://groups.google.com/group/comp.lang.fortran/browse_thread/thread...

At least one text describes an approach equivalent to your pre-
processing approach: Object-Oriented Programming via Fortran 90/95 by
Ed Akin, Cambridge University Press, 2003. I don't have my copy at
hand but I'm pretty sure it's in Chapter 6.

Also, FWIW, some would argue that Fortran 2003 supports at least some
form of generic programming via parametrized derived types (definitely
not the same capability as templates but generic programming
nonetheless) and then there is the class(*) to which Arjen alludes.

Damian

Gerry Ford

unread,
Apr 16, 2008, 2:10:31 AM4/16/08
to

"Arjen Markus" <arjen....@wldelft.nl> wrote in message
news:0a01dacb-fcd7-4279...@b5g2000pri.googlegroups.com...

In Fortran 2003 you can use the class facilities and unlimited
polymorphic variables to solve these problems, but material
showing how to do that is limited so far.

--->Would you consider these class facilities and unlimited polymorphic
variables as "object-oriented?"

--
"Shopping for toilets isn't the most fascinating way to spend a Saturday
afternoon. But it beats watching cable news."

~~ Booman


Arjen Markus

unread,
Apr 16, 2008, 2:50:46 AM4/16/08
to
On 16 apr, 08:10, "Gerry Ford" <ge...@nowhere.ford> wrote:
> "Arjen Markus" <arjen.mar...@wldelft.nl> wrote in message

>
> news:0a01dacb-fcd7-4279...@b5g2000pri.googlegroups.com...
>
> In Fortran 2003 you can use the class facilities and unlimited
> polymorphic variables to solve these problems, but material
> showing how to do that is limited so far.
>
> --->Would you consider these class facilities and unlimited polymorphic
> variables as "object-oriented?"
>

They can certainly be used in that way - for a suitable
definition of "object-oriented". But the main point is
that they give you the possibility to define interfaces
that allow a wide range of actual types.

Regards,

Arjen

Arjen Markus

unread,
Apr 16, 2008, 2:57:17 AM4/16/08
to
> Damian- Tekst uit oorspronkelijk bericht niet weergeven -
>
> - Tekst uit oorspronkelijk bericht weergeven -

I must have a look again at Akin's book - the details
have escaped me ;).

Note that templates in C++ are completely different
beasts than what we are talking about here:

If you define a template A with one parameter T in C++ and
then define in various pieces of your code the same instance,
say a type A< int> in ten different source files,
the compiler will have to create ten different copies of
the actual object code that implements that particular
type.

It is up to the linker to make sure that in the final
executable program there is only one copy (or at least that
there is no conflict between them).

So, in C++ the task is much more complicated than just
generating specific source code from generic source code.


In Fortran we would have to generate the specific source
code for each specific type. And as Michael indicates
there are several methods for that. But in the end,
for the compiler and linker it is just a larger set
of ordinary code!

Regards,

Arjen

Arjen Markus

unread,
Apr 16, 2008, 2:59:23 AM4/16/08
to

I often use Tcl scripts or Fortran programs to
generate the specific code, if other methods
fail - fairly ad hoc, but it works splendidly.

Regards,

Arjen

relaxmike

unread,
Apr 16, 2008, 4:15:40 AM4/16/08
to
On 16 avr, 07:25, Damian <dam...@rouson.net> wrote:
> At least one text describes an approach equivalent to your pre-
> processing approach: Object-Oriented Programming via Fortran 90/95 by
> Ed Akin, Cambridge University Press, 2003. I don't have my copy at
> hand but I'm pretty sure it's in Chapter 6.

I have the book in my hands : "Object oriented programming
via fortran 90/95" by Ed Akin.
In fact, I was very enthousiast about the title, but the content
was very disappointing to me.
For me, the single article "Object-based programming in Fortran 90"
by Gray and Roberts is much more interesting, because it gives
the essence of OO :
http://www.ccs.lanl.gov/CCS/CCS-4/pdf/obf90.pdf

About the problem of generic code, here is the only part where the
subject is detailed, in chapter 3, p. 51 :
"None of the OOP languages have all the features one might desire.
For example, the useful concept of <<template>>, which is standard in C
++,
is not in the F90 standard. Yet the author has found that a few
dozen lines of F90 code will define a preprocessor that allows
templates
to be defined in F90 and expanded in line at compile time. [...]"

That is all : I did not find any paragraph showing how to practically
implement templates in fortran, a subject which could have been at
the core of the book with such a title.
The second edition of "Fortran 90/95 explained" does not deal
with that subject neither (but it may be done in the current edition,
I don't know).

Comparing to the methods I allready know, Forpedo seems very
interesting.
The tool is easy to implement within a Makefile, as for m4 macros.
The generic source code really looks like fortran source code,
even with the macros like @T and <T>.

I did the following table to compare the methods
- standard fortran : is the method in the fortran norm
- debug : is the template source code easy to debug
- features : does the method allow to manage complex templates
(this is only my point of view)

Standard Debug Features
Fortran
-----------------------------------------------------------
include * ** *
Pre-processing - * **
m4 - * ****
forpedo - * ***
fortran 2003 * ? ?

Michaël

Damian

unread,
Apr 16, 2008, 10:23:32 AM4/16/08
to

Keep looking. There is a short passage with more details later in the
book. Again, I'm pretty sure it's chapter 6, but I'll send an exact
page later when I have my copy with me. The reason I found it similar
to what you wrote is your use of the leading underscore, which you
explained was because Fortran variables names cannot start with an
underscore. In his case, he uses a trailing dollar sign ($) for a
similar reason. Look for the string "Template$". The reason I'm so
certain it's there is because I use his technique in my own code. I
learned it from his book and haven't seen it anywhere else.

Damian

Damian

unread,
Apr 16, 2008, 10:25:45 AM4/16/08
to

Excellent points! Thanks for the clarification.

Damian

relaxmike

unread,
Apr 16, 2008, 11:14:52 AM4/16/08
to
You are right, I did not remember that part, but I read
it.
This is in chapter 6, p144, "Templates".
My point of view is different from the one of the
author about that problem and I hope that this chapter
will be rewritten in the second edition, if ever.

Michaël

gunnar

unread,
Apr 16, 2008, 1:36:41 PM4/16/08
to
On 15 Apr, 16:16, relaxmike <michael.bau...@gmail.com> wrote:
> Hi all fortran gurus !
>
>
> I currently know 3 ways of dealing with templates in fortran,
> even if none of them is included in any fortran norm
> (and none of them is detailed in a fortran book, to my knowledge) :
> - pre-processing macros,
> - clever use of the "include" statement,
> - m4 macros.
>


Hi,

From what I have experienced, after starting to program in Fortran
after 15 years, is that there is very boring to write the same control
logic all over the program, especially when you have to make a minor
change and need to do the same update in many places. Not to metion
the probability to introduce errors by doing so. I also find it boring
to write the logic needed to be able to select a certain algorithm
during program execution.

As I see it, parametrized derived types and CLASS(*) has the intention
to let you do this, but I think that you have to include something
more or extend the these constructs to really have something that
helps.

One idea that I have: Would it not be possible to also parametrize the
procedures?

Something like this

subroutine bbbbb(E,F,H,D)
type,type :: A=(INTEGER,REAL),B=(TypeX,typeY),C
character(10), length :: D=("AlgorithmA")
type(A),allocateable :: E(:)
class(B) :: B
type(C) :: H
.
.
.
end


The compiler could then use the type of E and F, and the value of D,
to select the specific subroutine to use, and I could use D to select
the algorithm to use during program execution.


Regards Gunnar


GaryScott

unread,
Apr 16, 2008, 4:39:46 PM4/16/08
to
On Apr 16, 12:36 pm, gunnar <bjorkm...@gmail.com> wrote:
> On 15 Apr, 16:16, relaxmike <michael.bau...@gmail.com> wrote:
>
<snip>

> From what I have experienced, after starting to program in Fortran
> after 15 years, is that there is very boring to write the same control
> logic all over the program, especially when you have to make a minor
> change and need to do the same update in many places. Not to metion
> the probability to introduce errors by doing so. I also find it boring
> to write the logic needed to be able to select a certain algorithm
> during program execution.

While I have experienced this code duplication "problem" on rare
occasion, I find that I fail to see that it is a great programming
concern. I don't oppose adding such capabilities to the language, but
please make it a formal part of the language and not some add-on that
will be poorly supported like a pre-processor that largely, probably,
poorly duplicates functionality of a pre-existing native JCL or macro
language or a more popular preprocessor for another language.

<snip>

Regards Gunnar

relaxmike

unread,
Apr 17, 2008, 7:58:16 AM4/17/08
to
On 16 avr, 22:39, GaryScott <garylsc...@sbcglobal.net> wrote:
> While I have experienced this code duplication "problem" on rare
> occasion, I find that I fail to see that it is a great programming
> concern.

That really depends on what kind of problems you have to solve.
"If you only have a hammer, everything you see is a nail."

Suppose you have to connect a fortran software
with a Java platform and transfer data between the two with
sockets. The subroutines which are used to manage the send/receive
are complex and have to be able to manage all basic
fortran data types : integer, real, character, and all possible
arrays from rank 1 to rank 7. Let us count 3 * 7 = 21 subroutines
to send data and 21 subroutines to receive data. But you also
have to manage errors (21 more), etc...
In short : the connection is not possible to do in fortran without
fortran templates. The conclusion is : "Aaahh, that fortran language
is crap, let's do it in C++ and the job will be done.".
You think that this case cannot happen ?
Yes it can happen, and it really happened 4 years ago in my last job.
So my colleagues used pre-processing, which is a very good idea,
and that worked.

In the case of Fox, that tool would simply not exist if fortran
templates did not exist (based on m4 in that case).

The fact that "fortran is crap" or that "fortran templates are not
a great concern" depends on what solutions you can put in front of
the problems you have.

Regards,

Michaël Baudin

Gary Scott

unread,
Apr 17, 2008, 8:27:21 AM4/17/08
to
relaxmike wrote:
> On 16 avr, 22:39, GaryScott <garylsc...@sbcglobal.net> wrote:
>
>>While I have experienced this code duplication "problem" on rare
>>occasion, I find that I fail to see that it is a great programming
>>concern.
>
>
> That really depends on what kind of problems you have to solve.
> "If you only have a hammer, everything you see is a nail."
>
> Suppose you have to connect a fortran software
> with a Java platform and transfer data between the two with
> sockets. The subroutines which are used to manage the send/receive
> are complex and have to be able to manage all basic
> fortran data types : integer, real, character, and all possible
> arrays from rank 1 to rank 7. Let us count 3 * 7 = 21 subroutines
> to send data and 21 subroutines to receive data. But you also
> have to manage errors (21 more), etc...
> In short : the connection is not possible to do in fortran without
> fortran templates. The conclusion is : "Aaahh, that fortran language
> is crap, let's do it in C++ and the job will be done.".

I've been doing sockets without templates for 15 years. I'm not
familiar enough with your above application comment on it though. I
didn't imply there was no use for templates. I was suggesting that if
we implement them, do it right, not through some half-a$$ preprocessor
method. I oppose preprocessors, because whatever gets put in the
standard will basically be inferior to tools I had on one or another
operating system in the past and I'll always know how inferior it is as
I try to force it to do things it wasn't designed for.

<snip>
>
> Regards,
>
> Michaël Baudin


--

Gary Scott
mailto:garylscott@sbcglobal dot net

Fortran Library: http://www.fortranlib.com

Support the Original G95 Project: http://www.g95.org
-OR-
Support the GNU GFortran Project: http://gcc.gnu.org/fortran/index.html

If you want to do the impossible, don't hire an expert because he knows
it can't be done.

-- Henry Ford

relaxmike

unread,
Apr 17, 2008, 11:18:58 AM4/17/08
to
On 17 avr, 14:27, Gary Scott <garylsc...@sbcglobal.net> wrote:
> I oppose preprocessors, because whatever gets put in the
> standard will basically be inferior to tools I had on one or another
> operating system in the past and I'll always know how inferior it is as
> I try to force it to do things it wasn't designed for.

You made a point here, because the pre-processing method is
difficult to debug directly, except if you debug the pre-processed
source code instead of the template one.
This can be easily solved though, via makefile commands such as :

%.f90 : %.F90
$(FPP) $(FPPOPTS) $< -P $@

%.o : %.f90
$(FC) $(F90OPTS) -c $< -o $@

That is to say that the .F90 contains the template to pre-process.
make start to compute the .f90 specific source code, depending
on the template one.
It is only the .f90 pre-processed source code which has to be
compiled. This allows to debug interactively the specific
source code, while maintaining only the template one.
The "-P" option of the fortran preprocessor may allow to remove
the lines introduced by the preprocessor in the preprocessed
file which points to the line number in the original template.

Now that I have revealed all my very secret tricks (!!!), I don't
see any main drawback to using pre-processing as a template method,
except that the features implemented in many pre-processor are very
poor, compared to what is possible in Python, m4 or Tcl for example.

Most compilers include a pre-processor inside the
compiler itself, for example Intel Fortran, gfortran, g95, Sun,
Ibm, etc...
Even if the pre-processor is not included in
your favorite compiler, or if you dislike the features
that it provides, you can use an external one, for example
the C pre-processor. If the C pre-processor is not available
on the machine (which has a very low probability, but is possible),
you can get the source of gcc and compile your own release.

In fact, the main reason why pre-processing macros are not used
widely in fortran is for historical reasons, I suppose, and by the
fact
that many fortran developers are developing only in fortran,
with high-level skills in science, most of the time, and, most
of the time, low-level skills in software : there are many,
many exceptions to this general fact, and these exceptions are
the core of this forum.
Put a C developer in front of a fortran source, and he will
introduce a new pre-processing macro every 3 minutes...

What I am curious is at the fortran 2003 features which
allows to do this. I think I'm going to read the standard
again.

Best regards,

Michaël

Richard Maine

unread,
Apr 17, 2008, 12:29:54 PM4/17/08
to
relaxmike <michael...@gmail.com> wrote:

> Suppose you have to connect a fortran software
> with a Java platform and transfer data between the two with
> sockets. The subroutines which are used to manage the send/receive
> are complex and have to be able to manage all basic
> fortran data types : integer, real, character, and all possible
> arrays from rank 1 to rank 7. Let us count 3 * 7 = 21 subroutines
> to send data and 21 subroutines to receive data. But you also
> have to manage errors (21 more), etc...
> In short : the connection is not possible to do in fortran without
> fortran templates.

I agree that something like templates can be useful in some situations.
But note that, depending on details, the particular scenario described
above can usually be handled better in other ways. Typically, all you
need is an equivalent of a C void pointer, along with the data size.
Then you can do it all in a single procedure, with no need for
templates. I've done this exact thing for a long time, as I also have
Fortran code that passes data via sockets. So have other people. This
one is done a lot.

You end up with a single procedure, rather than the multiplicity of
procedures that the template approach involves. The template approach
makes it more practical to create that multiplicity of procedures, but
they still do end up having to get created.

In Fortran 90/95, you have to use tricks that aren't technically
standard conforming, but still have decent portability in practice, in
order to do the needed "type cheating". In Fortran 2003, you can use
either the C interop features or the unlimitted polymorphic feature to
make such code standard conforming.

--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain

relaxmike

unread,
Apr 18, 2008, 6:07:30 AM4/18/08
to
On 17 avr, 18:29, nos...@see.signature (Richard Maine) wrote:
> I agree that something like templates can be useful in some situations.
> But note that, depending on details, the particular scenario described
> above can usually be handled better in other ways. Typically, all you
> need is an equivalent of a C void pointer, along with the data size.
> Then you can do it all in a single procedure, with no need for
> templates. I've done this exact thing for a long time, as I also have
> Fortran code that passes data via sockets. So have other people. This
> one is done a lot.
>
> You end up with a single procedure, rather than the multiplicity of
> procedures that the template approach involves. The template approach
> makes it more practical to create that multiplicity of procedures, but
> they still do end up having to get created.

I don't understand how to implement this way. What kind of fortran
data type corresponds to a C void pointer ?
Let's take a practical example or reading a data in a string.
The following source code read a data in an internal file
and returns the value read, depending on the type of variable.
It is over-simplified to make the point shorter.
It does not use fortran templates, but only standard fortran.

interface readfile_data
module procedure readfile_data_logical
module procedure readfile_data_integer
end interface readfile_data

subroutine readfile_data_logical ( current_data_line , myvalue ,
status )
implicit none
character ( len = 200 ) , intent(in) :: current_data_line
logical , intent (out) :: myvalue
integer, intent(out) :: status
status = 0
read ( current_data_line , * , err = 100 , end = 100) myvalue
return
100 continue
status = 1
end subroutine readfile_data_logical

subroutine readfile_data_integer ( current_data_line , myvalue ,
status )
implicit none
character ( len = 200 ) , intent(in) :: current_data_line
integer , intent (out) :: myvalue
integer, intent(out) :: status
status = 0
read ( current_data_line , * , err = 100 , end = 100) myvalue
return
100 continue
status = 1
end subroutine readfile_data_integer

It is very important to factorise this small code, because the
interface
could be extended way beyond the current state, for example reading
more
complex data types :

interface readfile_data
module procedure readfile_data_logical
module procedure readfile_data_integer
module procedure readfile_data_real
module procedure readfile_data_double
module procedure readfile_data_logical_array1
module procedure readfile_data_integer_array1
module procedure readfile_data_real_array1
module procedure readfile_data_double_array1
module procedure readfile_data_logical_array2
module procedure readfile_data_integer_array2
module procedure readfile_data_real_array2
module procedure readfile_data_double_array2
[etc... until the maximum rank size available in fortran is reached
= 7]
end interface readfile_data

with the following header for the readfile_data_double_array2
subroutine :

subroutine readfile_data_double_array2 ( current_data_line ,
myvalue , status )
implicit none
character ( len = 200 ) , intent(in) :: current_data_line
real(dp), dimension(:,:) , intent (out) :: myvalue
integer, intent(out) :: status

Developing the code with standard fortran is near to impossible,
although hiring a PhD with European funds to do this is allways
possible...
With fortran templates, it is really easy.
Suppose that the file "readfile_data_template.F90" contains the
following
source code, parametrized with the "_DATATYPE" macro.

subroutine readfile_data__DATATYPE ( current_data_line , myvalue ,
status )
implicit none
character ( len = 200 ) , intent(in) :: current_data_line
_DATATYPE , intent (out) :: myvalue
integer, intent(out) :: status
status = 0
read ( current_data_line , * , err = 100 , end = 100) myvalue
return
100 continue
status = 1
end subroutine readfile_data_integer

The previous source code would be using the template like this :

interface readfile_data
module procedure readfile_data_logical
module procedure readfile_data_integer
end interface readfile_data

#define _DATATYPE logical
#include "readfile_data_template.F90"

#define _DATATYPE integer
#include "readfile_data_template.F90"

And now, with only 6 lines of source code, extending it to
manage for several data types, and arrays of all possible
sizes is possible. The same problem would be even easier to
manage with Forpedo, I suppose.

Now, how do you suggest to solve the problem with your method
Richard ? What is the fortran type which corresponds to the
C void pointer ?

Regards,

Michaël

relaxmike

unread,
Apr 18, 2008, 8:51:03 AM4/18/08
to
I still try to experiment the idea, so here is a sample
full demonstration of the pre-processing way.
Here is the file "test2_template.F90" :

subroutine _READFILE_DATA_NAME ( current_data_line , myvalue ,
status )
implicit none
character ( len = * ) , intent(in) :: current_data_line


_DATATYPE , intent (out) :: myvalue
integer, intent(out) :: status
status = 0
read ( current_data_line , * , err = 100 , end = 100) myvalue
return
100 continue
status = 1

end subroutine _READFILE_DATA_NAME

And this is the test file :

module m_moduletest2


interface readfile_data
module procedure readfile_data_logical
module procedure readfile_data_integer
end interface readfile_data

contains

#define _DATATYPE logical
#define _READFILE_DATA_NAME readfile_data_logical
#include "test2_template.F90"
#undef _DATATYPE
#undef _READFILE_DATA_NAME

#define _DATATYPE integer
#define _READFILE_DATA_NAME readfile_data_integer
#include "test2_template.F90"
#undef _DATATYPE
#undef _READFILE_DATA_NAME

end module m_moduletest2

program test2
use m_moduletest2
character (len=200) :: current_data_line
integer :: myintvalue1
logical :: mylogicalvalue1
integer :: status
current_data_line = "1"
call readfile_data ( current_data_line , myintvalue1 , status )
write ( * , * ) "status:", status
write ( * , * ) "Integer : ", myintvalue1
current_data_line = ".true."
call readfile_data ( current_data_line , mylogicalvalue1 , status )
write ( * , * ) "status:", status
write ( * , * ) "Logical : ", mylogicalvalue1
end program test2

I tried to use "include" only statements without pre-processing,
and that lead to a source code which I am not proud of.
But the "C void pointer" idea leads to the following source.
The "void" idea is managed with a derived type containing all
possible fortran basic data types.
The "pointer" is the class name, a string containing "integer",
or "logical", depending on the type to manage.
This is the source code :

module m_moduletest4
type :: DATATYPE
integer :: value_integer
logical :: value_logical
character ( len = 200 ) :: classname
end type DATATYPE
contains
subroutine readfile_data ( current_data_line , classname , myvalue ,
status )
implicit none
character ( len = * ) , intent(in) :: current_data_line
character ( len = * ) , intent(in) :: classname
type ( DATATYPE ) , intent(out) :: myvalue


integer, intent(out) :: status
status = 0

myvalue % classname = classname
select case ( classname )
case ( "integer" )
read ( current_data_line , * , err = 100 , end = 100) myvalue %
value_integer
case ( "logical" )
read ( current_data_line , * , err = 100 , end = 100) myvalue %
value_logical
case default
write(6,*) "Bad classname."
end select


return
100 continue
status = 1
end subroutine readfile_data

subroutine printdata ( myvalue )
implicit none
type ( DATATYPE ) , intent(in) :: myvalue
write ( * , * ) trim(myvalue % classname) , ":"
select case ( myvalue % classname )
case ( "integer" )
write ( * , * ) myvalue % value_integer
case ( "logical" )
write ( * , * ) myvalue % value_logical
case default
write(6,*) "Bad classname."
end select
end subroutine printdata
end module m_moduletest4

program test4
use m_moduletest4
character (len=200) :: current_data_line
type ( DATATYPE ) :: myvalue
integer :: status
current_data_line = "1"
call readfile_data ( current_data_line , "integer", myvalue ,
status )
write ( * , * ) "status:", status
call printdata ( myvalue )
current_data_line = ".true."
call readfile_data ( current_data_line , "logical" , myvalue ,
status )
write ( * , * ) "status:", status
call printdata ( myvalue )
end program test4


Of course, dealing with that implementation consumes a lot of memory,
especially if we include arrays of integers, arrays of logicals,
etc...
To solve that, we can use fortran 90 pointers and allocate only the
type under use. Since the abstract data type is more complex, it
is now time for full OO.

module m_moduletest5
type :: DATATYPE
integer, pointer :: value_integer => NULL()
logical, pointer :: value_logical => NULL()
end type DATATYPE
contains
subroutine data_newfromstring ( this , current_data_line ,
classname , status )
implicit none
character ( len = * ) , intent(in) :: current_data_line
character ( len = * ) , intent(in) :: classname
type ( DATATYPE ) , intent(out) :: this


integer, intent(out) :: status
status = 0

select case ( classname )
case ( "integer" )
allocate ( this % value_integer )
read ( current_data_line , * , err = 100 , end = 100) this %
value_integer
case ( "logical" )
allocate ( this % value_logical )
read ( current_data_line , * , err = 100 , end = 100) this %
value_logical
case default
write(6,*) "Bad classname."
end select


return
100 continue
status = 1

end subroutine data_newfromstring
subroutine data_print ( this )
implicit none
type ( DATATYPE ) , intent(in) :: this
if ( associated ( this % value_integer ) ) then
write ( *, * ) "integer :"
write ( * , * ) this % value_integer
elseif ( associated ( this % value_logical ) ) then
write ( *, * ) "logical :"
write ( * , * ) this % value_logical
endif
end subroutine data_print
subroutine data_free ( this )
implicit none
type ( DATATYPE ) , intent(in) :: this
if ( associated ( this % value_integer ) ) then
deallocate ( this % value_integer )
elseif ( associated ( this % value_logical ) ) then
deallocate ( this % value_logical )
endif
end subroutine data_free

end module m_moduletest5

program test5
use m_moduletest5
character (len=200) :: current_data_line
type ( DATATYPE ) :: myvalue
integer :: status
current_data_line = "1"
call data_newfromstring ( myvalue , current_data_line , "integer",
status )
write ( * , * ) "status:", status
call data_print ( myvalue )
call data_free ( myvalue )
current_data_line = ".true."
call data_newfromstring ( myvalue , current_data_line , "logical" ,
status )
write ( * , * ) "status:", status
call data_print ( myvalue )
call data_free ( myvalue )
end program test5

I admit that the current code is still manageable, but have many
limitations
that the pre-processing version have not, including :
- the memory is managed at hand, which can lead to memory leaks.
This is easy to manage with only 2 basic types, but what if there are
21 ?
- the abstract data type cannot handle user-defined derived-types.
The pre-processing system allows to define whatever type you want to,
without any complication in the client source code.
This is not the case with the hand-crafted "pointer to everything"
abstract data type.
On the good side, the source code is very easy to debug and uses only
standard fortran 90 statements.

Still the current "pointer to everything" is a heavy solution.
What if we wand to define a "writetostring" method : another 21*3
block of source code. And what if we want to define an error
system for the class : another 21*3 block of source code !
This is not practical.

But you may suggest another way ?

Best regards,

Michaël

Richard Maine

unread,
Apr 18, 2008, 11:13:22 AM4/18/08
to
relaxmike <michael...@gmail.com> wrote:

> What kind of fortran
> data type corresponds to a C void pointer ?

Type(C_PTR) in the F2003 C interop stuff.

James Van Buskirk

unread,
Apr 18, 2008, 2:31:00 PM4/18/08
to
"relaxmike" <michael...@gmail.com> wrote in message
news:5be1dfe2-56d4-4a32...@k37g2000hsf.googlegroups.com...

> end module m_moduletest2

OK, how about:

C:\gfortran\clf\template_war>type test3_template.i90
private
public readfile_data
interface readfile_data
module procedure READFILE_DATA_NAME
end interface readfile_data
contains
subroutine READFILE_DATA_NAME(current_data_line,Qmyvalue,status)
character(len=*), intent(in) :: current_data_line
intent (out) :: Qmyvalue


integer, intent(out) :: status
status = 0

read(current_data_line, *, err=100 , end=100) Qmyvalue


return
100 continue
status = 1

end subroutine READFILE_DATA_NAME

C:\gfortran\clf\template_war>type test3.f90
module logical_mod
implicit logical (Q)
include 'test3_template.i90'
end module logical_mod

module integer_mod
implicit integer (Q)
include 'test3_template.i90'
end module integer_mod

module m_moduletest3
use logical_mod
use integer_mod
end module m_moduletest3

program test3
use m_moduletest3


character (len=200) :: current_data_line
integer :: myintvalue1
logical :: mylogicalvalue1
integer :: status
current_data_line = "1"
call readfile_data(current_data_line, myintvalue1, status)
write(*, *) "status:", status

write(*, *) "Integer : ", myintvalue1


current_data_line = ".true."
call readfile_data(current_data_line, mylogicalvalue1, status )
write(*, *) "status:", status

write(*, *) "Logical : ", mylogicalvalue1
end program test3

C:\gfortran\clf\template_war>C:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
te
st3.f90 -otest3

C:\gfortran\clf\template_war>test3
status: 0
Integer : 1
status: 0
Logical : T

Isn't this an exercise in shooting fish in a barrel?

C:\gfortran\clf\template_war>type test6.f90
module m_moduletest6
use ISO_C_BINDING, only: C_PTR, C_LOC, C_F_POINTER
implicit none
private C_PTR, C_LOC, C_F_POINTER
type :: DATATYPE
type(C_PTR) value


character ( len = 200 ) :: classname
end type DATATYPE
contains
subroutine readfile_data(current_data_line,classname,myvalue,status)
implicit none
character(len=*) , intent(in) :: current_data_line
character(len=*) , intent(in) :: classname

type(DATATYPE), intent(out) :: myvalue
integer, intent(out) :: status
integer, pointer :: pi4
logical, pointer :: pL4

status = 0
myvalue%classname = classname
select case(classname)
case ("integer")
allocate(pi4)
read(current_data_line, *, err=101, end=101) pi4
myvalue%value = C_LOC(pi4)
return
101 deallocate(pi4)
case("logical")
allocate(pL4)
read(current_data_line, *, err=102, end=102) pL4
myvalue%value = C_LOC(pL4)
return
102 deallocate(pL4)


case default
write(6,*) "Bad classname."
end select
return
100 continue
status = 1
end subroutine readfile_data
subroutine printdata(myvalue)
implicit none
type(DATATYPE), intent(in) :: myvalue

integer, pointer :: pi4
logical, pointer :: pL4

write(* ,*) trim(myvalue%classname), ":"
select case(myvalue%classname)
case("integer")
call C_F_POINTER(myvalue%value, pi4)
write(* ,*) pi4
case("logical")
call C_F_POINTER(myvalue%value, pL4)
write(* ,*) pL4


case default
write(6,*) "Bad classname."
end select
end subroutine printdata

end module m_moduletest6

program test6
use m_moduletest6
character (len=200) :: current_data_line
type(DATATYPE) :: myvalue


integer :: status
current_data_line = "1"
call readfile_data(current_data_line, "integer", myvalue, status)
write(*, *) "status:", status
call printdata(myvalue)
current_data_line = ".true."
call readfile_data(current_data_line, "logical", myvalue, status)
write(* ,*) "status:", status
call printdata(myvalue)

end program test6

C:\gfortran\clf\template_war>C:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
te
st6.f90 -otest6

C:\gfortran\clf\template_war>test6
status: 0
integer:
1
status: 0
logical:
T

> - the abstract data type cannot handle user-defined derived-types.
> The pre-processing system allows to define whatever type you want to,
> without any complication in the client source code.
> This is not the case with the hand-crafted "pointer to everything"
> abstract data type.

Not the case, see my examples above.

> But you may suggest another way ?

Naturally.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


Charles Coldwell

unread,
Apr 20, 2008, 3:10:02 PM4/20/08
to
relaxmike <michael...@gmail.com> writes:

> One example of the problems solved by C++ templates is to
> have a sorting source code which is able to manage for any
> data type, including integers, reals, or even abstract data
> types.
>
> I currently know 3 ways of dealing with templates in fortran,
> even if none of them is included in any fortran norm
> (and none of them is detailed in a fortran book, to my knowledge) :
> - pre-processing macros,
> - clever use of the "include" statement,
> - m4 macros.

You forgot "callbacks":

! Fortran 95 implementation of the quicksort algorithm. This
! subroutine does not directly touch the array it sorts; rather it
! relies on the two callbacks "compare" and "exchange" for that. Code
! inspired by R. Sedgewick, "Algorithms in C" and correspondence with
! Glen Herrmannsfeldt.

subroutine quicksort(n, compare, exchange)
implicit none
integer, intent(in) :: n ! the length of the implied array

! The compare function must return an integer that is
! greater than zero if element(i) > element(j)
! equal to zero if element(i) = element(j)
! less than zero if element(i) < element(j)
interface
integer function compare(i,j)
integer, intent(in) :: i, j
end function compare
end interface

! The exchange subroutine exchanges the elements at locations i and j.
interface
subroutine exchange(i,j)
integer, intent(in) :: i, j
end subroutine exchange
end interface

integer, dimension(2,n) :: stack
integer :: sptr, left, right, pivot

sptr = 1
call push(1, n)

do while(pop(left, right))
if (left .ge. right) cycle
pivot = partition(left, right)
if (pivot .gt. (left+right)/2) then
call push(left,pivot-1)
call push(pivot+1,right)
else
call push(pivot+1,right)
call push(left,pivot-1)
end if
end do
return

contains

subroutine push(l, r)
integer, intent(in) :: l, r

stack(1, sptr) = l
stack(2, sptr) = r
sptr = sptr + 1
end subroutine push

logical function pop(l, r)
integer, intent(out) :: l, r

if (sptr .gt. 1) then
sptr = sptr - 1
l = stack(1, sptr)
r = stack(2, sptr)
pop = .true.
else
pop = .false.
end if
end function pop

integer function partition(l,r)
integer, intent(in) :: l, r
integer :: i, j

i = l
j = r - 1
do
do while (compare(i,r) .lt. 0)
i = i+1
end do
do while (compare(j,r) .gt. 0)
if (j .eq. l) exit
j = j-1
end do
if (i .ge. j) exit
call exchange(i,j)
end do
call exchange(i,r)
partition = i
end function partition
end subroutine quicksort

Chip

--
Charles M. "Chip" Coldwell
"Turn on, log in, tune out"
GPG Key ID: 852E052F
GPG Key Fingerprint: 77E5 2B51 4907 F08A 7E92 DE80 AFA9 9A8F 852E 052F

relaxmike

unread,
Apr 21, 2008, 8:27:30 AM4/21/08
to
Thanks for these suggestions.

I must say that the method based on implicit statements is
very clever, and, based only on fortran statements, allows
to minimize the code duplication. But is it possible to
declare something like this :

module data_mod


type MYDATA
character(len=20) :: string
end type MYDATA

implicit type(MYDATA) (Q)
include 'test3_template.i90'
end module data_mod

I don't think so, which shows that the method cannot be extended
to abstract data types. But funny though.

The method based on "ISO_C_BINDING, only: C_PTR, C_LOC, C_F_POINTER"
is interesting but leads to code duplication.
If I where to use "ISO_C_BINDING", I think that I would use the
C++ templates, and only define in fortran the interface to the
C++ source code.

All in all, it would be much more simpler if the fortran language
include a "template" feature in the core, which may be done,
in 2043, may be...

Michaël

Richard Maine

unread,
Apr 21, 2008, 10:52:18 AM4/21/08
to
relaxmike <michael...@gmail.com> wrote:

> I must say that the method based on implicit statements is
> very clever, and, based only on fortran statements, allows
> to minimize the code duplication. But is it possible to
> declare something like this :
>
> module data_mod
> type MYDATA
> character(len=20) :: string
> end type MYDATA
> implicit type(MYDATA) (Q)
> include 'test3_template.i90'
> end module data_mod
>
> I don't think so, which shows that the method cannot be extended
> to abstract data types.

Yes, it is possible to do something like that. The only problem with the
above is one of ordering. Implicit statements have to come quite early
in a scoping unit, before pretty much anything other than USE
statements. I forget whether it would be ok to just swap the order or
whether a different ordering constraint prevents that. If that doesn't
work, defining the derived type in a separate module, which you USE
here, would work. Derived types are allowed in implicit statements, I'm
sure.

However, I can't really recommend that approach. The well-known
error-proneness of implicit typing is greatly magnified in the presence
of derived types, modules, and host association.

James Van Buskirk

unread,
Apr 21, 2008, 12:43:05 PM4/21/08
to
"relaxmike" <michael...@gmail.com> wrote in message
news:ffd371af-6da7-44b4...@c19g2000prf.googlegroups.com...

> I must say that the method based on implicit statements is
> very clever, and, based only on fortran statements, allows
> to minimize the code duplication. But is it possible to
> declare something like this :

> module data_mod
> type MYDATA
> character(len=20) :: string
> end type MYDATA
> implicit type(MYDATA) (Q)
> include 'test3_template.i90'
> end module data_mod

> I don't think so, which shows that the method cannot be extended
> to abstract data types. But funny though.

The method is indeed clever. Wish I had thought of it myself,
but someone else suggested it in clf.

As Richard pointed out of course you can do this:

C:\gfortran\clf\template_war>type test3_template.i90
private
public readfile_data
interface readfile_data
module procedure READFILE_DATA_NAME
end interface readfile_data
contains
subroutine READFILE_DATA_NAME(current_data_line,Qmyvalue,status)
character(len=*), intent(in) :: current_data_line
intent (out) :: Qmyvalue
integer, intent(out) :: status
status = 0
read(current_data_line, *, err=100 , end=100) Qmyvalue
return
100 continue
status = 1
end subroutine READFILE_DATA_NAME

C:\gfortran\clf\template_war>type test7.f90
module data_mod
implicit none


type MYDATA
character(len=20) string
end type MYDATA

end module data_mod

module logical_mod
implicit logical (Q)
include 'test3_template.i90'
end module logical_mod

module integer_mod
implicit integer (Q)
include 'test3_template.i90'
end module integer_mod

module MYDATA_mod
use data_mod


implicit type(MYDATA) (Q)
include 'test3_template.i90'

end module MYDATA_mod

module m_moduletest3
use logical_mod
use integer_mod

use MYDATA_mod
end module m_moduletest3

program test3
use data_mod


use m_moduletest3
character (len=200) :: current_data_line
integer :: myintvalue1
logical :: mylogicalvalue1

type(MYDATA) :: myMYDATAvalue1


integer :: status
current_data_line = "1"
call readfile_data(current_data_line, myintvalue1, status)
write(*, *) "status:", status
write(*, *) "Integer : ", myintvalue1
current_data_line = ".true."
call readfile_data(current_data_line, mylogicalvalue1, status )
write(*, *) "status:", status
write(*, *) "Logical : ", mylogicalvalue1

current_data_line = "'Hello from James'"
call readfile_data(current_data_line, myMYDATAvalue1, status )


write(*, *) "status:", status

write(*, *) "MYDATA : ", myMYDATAvalue1
end program test3

C:\gfortran\clf\template_war>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
te
st7.f90 -otest7

C:\gfortran\clf\template_war>test7


status: 0
Integer : 1
status: 0
Logical : T

status: 0
MYDATA : Hello from James

James Van Buskirk

unread,
Apr 21, 2008, 1:00:38 PM4/21/08
to
"Richard Maine" <nos...@see.signature> wrote in message
news:1ifqlj2.1nemrlg1h3158cN%nos...@see.signature...

> relaxmike <michael...@gmail.com> wrote:

>> I must say that the method based on implicit statements is
>> very clever, and, based only on fortran statements, allows
>> to minimize the code duplication. But is it possible to
>> declare something like this :

>> module data_mod
>> type MYDATA
>> character(len=20) :: string
>> end type MYDATA
>> implicit type(MYDATA) (Q)
>> include 'test3_template.i90'
>> end module data_mod

>> I don't think so, which shows that the method cannot be extended
>> to abstract data types.

> Yes, it is possible to do something like that. The only problem with the
> above is one of ordering. Implicit statements have to come quite early
> in a scoping unit, before pretty much anything other than USE
> statements. I forget whether it would be ok to just swap the order or
> whether a different ordering constraint prevents that. If that doesn't
> work, defining the derived type in a separate module, which you USE
> here, would work. Derived types are allowed in implicit statements, I'm
> sure.

The problem of ordering should not arise in practice because the
normal thing to do given what is forced upon the user by F90 rules is
to declare the user-defined type in one module along with assignment(=)
and any operators or procedures you want to overload for the type.
Then this module is USEd in the module that defines the template
procedures, which after all need assignment(=) and operators and
procedures overloaded for the type in the body of any template
procedures.

> However, I can't really recommend that approach. The well-known
> error-proneness of implicit typing is greatly magnified in the presence
> of derived types, modules, and host association.

This objection shows that you have never done template programming
in C++. It turns out to be so gawdawful that the normal way to
develop template code is to start with a non-template class and
methods, get that working OK, then convert to templates. Working
in an analogous fashion in Fortran, the non-template version would
have IMPLICIT NONE in force and only when it's working OK would
implicit typing be introduced. Besides IMPLICIT NONE won't help
with host association until the committee gives us some way to
restrict host association.

The really ugly part is: what are you supposed to do when you
need more than 26 template types? Maybe a Chinese Fortran
compiler would become popular at this point. Code with > 26
template types would be harder to read than Chinese, for sure.

glen herrmannsfeldt

unread,
Apr 21, 2008, 4:01:41 PM4/21/08
to
James Van Buskirk wrote:
(snip)

> The really ugly part is: what are you supposed to do when you
> need more than 26 template types? Maybe a Chinese Fortran
> compiler would become popular at this point. Code with > 26
> template types would be harder to read than Chinese, for sure.

Java has all the alphabetic unicode characters legal for
variable names. Confusingly, there are some that look exactly
like normal letters but are different places in the unicode
character set.

I would probably keep one letter each for integer and real.

-- glen

Craig Powers

unread,
Apr 21, 2008, 4:16:07 PM4/21/08
to
glen herrmannsfeldt wrote:
> James Van Buskirk wrote:
> (snip)
>
>> The really ugly part is: what are you supposed to do when you
>> need more than 26 template types? Maybe a Chinese Fortran
>> compiler would become popular at this point. Code with > 26
>> template types would be harder to read than Chinese, for sure.
>
> Java has all the alphabetic unicode characters legal for
> variable names. Confusingly, there are some that look exactly
> like normal letters but are different places in the unicode
> character set.

That could take obfuscated code to a whole new level.

relaxmike

unread,
Apr 22, 2008, 8:08:34 AM4/22/08
to
On 20 avr, 21:10, Charles Coldwell <coldw...@gmail.com> wrote:
> You forgot "callbacks":
>
> ! Fortran 95 implementation of the quicksort algorithm. This
> ! subroutine does not directly touch the array it sorts; rather it
> ! relies on the two callbacks "compare" and "exchange" for that. Code
> ! inspired by R. Sedgewick, "Algorithms in C" and correspondence with
> ! Glen Herrmannsfeldt.
>
> subroutine quicksort(n, compare, exchange)
> implicit none
> integer, intent(in) :: n ! the length of the implied array

Callbacks are a powerful way of uncoupling the sorting
algorithm from the type of data to sort, that is true.
That moves the data type from the sorting method
to the client code, where you have to implement a "quicksort"
method which is based on a particular type of data and pass the
type-specific "compare" method to the generic quicksort method.

The final implementation with callbacks is a move of
the code duplication from the sorting algorithm to the
client of the sorting algorithm, because, in the end,
you need a method to sort all types of arrays of data :
there is a benefit to use the quicksort algorithm
the way you showed it, but the problem appears in the client
of the quicksort.

My point of view is that callbacks cannot be used to manage
templates, but are a powerful way to configure a method for a
particular processing. In the sorting example, one can imagine
several "compare" methods to sort a list of strings :
- compare strings so that the sort is done in increasing order,
- compare strings so that the sort is done in decreasing order,
- compare strings so that the strings are converted to
integers and that the sort is done in increasing order,
- compare strings so that the strings are converted to
integers and that the sort is done in decreasing order,
- compare strings so that the strings are converted to
real and that the sort is done in increasing order,
- compare strings so that the strings are converted to
real and that the sort is done in decreasing order,
etc...

Regards,

Michaël

glen herrmannsfeldt

unread,
Apr 22, 2008, 2:42:12 PM4/22/08
to
relaxmike wrote:
(snip)

> Callbacks are a powerful way of uncoupling the sorting
> algorithm from the type of data to sort, that is true.
> That moves the data type from the sorting method
> to the client code, where you have to implement a "quicksort"
> method which is based on a particular type of data and pass the
> type-specific "compare" method to the generic quicksort method.

Not only that, but it is independent of the actual sort
algorithm. Rumor is that C's qsort() often isn't
quicksort despite the name.
(snip)

> My point of view is that callbacks cannot be used to manage
> templates, but are a powerful way to configure a method for a
> particular processing. In the sorting example, one can imagine
> several "compare" methods to sort a list of strings :

One I did some years ago had to compare strings with TeX
accents removed.

-- glen

Charles Coldwell

unread,
Apr 23, 2008, 6:06:43 AM4/23/08
to
relaxmike <michael...@gmail.com> writes:

> On 20 avr, 21:10, Charles Coldwell <coldw...@gmail.com> wrote:
>> You forgot "callbacks":
>>
>> ! Fortran 95 implementation of the quicksort algorithm. This
>> ! subroutine does not directly touch the array it sorts; rather it
>> ! relies on the two callbacks "compare" and "exchange" for that. Code
>> ! inspired by R. Sedgewick, "Algorithms in C" and correspondence with
>> ! Glen Herrmannsfeldt.
>>
>> subroutine quicksort(n, compare, exchange)
>> implicit none
>> integer, intent(in) :: n ! the length of the implied array
>
> Callbacks are a powerful way of uncoupling the sorting
> algorithm from the type of data to sort, that is true.
> That moves the data type from the sorting method
> to the client code, where you have to implement a "quicksort"
> method which is based on a particular type of data and pass the
> type-specific "compare" method to the generic quicksort method.

I don't understand. Let me repost the top of the function above:

subroutine quicksort(n, compare, exchange)
implicit none
integer, intent(in) :: n ! the length of the implied array

! The compare function must return an integer that is


! greater than zero if element(i) > element(j)
! equal to zero if element(i) = element(j)
! less than zero if element(i) < element(j)
interface
integer function compare(i,j)
integer, intent(in) :: i, j
end function compare
end interface

! The exchange subroutine exchanges the elements at locations i and j.
interface
subroutine exchange(i,j)
integer, intent(in) :: i, j
end subroutine exchange
end interface

Client code only implements "compare" and "exchange". So a generic
"sortable" type is any type that implements these. Here's how to make
reals "sortable". Think of this as doing something like a template
specialization for sortable<real> in C++:

module sortable_real

real, dimension(:), allocatable :: data

contains

integer function compare(i, j)
integer, intent(in) :: i, j
real :: r

r = data(i) - data(j)
if (r .eq. 0.0) then
compare = 0
else
compare = 2*r - 1
end if
end function compare

subroutine exchange(i, j)
integer, intent(in) :: i, j
real :: tmp

tmp = data(i)
data(i) = data(j)
data(j) = tmp
end subroutine exchange
end module sortable_real

program qsort_example
use sortable_real
integer, parameter :: n = 100

allocate(data(n))

! Do something to fill data

call quicksort(n, compare, exchange)
end program qsort_example

Obviously, one could also write template specializations for integers
and strings along the same lines, and use the same, identical
quicksort subroutine to sort them. So, a generic sortable in Fortran
is any type for which you can write "compare" and "exchange"
callbacks. Note that the dummy arguments passed to these functions
are always two integers, irregardless of the type of data being
sorted.

> The final implementation with callbacks is a move of the code
> duplication from the sorting algorithm to the client of the sorting
> algorithm, because, in the end, you need a method to sort all types
> of arrays of data : there is a benefit to use the quicksort
> algorithm the way you showed it, but the problem appears in the
> client of the quicksort.

I don't see any code duplication in the client. The implementations
of "compare" and "exchange" must be different for different types.
It's like in C++, where "operator<" (the analog of "compare") must be
defined for any type in an STL container.

Ken.Fa...@gmail.com

unread,
Apr 23, 2008, 12:03:50 PM4/23/08
to
One small nit:

On Apr 23, 3:06 am, Charles Coldwell <coldw...@gmail.com> wrote:
[...]


> module sortable_real
>
> real, dimension(:), allocatable :: data
>
> contains
>
> integer function compare(i, j)
> integer, intent(in) :: i, j
> real :: r
>
> r = data(i) - data(j)
> if (r .eq. 0.0) then
> compare = 0
> else
> compare = 2*r - 1
> end if
> end function compare

[...]

"compare = 2*r - 1" certainly makes some assumptions
about the values of data(i), eh? A typo? I don't understand
hoe that would work at all in this case.

Perhaps something more along the lines of:

compare = Int(Sign(1.0, r))

Granted, the Int(...) is not required in this assignment,
but may clarify to the reader what is desired...

-Ken

Tobias Burnus

unread,
Apr 23, 2008, 1:22:13 PM4/23/08
to
On Apr 23, 12:06 pm, Charles Coldwell <coldw...@gmail.com> wrote:
> subroutine quicksort(n, compare, exchange)

>   integer, intent(in) :: n ! the length of the implied array
>   interface
>      integer function compare(i,j)
>        integer, intent(in) :: i, j

If I see this (and not completely remembering the initial question),
I'm wondering whether one could not implement this using something
like the following (I'm sure I got the syntax wrong somewhere). This
uses Fortran 2003's OOP, which most compilers have presumably not
implemented, and is completely untested. Whether this code is elegant
is in the eye of the beholder.

Tobias


module sort
abstract interface
integer function comp_func(i,j, this)
integer :: i, j
type(sort_type) :: this
end function comp_func
subroutine exch_func(i,j, this)
integer :: i, j
type(sort_type) :: this
end subroutine exch_f
end interface
type, abstract :: sortType
integer :: n
contains
procedure(comp_func), deferred :: compare
procedure(exch_func), deferred :: exchange
end type
contains
subroutine quicksort(data)
class(sortType) :: data
! ...
if (data%compare(i,j) > 0) call data%exchange(i,j)
! ...
end subroutine quicksort
end module sort

module myDataType
use sort
type, extends(sortType) :: myType
private
integer, allocatable :: i()
contains
procedure set
procedure compare
procedure exchange
end type

contains
subroutine set(array, this)
integer :: array(:)
type(myType) :: this
this%i = data
this%n = size(array)
end subroutine set
subroutine exchange(i, j, this)
integer :: i, j
type(myType) :: this
integer :: tmp
tmp = this%i(i)
this%i(i) = this%i(j)
this%i(j) = tmp
end subroutine exchange
! ...
end myDataType

Charles Coldwell

unread,
Apr 23, 2008, 5:21:28 PM4/23/08
to
Ken.Fa...@gmail.com writes:

>
> "compare = 2*r - 1" certainly makes some assumptions
> about the values of data(i), eh? A typo? I don't understand
> hoe that would work at all in this case.

It wouldn't. More like a brain-o.

In C, I would write this

int compare(int i, int j)
{
double r = data[i] - data[j];

return (2*(r>0) - 1) & -(r != 0);
}

This is a clever thing you can do in C because there is no distinction
between LOGICAL and INTEGER. So if (r>0) is true, it has the integer
value 1, then (2*(r>0) - 1) has the integer values plus or minus 1,
and finally to deal with the case of r equal to zero the final
*bitwise* and with -(r != 0). If r is equal to zero, (r != 0) is
false, or integer value 0, thus -(r != 0) is still 0 and bitwise and 0
with anything give zero back. If (r != 0) is true, then -(r != 0)
equals integer value -1, and bitwise and -1 with anything (on 2s
complement arch) and you get the same thing back.

Very clever, not translatable into Fortran. In Fortran, the right
thing to do is

integer function compare(i, j)
integer, intent(in) :: i, j

if (data(i) .lt. data(j)) then
compare = -1
else if (data(i) .gt. data(j)) then
compare = 1
else
compare = 0
end function compare

Sorry about the silliness.

glen herrmannsfeldt

unread,
Apr 23, 2008, 5:56:00 PM4/23/08
to
Charles Coldwell wrote:
(snip)

> In C, I would write this

> int compare(int i, int j)
> {
> double r = data[i] - data[j];
> return (2*(r>0) - 1) & -(r != 0);
> }

I think I more often see:

return (r>0)-(r<0);

For character compare it is often done as:

return i-j;

(or, for a string compare, subtract the character values
of the first character that doesn't compare equal.)

Otherwise, how about

compare=count((/ data[i].ge.data[j], data[i].gt.data[j] /))-1

-- glen

Charles Coldwell

unread,
Apr 23, 2008, 6:43:29 PM4/23/08
to
glen herrmannsfeldt <g...@ugcs.caltech.edu> writes:

> Charles Coldwell wrote:
> (snip)
>
>> In C, I would write this
>
>> int compare(int i, int j)
>> {
>> double r = data[i] - data[j];
>> return (2*(r>0) - 1) & -(r != 0);
>> }
>
> I think I more often see:
>
> return (r>0)-(r<0);

Oooh. I like that.

> Otherwise, how about
>
> compare=count((/ data[i].ge.data[j], data[i].gt.data[j] /))-1

Very cute, I like it, too. I have more confidence in the compiler's
ability to generate tight code for the C version, above, though.

relaxmike

unread,
Apr 24, 2008, 4:01:37 AM4/24/08
to
On 23 avr, 12:06, Charles Coldwell <coldw...@gmail.com> wrote:
> I don't see any code duplication in the client. The implementations
> of "compare" and "exchange" must be different for different types.
> It's like in C++, where "operator<" (the analog of "compare") must be
> defined for any type in an STL container.

I think I now understand what you mean, and let me explain what
I mean : there is no duplication in the client, simply because there
is only one client, of course ! This is because you transformed
the problem of "having a method which sorts whatever type of array"
into "having a method which sorts the client-side array".
The consequence is that you define only one "compare" method, and
you conclude that there is no duplication.

Now let's define what should be to do to have a method to
sort arrays of logical, integers, real and double :
you, of course, have to define a "compare_logical" ,
"compare_integer",
"compare_real", "compare_double". This is because the compare method
have to explicitly define the type of data to compare.
In the end, there are 4 implementations of the same "template"
compare, all using the same (and basic) "<".
That is true that if you want to compare characters, you may have to
define "<" (but I am not sure).
Isn'it duplication ?

I looked in the internet to search for how the C++ STL solved that
problem. I found a web page, but I do not known if it is of high
quality :
http://www.cppreference.com/cppalgorithm/sort.html
See how the sort of an array of integers is solved :

vector<int> v;
v.push_back( 23 );
v.push_back( -1 );
v.push_back( 9999 );
v.push_back( 0 );
v.push_back( 4 );
sort( v.begin(), v.end() );

I don't pretend that fortran can achieve such a level of
simplicity, but that's the goal isn't it ?
I think that callbacks are a good and efficient way
to uncouple the sorting algorithm from the type of
data to sort, but that is only half of the work to do
to have a routine which can sort <whatever> type of
array.

Regards,

Michaël

Greg Lindahl

unread,
Apr 24, 2008, 4:14:24 AM4/24/08
to
In article <9b4f7b86-02ae-4717...@l64g2000hse.googlegroups.com>,
relaxmike <michael...@gmail.com> wrote:

>See how the sort of an array of integers is solved :
>
> vector<int> v;
> v.push_back( 23 );
> v.push_back( -1 );
> v.push_back( 9999 );
> v.push_back( 0 );
> v.push_back( 4 );
> sort( v.begin(), v.end() );
>
>I don't pretend that fortran can achieve such a level of
>simplicity, but that's the goal isn't it ?

You've got to be kidding! Variables that lack a constructor syntax are
not that simple. And the STL is well known to have a pretty steep
learning curve.

-- greg


Charles Coldwell

unread,
Apr 24, 2008, 6:36:44 AM4/24/08
to
relaxmike <michael...@gmail.com> writes:

> On 23 avr, 12:06, Charles Coldwell <coldw...@gmail.com> wrote:
>> I don't see any code duplication in the client. The implementations
>> of "compare" and "exchange" must be different for different types.
>> It's like in C++, where "operator<" (the analog of "compare") must be
>> defined for any type in an STL container.
>
> I think I now understand what you mean, and let me explain what
> I mean : there is no duplication in the client, simply because there
> is only one client, of course ! This is because you transformed
> the problem of "having a method which sorts whatever type of array"
> into "having a method which sorts the client-side array".
> The consequence is that you define only one "compare" method, and
> you conclude that there is no duplication.
>
> Now let's define what should be to do to have a method to
> sort arrays of logical, integers, real and double :
> you, of course, have to define a "compare_logical" ,
> "compare_integer",
> "compare_real", "compare_double". This is because the compare method
> have to explicitly define the type of data to compare.
> In the end, there are 4 implementations of the same "template"
> compare, all using the same (and basic) "<".
> That is true that if you want to compare characters, you may have to
> define "<" (but I am not sure).
> Isn'it duplication ?

Yes and no. It's more duplication than what you find in C++. For
example, to use an ordered STL container, a class must define
"operator<". For POD types (int, long, double, float), you can use
the built-in compiler provided operator, you don't have to write any
code, whereas in my Fortran quicksort-with-callbacks, you must define
"compare", even for the built-in types (e.g. INTEGER, REAL, LOGICAL).
However, if you want to sort *strings*, you will have to define
"operator<" for the string class. Similarly for any user-defined
types. Now, admittedly, "operator<" for std::string comes with the
standard library, so you would be crazy to actually write it yourself,
but you still have to define it for all the types you design yourself.

> I looked in the internet to search for how the C++ STL solved that
> problem. I found a web page, but I do not known if it is of high
> quality :
> http://www.cppreference.com/cppalgorithm/sort.html
> See how the sort of an array of integers is solved :
>
> vector<int> v;
> v.push_back( 23 );
> v.push_back( -1 );
> v.push_back( 9999 );
> v.push_back( 0 );
> v.push_back( 4 );
> sort( v.begin(), v.end() );
>
> I don't pretend that fortran can achieve such a level of
> simplicity, but that's the goal isn't it ?

Sure. What I wrote is a Fortran function that implements the "sort"
that comes with the C++ standard library. You could argue that
Fortran should supply a sort intrinsic so that the code isn't
duplicated by many implementors, but that's different from generic
programming. In terms of simplicity, if I wrote a module implementing
compare and exchange for integers, the C++ code above in Fortran would
read

use my_module
data(1) = 23
data(2) = -1
data(3) = 9999
data(4) = 0
data(5) = 4
call quicksort(5, compare, exchange)

The difference with C++ in this example is that the standard library
is supplying the sort routine, the vector<> template and it's
associated "operator<" (which in this case is just the built-in
"operator<" for type int). The quicksort routine I posted previously
is almost as general as the one that the C++ standard library
supplies: it can work on an array of any type, including user-defined
types, hence my claim that it demonstrates "generic programming" in
Fortran. The C++ STL sort is more generic because it works on
containers of any type, including lists, deques, arrays, vectors, maps
and multimaps.

> I think that callbacks are a good and efficient way
> to uncouple the sorting algorithm from the type of
> data to sort,

Hence my claim that it is "generic".

> but that is only half of the work to do to have a routine which can
> sort <whatever> type of array.

Fair enough, but at least for user-defined types, you have to do the
other half of the work yourself whether you use C++ or Fortran.

relaxmike

unread,
Apr 24, 2008, 7:28:01 AM4/24/08
to
On 24 avr, 12:36, Charles Coldwell <coldw...@gmail.com> wrote:
> Yes and no. It's more duplication than what you find in C++. For
> example, to use an ordered STL container, a class must define
> "operator<". For POD types (int, long, double, float), you can use
> the built-in compiler provided operator, you don't have to write any
> code, whereas in my Fortran quicksort-with-callbacks, you must define
> "compare", even for the built-in types (e.g. INTEGER, REAL, LOGICAL).
> However, if you want to sort *strings*, you will have to define
> "operator<" for the string class. Similarly for any user-defined
> types. Now, admittedly, "operator<" for std::string comes with the
> standard library, so you would be crazy to actually write it yourself,
> but you still have to define it for all the types you design yourself.

Yes of course, you are right.
I cannot consider that defining three compare routines is such a
big deal that one should use a template for that.
We agree that the callback method is more closed to generic
programming
that it is to templates.

> The difference with C++ in this example is that the standard library
> is supplying the sort routine, the vector<> template and it's
> associated "operator<" (which in this case is just the built-in
> "operator<" for type int). The quicksort routine I posted previously
> is almost as general as the one that the C++ standard library
> supplies: it can work on an array of any type, including user-defined
> types, hence my claim that it demonstrates "generic programming" in
> Fortran. The C++ STL sort is more generic because it works on
> containers of any type, including lists, deques, arrays, vectors, maps
> and multimaps.

Obviously, I don't know the STL sufficiently to compare fortran and
C++ at that level. And I do not master the C++ iterators.

> Fair enough, but at least for user-defined types, you have to do the
> other half of the work yourself whether you use C++ or Fortran.

You are right and there is no other way to do it.
Anyway, the algorithm you provided is an interesting solution for the
generic sorting problem.

Michaël

James Van Buskirk

unread,
Apr 24, 2008, 1:33:13 PM4/24/08
to
"Charles Coldwell" <cold...@gmail.com> wrote in message
news:rzpy773...@gmail.com...

> Sure. What I wrote is a Fortran function that implements the "sort"
> that comes with the C++ standard library. You could argue that
> Fortran should supply a sort intrinsic so that the code isn't
> duplicated by many implementors, but that's different from generic
> programming. In terms of simplicity, if I wrote a module implementing
> compare and exchange for integers, the C++ code above in Fortran would
> read

> use my_module
> data(1) = 23
> data(2) = -1
> data(3) = 9999
> data(4) = 0
> data(5) = 4
> call quicksort(5, compare, exchange)


What I don't understand is why you're always passing around these
compare and exchange procedures. In Fortran, wouldn't it be more
normal to do something like:

C:\gfortran\clf\qsort>type qsort.i90
recursive subroutine qsort_sub(Qarray)
dimension Qarray(0:)
real harvest
integer pivot
integer left
integer right

if(size(Qarray) < 2) return
call random_number(harvest)
pivot = min(int(harvest*size(Qarray)),ubound(Qarray,1))
if(pivot > 0) Qarray(0:pivot:pivot) = Qarray(pivot:0:-pivot)
left = 1
right = ubound(Qarray,1)
do while(left <= right)
do
if(left > right) exit
if(Qarray(0) < Qarray(left)) exit
left = left+1
end do
do
if(left > right) exit
if(.NOT.(Qarray(0) < Qarray(right))) exit
right = right-1
end do
if(left < right) then
Qarray(left:right:right-left) = Qarray(right:left:left-right)
left = left+1
right = right-1
end if
end do
if(left > 1) Qarray(0:left-1:left-1) = Qarray(left-1:0:1-left)
call qsort(Qarray(0:left-2))
call qsort(Qarray(left:))
end subroutine qsort_sub

C:\gfortran\clf\qsort>type qtest.f90
module mytype_define
implicit none
private
public mytype, assignment(=), operator(<)
type mytype
integer data
end type mytype
interface assignment(=)
module procedure assign
end interface assignment(=)
interface operator(<)
module procedure less_than
end interface operator(<)
contains
elemental subroutine assign(x,y)
type(mytype), intent(out) :: x
type(mytype), intent(in) :: y

x%data = y%data
end subroutine assign
elemental function less_than(x,y)
type(mytype), intent(in) :: x
type(mytype), intent(in) :: y
logical less_than

less_than = x%data < y%data
end function less_than
end module mytype_define

module i4_mod
implicit integer(Q)
interface qsort
module procedure qsort_sub
end interface qsort
contains
include 'qsort.i90'
end module i4_mod

module r4_mod
implicit real(Q)
interface qsort
module procedure qsort_sub
end interface qsort
contains
include 'qsort.i90'
end module r4_mod

module mytype_mod
use mytype_define
implicit type(mytype) (Q)
interface qsort
module procedure qsort_sub
end interface qsort
contains
include 'qsort.i90'
end module mytype_mod

module generic_sort
use i4_mod, only: qsort
use r4_mod, only: qsort
use mytype_mod, only: qsort
implicit none
end module generic_sort

program test
use generic_sort
use mytype_define
implicit none
integer, allocatable :: i4x(:)
real, allocatable :: r4x(:)
type(mytype), allocatable :: mtx(:)
integer n
integer i
real harvest

call random_seed()
n = 10
allocate(i4x(n))
allocate(r4x(n))
allocate(mtx(n))
do i = 1, n
call random_number(harvest)
i4x(i) = harvest*n**2
call random_number(harvest)
r4x(i) = harvest
call random_number(harvest)
mtx(i) = mytype(harvest*n**2)
end do
call qsort(i4x)
call qsort(r4x)
call qsort(mtx)
do i = 1, n
write(*,*) i, i4x(i), r4x(i), mtx(i)
end do
end program test

C:\gfortran\clf\qsort>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
qtest.f90 -
oqtest

C:\gfortran\clf\qsort>qtest

C:\gfortran\clf\qsort>ifort qtest.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.

-out:qtest.exe
-subsystem:console
qtest.obj

C:\gfortran\clf\qsort>qtest
1 5 6.2769456E-03 10
2 14 3.7755325E-02 10
3 26 0.1845653 10
4 31 0.2239807 10
5 33 0.4793749 53
6 35 0.6826791 10
7 48 0.8232406 10
8 62 0.9100211 10
9 77 0.9194745 53
10 97 0.9970329 61

C:\gfortran\clf\qsort>ftn95 qtest.f90 /link
[FTN95/Win32 Ver. 5.0.0 Copyright (c) Silverfrost Ltd 1993-2006]
PROCESSING MODULE [<MYTYPE_DEFINE> FTN95/Win32 v5.0.0]
NO ERRORS [<ASSIGN> FTN95/Win32 v5.0.0]
NO ERRORS [<LESS_THAN> FTN95/Win32 v5.0.0]
NO ERRORS [<MYTYPE_DEFINE> FTN95/Win32 v5.0.0]
PROCESSING MODULE [<I4_MOD> FTN95/Win32 v5.0.0]
NO ERRORS [<QSORT_SUB> FTN95/Win32 v5.0.0]
NO ERRORS [<I4_MOD> FTN95/Win32 v5.0.0]
PROCESSING MODULE [<R4_MOD> FTN95/Win32 v5.0.0]
NO ERRORS [<QSORT_SUB> FTN95/Win32 v5.0.0]
NO ERRORS [<R4_MOD> FTN95/Win32 v5.0.0]
PROCESSING MODULE [<MYTYPE_MOD> FTN95/Win32 v5.0.0]
NO ERRORS [<QSORT_SUB> FTN95/Win32 v5.0.0]
NO ERRORS [<MYTYPE_MOD> FTN95/Win32 v5.0.0]
PROCESSING MODULE [<GENERIC_SORT> FTN95/Win32 v5.0.0]
NO ERRORS [<GENERIC_SORT> FTN95/Win32 v5.0.0]
0066) use generic_sort
*** Specific procedure QSORT_SUB of type SUBROUTINE is too similar to
QSORT_SUB of type SUBROUTINE for overload QSORT in module R4_MOD. The
arguments are too similar.
0091) call qsort(mtx)
*** No matching specific procedure for generic overloaded name QSORT
2 ERRORS [<TEST> FTN95/Win32 v5.0.0]
*** Compilation failed

See how much more simple? And works great!

James Van Buskirk

unread,
Apr 24, 2008, 6:19:33 PM4/24/08
to
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:-_-dnbCWBbjHWY3V...@comcast.com...

> And works great!

Well, I tried one small change:

! public mytype, assignment(=), operator(<)
public mytype, operator(<)!, assignment(=)

And named the modified file qtest1.f90 .

C:\gfortran\clf\qsort>ifort qtest.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.

-out:qtest.exe
-subsystem:console
qtest.obj

C:\gfortran\clf\qsort>qtest

1 0 2.0906582E-02 13
2 0 8.1042856E-02 13
3 8 8.6839773E-02 15
4 25 0.1941290 15
5 33 0.2825256 15
6 57 0.3655469 13
7 65 0.4545409 36
8 75 0.4681000 36
9 84 0.4968905 33
10 93 0.5864838 33

C:\gfortran\clf\qsort>ifort qtest1.f90


Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.

-out:qtest1.exe
-subsystem:console
qtest1.obj

C:\gfortran\clf\qsort>qtest1
1 14 0.2617609 9
2 28 0.3624567 9
3 31 0.3778799 11
4 34 0.4567579 19
5 64 0.4780824 20
6 73 0.5333692 24
7 79 0.5638247 36
8 89 0.6962512 41
9 98 0.7379416 72
10 98 0.8467951 77

C:\gfortran\clf\qsort>C:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran -v
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:
../../trunk/configure --prefix=/home/FX/irun64 --build=i586-pc-
mingw32 --target=x86_64-pc-mingw32 --with-gmp=/home/FX/local --enable-languages=
c,fortran --disable-werror --disable-nls --enable-threads
Thread model: win32
gcc version 4.4.0 20080421 (experimental) [trunk revision 134506] (GCC)

C:\gfortran\clf\qsort>C:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
qtest.f90
-oqtest

C:\gfortran\clf\qsort>qtest

C:\gfortran\clf\qsort>C:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
qtest1.f9
0 -oqtest1

C:\gfortran\clf\qsort>qtest1
1 7 5.35517931E-03 13
2 33 1.61082745E-02 34
3 34 0.10038292 44
4 40 0.20687431 45
5 59 0.21795171 48
6 64 0.32298726 65
7 66 0.36739087 75
8 74 0.38676596 85
9 90 0.56682467 96
10 99 0.67298073 96

C:\gfortran\clf\qsort>C:\gcc_equation\bin\x86_64-pc-mingw32-gfortran -v
Built by Equation Solution (http://www.Equation.com).
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:
../gcc-4.4-20080418-mingw/configure --host=x86_64-pc-mingw32 --
build=x86_64-unknown-linux-gnu --target=x86_64-pc-mingw32 --prefix=/home/gfortra
n/gcc-home/binary/mingw32/native/x86_64/gcc/4.4-20080418 --with-gmp=/home/gfortr
an/gcc-home/binary/mingw32/native/x86_64/gmp --with-mpfr=/home/gfortran/gcc-home
/binary/mingw32/native/x86_64/mpfr --with-sysroot=/home/gfortran/gcc-home/binary
/mingw32/cross/x86_64/gcc/4.4-20080418 --with-gcc --with-gnu-ld --with-gnu-as
--
disable-shared --disable-nls --disable-tls --enable-languages=c,fortran --enable
-threads=win32 --enable-libgomp --disable-win32-registry
Thread model: win32
gcc version 4.4.0 20080418 (experimental) (GCC)

C:\gfortran\clf\qsort>C:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
qtest.f90 -
oqtest

C:\gfortran\clf\qsort>qtest

C:\gfortran\clf\qsort>C:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
qtest1.f90
-oqtest1

C:\gfortran\clf\qsort>qtest1
1 7 5.35517931E-03 13
2 33 1.61082745E-02 34
3 34 0.10038292 44
4 40 0.20687431 45
5 59 0.21795171 48
6 64 0.32298726 65
7 66 0.36739087 75
8 74 0.38676596 85
9 90 0.56682467 96
10 99 0.67298073 96

So we can see that old ifort and the latest versions of gfortran
(both styles) can cope if intrinsic rather than defined
assignment is used. Kind of strange because gfortran can
handle
http://home.comcast.net/~kmbtib/Fortran_stuff/elem_assign.f90
but this stuff with array triplets rather than vector subscripts
is too much for it. Also one might have thought that ifort
would have gotten it fixed after
http://groups.google.com/group/comp.lang.fortran/msg/a76a9ebd13bd04ab
but even after issue #345003 +590 days I can't state with confidence
that it has been fixed. I sure haven't been notified if so.
Salford/Silverfrost barfs at the same place in the modified version.

Charles Coldwell

unread,
Apr 24, 2008, 8:12:36 PM4/24/08
to
"James Van Buskirk" <not_...@comcast.net> writes:

> "Charles Coldwell" <cold...@gmail.com> wrote in message
> news:rzpy773...@gmail.com...
>
>> Sure. What I wrote is a Fortran function that implements the "sort"
>> that comes with the C++ standard library. You could argue that
>> Fortran should supply a sort intrinsic so that the code isn't
>> duplicated by many implementors, but that's different from generic
>> programming. In terms of simplicity, if I wrote a module implementing
>> compare and exchange for integers, the C++ code above in Fortran would
>> read
>
>> use my_module
>> data(1) = 23
>> data(2) = -1
>> data(3) = 9999
>> data(4) = 0
>> data(5) = 4
>> call quicksort(5, compare, exchange)
>
>
> What I don't understand is why you're always passing around these
> compare and exchange procedures. In Fortran, wouldn't it be more
> normal to do something like:

[ snip ]

Well, if you don't object to implicit typing, it's a pretty clever
solution, although the "includes" sort of imply a certain amount of
code duplication, although I guess I wouldn't really call it that.
What you are doing here is very close to what would be called
"template instantiation" in C++, and like template instantiation, your
includes are evaluated at compile time.

> 0091) call qsort(mtx)
> *** No matching specific procedure for generic overloaded name QSORT
> 2 ERRORS [<TEST> FTN95/Win32 v5.0.0]
> *** Compilation failed

Heh. You are exploring the dusty corners of the F95 spec. It may be
conforming code, but you have confused at least one compiler.

> See how much more simple? And works great!

Chip

Jim Xia

unread,
Apr 30, 2008, 11:07:37 AM4/30/08
to
The correct syntax for your code should be something like this

module sort


type, abstract :: sortType
integer :: n
contains

procedure(comp_func), deferred, pass(this) :: compare
procedure(exch_func), deferred, pass(this) :: exchange
end type

abstract interface
integer function comp_func(i,j, this)

import
integer :: i, j
class(sorttype) :: this
end function comp_func

subroutine exch_func(i,j, this)
import
integer :: i, j
class(sorttype) :: this
end subroutine exch_func
end interface


contains
subroutine quicksort(data)
class(sortType) :: data
! ...
if (data%compare(i,j) > 0) call data%exchange(i,j)
! ...
end subroutine quicksort
end module sort


module myDataType
use sort
type, extends(sortType) :: myType
private

integer, allocatable :: i(:)
contains
procedure, pass(this) ::set
procedure, pass(this) ::compare
procedure, pass(this) ::exchange
end type


contains
subroutine set(array, this)
integer :: array(:)

class(myType) :: this


this%i = data
this%n = size(array)
end subroutine set
subroutine exchange(i, j, this)
integer :: i, j

class(myType) :: this


integer :: tmp
tmp = this%i(i)
this%i(i) = this%i(j)
this%i(j) = tmp
end subroutine exchange

integer function compare(i,j, this)
integer :: i, j
class(myType) :: this
! ...
end function compare
! ...
end module myDataType


If we want to introduce templates in Fortran, one way I see it is to
do it in the module procedures. The problem is not whether we can do
it, it's rather if there are enough demands for it.

Cheers,

Damian

unread,
Apr 30, 2008, 8:54:17 PM4/30/08
to

Jim,

The question of whether there is demand for templates is an
interesting one. As a parallel, consider my correspondence with one
Fortran vendor who indicated that the reason they didn't see the
implementation of the object-oriented features in Fortran as a high
priority is that few of their users had requested it. I commented
that might be because most people who wanted object-orientation
abandoned Fortran long ago. Possibly the same is true of templates.
Is there a role for the standards committee to provide leadership in
guiding the programming community toward new paradigms (generic
programming) as opposed to responding to existing demand? Possibly
leading is a bad idea compared to responding to demands, but I would
put co-arrays in the category of leadership, and I'm hoping co-array
Fortran eventually displace existing approaches to parallel
programming.

Damian

Damian

unread,
May 2, 2008, 12:55:50 AM5/2/08
to

Amongst the last 300 threads in this newsgroup, this thread ranks
approximately 6th in terms of the number of postings. That might be
one indicator of high interest in templates -- although I suppose it's
also a measure of how complicated the topic is!

Damian

Ron Ford

unread,
May 2, 2008, 1:17:32 AM5/2/08
to

"Damian" <dam...@rouson.net> wrote in message
news:80d7d71a-7e93-4c42...@q27g2000prf.googlegroups.com...

I'd like to see that script. The subjects in this ng are a curious
phenomenon.

I think the call for "templates" is identical to "generic programming." I
like the emphasis on algorithms. The only thing that excites me about C is
that Heathfield has a book on C algorithms now.

For me it all goes back to Knuth. Unbelievable that Condie Rice hails from
whatever alma doesn't matter.

--
"Life in Lubbock, Texas, taught me two things: One is that God loves you
and you're going to burn in hell. The other is that sex is the most
awful, filthy thing on earth and you should save it for someone you love."

~~ Butch Hancock


Jim Xia

unread,
May 2, 2008, 11:16:26 AM5/2/08
to
> Damian- Hide quoted text -
>
> - Show quoted text -


Damian

The Fortran committee recognizes the need for generic programming in
the language. There was an effort in 2008 standard to solve the
problem by MACRO. But that feature was voted out due to its
complexity and also the fact it is based on a rather outdated
technology. I personally dislike the MACRO approach at all: it's a
clumsy tool that generates very convoluted code -- the readability is
very poor. I think there must be a better way to achieve that.

Resistance to templates in Fortran is likely by the standard body
since it'll hurt performance due to code-bloat. Fortran is known to
generate fast code that C++ can even dream of. Any feature
endangering that capability is not going to be an easy sell. That
also explains why coarrays are welcomed to the language while the
templates are still being talked about.

If we shoot for a limited version of templates (as compared to C++)
that minimizes its damage to speed, then I think we have a chance.

Cheers,

Jim

Craig Powers

unread,
May 2, 2008, 1:47:58 PM5/2/08
to
Jim Xia wrote:
>
> The Fortran committee recognizes the need for generic programming in
> the language. There was an effort in 2008 standard to solve the
> problem by MACRO. But that feature was voted out due to its
> complexity and also the fact it is based on a rather outdated
> technology. I personally dislike the MACRO approach at all: it's a
> clumsy tool that generates very convoluted code -- the readability is
> very poor. I think there must be a better way to achieve that.
>
> Resistance to templates in Fortran is likely by the standard body
> since it'll hurt performance due to code-bloat. Fortran is known to
> generate fast code that C++ can even dream of. Any feature
> endangering that capability is not going to be an easy sell. That
> also explains why coarrays are welcomed to the language while the
> templates are still being talked about.
>
> If we shoot for a limited version of templates (as compared to C++)
> that minimizes its damage to speed, then I think we have a chance.

Bloat is an orthogonal problem to speed (which is not to say that it's
not something to be concerned about, only that it's largely independent
of performance), and as far as I know, C++ comes closest to Fortran for
performance in heavily templated code.

From a C++ point of view, my sense is that the biggest issues with
templates revolve around the complexity of the feature. I don't think
there's all that much concern with code size, as I don't think user code
typically sees more than one or two instantiations of a particular template.

James Van Buskirk

unread,
May 2, 2008, 5:33:31 PM5/2/08
to
"Jim Xia" <jim...@hotmail.com> wrote in message
news:c81e0a74-b847-4aa7...@m73g2000hsh.googlegroups.com...

> Resistance to templates in Fortran is likely by the standard body
> since it'll hurt performance due to code-bloat. Fortran is known to
> generate fast code that C++ can even dream of. Any feature
> endangering that capability is not going to be an easy sell. That
> also explains why coarrays are welcomed to the language while the
> templates are still being talked about.

Fortran has language features that C++ doesn't have that make
templates more awkward. One of my favorites is that of having 3
different classes of expressions: ordinary expressions,
initialization expressions, and specification expressions.
Consider the following example, where template code is used for
9 different instances: all 3 real KINDs times 3 different functions.
Since in f03 pretty much any intrinsic elemental function may take
part in an initialization epression, in our template code we must
be ready to handle named constants whose type is unknown to us as
we write the template code. I don't see how to handle this in
the general case, but if the type is constrained to be numeric it
can be done! Consider the following example:

C:\gfortran\test\num_init>type finit.s
.text
.globl _xfinit_
.def _xfinit_; .scl 2; .type 32; .endef
_xfinit_:
finit
ret

C:\gfortran\test\num_init>type num_init.i90
subroutine initr1_template(Q)
parameter(Qarg1 = x)
integer, parameter :: k2 = kind(fun(Qarg1))
integer, parameter :: is_int = 1-1/(2+0*fun(Qarg1))*2
integer, parameter :: kind_if_real = &
(1-is_int)*k2+is_int*kind(1.0)
integer, parameter :: is_cmplx = -((1-is_int)* &
transfer(cmplx(0,1,kind_if_real),fun(Qarg1)))**2
integer, parameter :: is_real = 1-is_int-is_cmplx
integer, parameter :: kind_if_int = is_int*k2+(1-is_int)*kind(1)
integer(kind_if_int), parameter :: iparam = is_int*fun(Qarg1)
real(kind_if_real), parameter :: rparam = is_real*fun(Qarg1)
complex(kind_if_real), parameter :: cparam = is_cmplx*fun(Qarg1)
character(7), parameter :: type_label_array(3) = &
['INTEGER','REAL ','COMPLEX']
character(*), parameter :: type_label = &
trim(type_label_array(is_int+2*is_real+3*is_cmplx))
character(100) fmt
integer s

write(*,'(a)') 'Test of '//name//' intrinsic'
write(*,'(a,i0)') 'k2 = ', k2
write(*,'(a,i0)') 'is_int = ', is_int
write(*,'(a,i0)') 'kind_if_real = ', kind_if_real
write(*,'(a,i0)') 'is_cmplx = ', is_cmplx
write(*,'(a,i0)') 'is_real = ', is_real
write(*,'(a,i0)') 'kind_if_int = ', kind_if_int
write(*,'(2a)') 'type_label = ', type_label
if(is_int == 1) then
write(*,'(a,i0)') name//'('//arg1//') = ', iparam
else if(is_real == 1) then
s = precision(rparam)
write(fmt,'(3(a,i0))') '(a,g',s+9,'.',s+2,')'
write(*,fmt) name//'('//arg1//') = ', rparam
else if(is_cmplx == 1) then
s = precision(cparam)
write(fmt,'(3(a,i0))') '(3(a,g',s+9,'.',s+2,'))'
write(*,fmt) name//'('//arg1//') = (', &
real(cparam), ',', aimag(cparam), ')'
end if
end subroutine initr1_template

C:\gfortran\test\num_init>type initr1.f90
module mykinds
implicit none
integer, parameter :: ck1 = kind('x')
integer, parameter :: ik1 = selected_int_kind(2)
integer, parameter :: ik2 = selected_int_kind(4)
integer, parameter :: ik4 = selected_int_kind(9)
integer, parameter :: ik8 = selected_int_kind(18)
integer, parameter :: Lk1 = ik1
integer, parameter :: Lk2 = ik2
integer, parameter :: Lk4 = ik4
integer, parameter :: Lk8 = ik8
integer, parameter :: sp = selected_real_kind(6,37)
integer, parameter :: dp = selected_real_kind(15,307)
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
type mytype
integer(ik4) data
end type mytype
end module mykinds

module name_mod
use mykinds, only: qp
private qp
intrinsic INT
intrinsic ABS
intrinsic CMPLX
!DEC$ IF DEFINED(.FALSE.)
intrinsic BESSEL_J0
!DEC$ ENDIF
intrinsic ERF
intrinsic SPACING
character(*), parameter :: int_name = 'INT'
character(*), parameter :: int_arg1 = 'X'
real(qp), parameter :: int_X = 1.795195802051310421978653361874_qp
character(*), parameter :: abs_name = 'ABS'
character(*), parameter :: abs_arg1 = 'X'
real(qp), parameter :: abs_X = 1.795195802051310421978653361874_qp
character(*), parameter :: cmplx_name = 'CMPLX'
character(*), parameter :: cmplx_arg1 = 'X'
real(qp), parameter :: cmplx_X = 1.795195802051310421978653361874_qp
end module name_mod

module int_mod
use name_mod, only: fun=> INT
use name_mod, only: name => int_name
use name_mod, only: arg1 => int_arg1
use name_mod, only: X => int_X
end module int_mod

module abs_mod
use name_mod, only: fun=> ABS
use name_mod, only: name => abs_name
use name_mod, only: arg1 => abs_arg1
use name_mod, only: X => abs_X
end module abs_mod

module cmplx_mod
use name_mod, only: fun=> CMPLX
use name_mod, only: name => cmplx_name
use name_mod, only: arg1 => cmplx_arg1
use name_mod, only: X => cmplx_X
end module cmplx_mod

module initr1_sp_int
use mykinds
use int_mod
implicit real(sp) (Q)
private
public initr1_int
interface initr1_int
module procedure initr1_template
end interface initr1_int
contains
include 'num_init.i90'
end module initr1_sp_int

module initr1_dp_int
use mykinds
use int_mod
implicit real(dp) (Q)
private
public initr1_int
interface initr1_int
module procedure initr1_template
end interface initr1_int
contains
include 'num_init.i90'
end module initr1_dp_int

module initr1_qp_int
use mykinds
use int_mod
implicit real(qp) (Q)
private
public initr1_int
interface initr1_int
module procedure initr1_template
end interface initr1_int
contains
include 'num_init.i90'
end module initr1_qp_int

module initr1_sp_abs
use mykinds
use abs_mod
implicit real(sp) (Q)
private
public initr1_abs
interface initr1_abs
module procedure initr1_template
end interface initr1_abs
contains
include 'num_init.i90'
end module initr1_sp_abs

module initr1_dp_abs
use mykinds
use abs_mod
implicit real(dp) (Q)
private
public initr1_abs
interface initr1_abs
module procedure initr1_template
end interface initr1_abs
contains
include 'num_init.i90'
end module initr1_dp_abs

module initr1_qp_abs
use mykinds
use abs_mod
implicit real(qp) (Q)
private
public initr1_abs
interface initr1_abs
module procedure initr1_template
end interface initr1_abs
contains
include 'num_init.i90'
end module initr1_qp_abs

module initr1_sp_cmplx
use mykinds
use cmplx_mod
implicit real(sp) (Q)
private
public initr1_cmplx
interface initr1_cmplx
module procedure initr1_template
end interface initr1_cmplx
contains
include 'num_init.i90'
end module initr1_sp_cmplx

module initr1_dp_cmplx
use mykinds
use cmplx_mod
implicit real(dp) (Q)
private
public initr1_cmplx
interface initr1_cmplx
module procedure initr1_template
end interface initr1_cmplx
contains
include 'num_init.i90'
end module initr1_dp_cmplx

module initr1_qp_cmplx
use mykinds
use cmplx_mod
implicit real(qp) (Q)
private
public initr1_cmplx
interface initr1_cmplx
module procedure initr1_template
end interface initr1_cmplx
contains
include 'num_init.i90'
end module initr1_qp_cmplx

module initr1_mod
use initr1_sp_int
use initr1_dp_int
use initr1_qp_int
use initr1_sp_abs
use initr1_dp_abs
use initr1_qp_abs
use initr1_sp_cmplx
use initr1_dp_cmplx
use initr1_qp_cmplx
end module initr1_mod

program initr1_test
use mykinds
use initr1_mod
use name_mod
implicit none
real(sp) xr
real(dp) xd
real(qp) xq

!DEC$ IF DEFINED(.FALSE.)
call xfinit ! Workaround for gfortran bug
!DEC$ ENDIF
xr = int_X
xd = int_X
xq = int_X
call initr1_int(xr)
call initr1_int(xd)
call initr1_int(xq)
xr = abs_X
xd = abs_X
xq = abs_X
call initr1_abs(xr)
call initr1_abs(xd)
call initr1_abs(xq)
xr = cmplx_X
xd = cmplx_X
xq = cmplx_X
call initr1_cmplx(xr)
call initr1_cmplx(xd)
call initr1_cmplx(xq)
end program initr1_test

C:\gfortran\test\num_init>gfortran initr1.f90 finit.s -oinitr1

C:\gfortran\test\num_init>initr1
Test of INT intrinsic
k2 = 4
is_int = 1
kind_if_real = 4
is_cmplx = 0
is_real = 0
kind_if_int = 4
type_label = INTEGER
INT(X) = 1
Test of INT intrinsic
k2 = 4
is_int = 1
kind_if_real = 4
is_cmplx = 0
is_real = 0
kind_if_int = 4
type_label = INTEGER
INT(X) = 1
Test of INT intrinsic
k2 = 4
is_int = 1
kind_if_real = 4
is_cmplx = 0
is_real = 0
kind_if_int = 4
type_label = INTEGER
INT(X) = 1
Test of ABS intrinsic
k2 = 4
is_int = 0
kind_if_real = 4
is_cmplx = 0
is_real = 1
kind_if_int = 4
type_label = REAL
ABS(X) = 1.7951958
Test of ABS intrinsic
k2 = 8
is_int = 0
kind_if_real = 8
is_cmplx = 0
is_real = 1
kind_if_int = 4
type_label = REAL
ABS(X) = 1.7951958020513104
Test of ABS intrinsic
k2 = 10
is_int = 0
kind_if_real = 10
is_cmplx = 0
is_real = 1
kind_if_int = 4
type_label = REAL
ABS(X) = 1.7951958020513104220
Test of CMPLX intrinsic
k2 = 4
is_int = 0
kind_if_real = 4
is_cmplx = 1
is_real = 0
kind_if_int = 4
type_label = COMPLEX
CMPLX(X) = ( 1.7951958 , 0.0000000 )
Test of CMPLX intrinsic
k2 = 4
is_int = 0
kind_if_real = 4
is_cmplx = 1
is_real = 0
kind_if_int = 4
type_label = COMPLEX
CMPLX(X) = ( 1.7951958 , 0.0000000 )
Test of CMPLX intrinsic
k2 = 4
is_int = 0
kind_if_real = 4
is_cmplx = 1
is_real = 0
kind_if_int = 4
type_label = COMPLEX
CMPLX(X) = ( 1.7951958 , 0.0000000 )

The xfinit subroutine in finit.s was only there to work around a bug
in gfortran for Win64 where it sets up the f.p. control word
improperly and is not essential for understanding this example.
Look at how carefully we have to avoid the consequences of an
unanticipated numeric type or integer overflow in num_init.i90.
This example works for ifort too, with the omission of finit.s.

Now this kind of code is near the edge for both gfortran and
ifort. With gfortran, if line 53 of initr1.f90 is changed from

use name_mod, only: fun=> ABS

to

use name_mod, only: fun=> BESSEL_J0

we get:

C:\gfortran\test\num_init>gfortran initr1.f90 finit.s -oinitr1
num_init.i90:8.25:
Included at initr1.f90:115:

transfer(cmplx(0,1,kind_if_real),fun(Qarg1)))**2
1
Error: 'kind' argument of 'cmplx' intrinsic at (1) must be a constant
num_init.i90:9.52:
Included at initr1.f90:115:

integer, parameter :: is_real = 1-is_int-is_cmplx
1
Error: Parameter 'is_cmplx' at (1) has not been declared or is a variable,
which
does not reduce to a constant expression
num_init.i90:11.22:
Included at initr1.f90:115:

integer(kind_if_int), parameter :: iparam = is_int*fun(Qarg1)
1
Error: Constant expression required at (1)

[snip many more error messages]

So gfortran gets all messed up by that BESSEL_J0. But it can
handle BESSEL_J0 in initialization expressions:

C:\gfortran\test\num_init>type testb.f90
program testb
implicit none
real, parameter :: x = 1.7
real, parameter :: j = BESSEL_J0(x)

write(*,*) x, j
end program testb

C:\gfortran\test\num_init>gfortran testb.f90 -otestb

C:\gfortran\test\num_init>testb
1.7000000 0.39798483

Also it can handle weird new functions in a similar context.
If we change the fateful line 53 to:

use name_mod, only: fun=> ERF

It compiles OK with gfortran. But with that change ifort
gets unhappy:

C:\gfortran\test\num_init>ifort initr1.f90


Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.

num_init.i90(4) : Severe: Please report this error along with the
circumstances
in which it occurred in a Software Problem Report
integer, parameter :: is_int = 1-1/(2+0*fun(Qarg1))*2
^
[ Aborting due to internal error. ]
compilation aborted for initr1.f90 (code 1)

Maybe that's due to an old version of ifort: allowing intrinsics
with REAL inputs and outputs is new to f03, ifort may have caught
up by now. It doesn't take a very large program to test for that
event:

C:\gfortran\test\num_init>type teste.f90
program teste
implicit none
real, parameter :: x = 1.7
real, parameter :: j = ERF(x)

write(*,*) x, j
end program teste

C:\gfortran\test\num_init>ifort teste.f90


Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.

teste.f90(4) : Severe: Please report this error along with the circumstances
in
which it occurred in a Software Problem Report
real, parameter :: j = ERF(x)
^
[ Aborting due to internal error. ]
compilation aborted for teste.f90 (code 1)

What would be an example without ICE? :)

James Van Buskirk

unread,
May 2, 2008, 7:56:49 PM5/2/08
to
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:w5KdnUgg3rA0FYbV...@comcast.com...

> Error: 'kind' argument of 'cmplx' intrinsic at (1) must be a constant

I have succeeded in getting down to a much smaller version:

C:\gfortran\test\num_init>type bug3.f90
module funcs
! These fail
intrinsic BESSEL_J0
intrinsic BESSEL_J1
intrinsic BESSEL_Y0
intrinsic BESSEL_Y1
intrinsic BESJ0

! These work
intrinsic ERF
intrinsic ACOSH
end module funcs

program bug3
use funcs, only: fun => BESSEL_J0
implicit none
real, parameter :: Qarg1 = 1.7


integer, parameter :: k2 = kind(fun(Qarg1))
integer, parameter :: is_int = 1-1/(2+0*fun(Qarg1))*2
integer, parameter :: kind_if_real = &
(1-is_int)*k2+is_int*kind(1.0)
integer, parameter :: is_cmplx = -((1-is_int)* &
transfer(cmplx(0,1,kind_if_real),fun(Qarg1)))**2

end program bug3

C:\gfortran\test\num_init>gfortran bug3.f90 -obug3
bug3.f90:23.25:

transfer(cmplx(0,1,kind_if_real),fun(Qarg1)))**2
1
Error: 'kind' argument of 'cmplx' intrinsic at (1) must be a constant

C:\gfortran\test\num_init>gfortran -v


Built by Equation Solution (http://www.Equation.com).
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:

../gcc-4.4-20080425-mingw/configure --host=x86_64-pc-mingw32 --
build=x86_64-unknown-linux-gnu --target=x86_64-pc-mingw32 --prefix=/home/gfortra
n/gcc-home/binary/mingw32/native/x86_64/gcc/4.4-20080425 --with-gmp=/home/gfortr
an/gcc-home/binary/mingw32/native/x86_64/gmp --with-mpfr=/home/gfortran/gcc-home
/binary/mingw32/native/x86_64/mpfr --with-sysroot=/home/gfortran/gcc-home/binary
/mingw32/cross/x86_64/gcc/4.4-20080425 --with-gcc --with-gnu-ld --with-gnu-as

--
disable-shared --disable-nls --disable-tls --enable-languages=c,fortran --enable
-threads=win32 --enable-libgomp --disable-win32-registry
Thread model: win32

gcc version 4.4.0 20080425 (experimental) (GCC)

This version fails if fun is one of the Bessel functions listed
but compiles if, for example, ERF or ACOSH is used. Also it
compiles if the offending statement is commented out, which
indicates that the compiler has lost its way.

Tobias Burnus

unread,
May 3, 2008, 8:38:33 AM5/3/08
to
On May 3, 1:56 am, "James Van Buskirk" <not_va...@comcast.net> wrote:
> > Error: 'kind' argument of 'cmplx' intrinsic at (1) must be a constant
> I have succeeded in getting down to a much smaller version:
>
> C:\gfortran\test\num_init>type bug3.f90
> module funcs
> ! These fail
>    intrinsic BESSEL_J0
>    intrinsic BESSEL_J1
>    intrinsic BESSEL_Y0
>    intrinsic BESSEL_Y1
>    intrinsic BESJ0

The problem is that the front-end of gfortran currently does not
simplify the calls (with constant arguments) to Bessel functions.
(However, doing so is relatively simple.) I think they are later
simplified by the middle end, but that is too late in this case.

I filled a bug report http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36117
(Missed optimization and valid Fortran 2008 is rejected.)

Thanks for reporting this bug*.

Tobias

* If you want to be sure that a bug does not get forgotten, please use
Bugzilla. I think most of the time someone will find the bug report
here, but I'm sure some of the reports are missed.