Passing complex arrays as real ones (fftpack)
Frank Wiedmann
Hello,
When looking at the routines cfftf and cfftb of fftpack from netlib,
I found that in the documentation and in the test program, the array c
is declared complex. In the subroutines, however, there is only
a dimension statement, and so there the array is supposed to be real.
This leads me to the following questions:
1) Why does this work? [I hope it does :-)]
2) What advantages (if any) does this have?
3) Is this conform with the standard?
4) Will this work with f2c?
Thank you for any enlightenment,
Frank
E-Mail: wied...@ihf.uni-stuttgart.de
Tony Willis
>Hello,
>When looking at the routines cfftf and cfftb of fftpack from netlib,
>I found that in the documentation and in the test program, the array c
>is declared complex. In the subroutines, however, there is only
>a dimension statement, and so there the array is supposed to be real.
>This leads me to the following questions:
>1) Why does this work? [I hope it does :-)]
Yep, FORTRAN just passes everything by address, so as long as you are
passing around the right amount of numbers it will work, ie. one
complex number just equals an array of two reals.
>2) What advantages (if any) does this have?
Probably none, if you're concerned about a programming language
having strong type checking.
>3) Is this conform with the standard?
I don't know that there is any standard on this point.
>4) Will this work with f2c?
Nope: f2c will let a few real nasties in FFTPACK through (i.e. like
passing an array of ints to a subroutine, which then starts using them
as reals - or vice versa - can't remember which). When you
try to compile with your C or C++ compiler you will indeed get errors
about passing to the wrong data type. You will have to do some hand
editing of your f2c'd code in order to get it compiled.
Tony
___________
Tony Willis
Internet : twi...@drao.nrc.ca
Snailnet : Dominion Radio Astrophysical Observatory
P.O. Box 248, Penticton, BC, Canada V2A 6K3
BC Tel net: (604) 493-2277 Faxnet : (604) 493-7767
voicemailnet: (604) 490-4343 Localnet : ext 343
robin
>In <2ii9rp$1b...@info2.rus.uni-stuttgart.de> wied...@ihfrs4.ihf.uni-stuttgart.de (Frank Wiedmann) writes:
>When looking at the routines cfftf and cfftb of fftpack from netlib,
>I found that in the documentation and in the test program, the array c
>is declared complex. In the subroutines, however, there is only
>a dimension statement, and so there the array is supposed to be real.
>This leads me to the following questions:
>1) Why does this work? [I hope it does :-)]
>2) What advantages (if any) does this have?
>3) Is this conform with the standard?
>4) Will this work with f2c?
Fortran passes by reference.
Here's an example:
program myprog
real a(10)
call mysub(a(1),5)
end
subroutine mysub(a,n)
complex a(n)
..
.
return
end
Note: A complex number in fortran is essentially stored as 2 real numbers.
Hence, in my main program I can store 10 real numbers. However, in the
subroutine I can only store 5 complex numbers.
---I understand that this is not standard, and know that some systems will
reject this code because the types don't match (real /= complex, by any
stretch of the imagination !!) In any case, it's just bad practice.
Michel Olagnon, Ifremer DITI GO, 98.22.41.44
>
>Hello,
>
>When looking at the routines cfftf and cfftb of fftpack from netlib,
>I found that in the documentation and in the test program, the array c
>is declared complex. In the subroutines, however, there is only
>a dimension statement, and so there the array is supposed to be real.
>
>This leads me to the following questions:
>
>1) Why does this work? [I hope it does :-)]
>2) What advantages (if any) does this have?
Apart from the reasons posted by others about separate use of real and
imaginary parts to avoid duplicating operations, there may well be
some historical reason: the Fortran 77 standard consists of a full
standard, with complex data type, and of a subset, without the
complex data type. Anyway, double complex were not included in the
77 standard.
So, once upon a time (10-15 years ago), library routines using complex
were written with real (easy to convert to double precision) type
arrays having an additional leading dimension of 2. They could thus
run on any system, and required no modification, whether or not the
calling program had complex type. They were usually tested on systems
allowing (double) complex types, so the test programs need not be as
portable.
>3) Is this conform with the standard?
As long as no scalars are implied, and I guess that it is the case,
yes. I even think that backward compatibility may be the reason
why the standard specifies how complex type should be represented
(real and imaginary part stored as two REAL type values)
Michel
---
| Michel OLAGNON email : Michel....@ifremer.fr|
| IFREMER: Institut Francais de Recherches pour l'Exploitation de la Mer|
Shiva Shenoy
>Hello,
>When looking at the routines cfftf and cfftb of fftpack from netlib,
>I found that in the documentation and in the test program, the array c
>is declared complex. In the subroutines, however, there is only
>a dimension statement, and so there the array is supposed to be real.
>This leads me to the following questions:
>1) Why does this work? [I hope it does :-)]
>2) What advantages (if any) does this have?
>3) Is this conform with the standard?
>4) Will this work with f2c?
Fortran passes by reference.
Here's an example:
program myprog
real a(10)
call mysub(a(1),5)
end
subroutine mysub(a,n)
complex a(n)
..
.
return
end
Note: A complex number in fortran is essentially stored as 2 real numbers.
Hence, in my main program I can store 10 real numbers. However, in the
subroutine I can only store 5 complex numbers.
--
----------------------------------------------------------------------
Shiva Shenoy | e-mail: she...@iastate.edu
2021 Black,Dept of AEEM,ISU,Ames,IA 50011 | Office: (515)-294-0092
Alan M. McKenney
bur...@mole.gnu.ai.mit.edu (Craig Burley) writes:
>In article <2ii9rp$1b...@info2.rus.uni-stuttgart.de>
wied...@ihfrs4.ihf.uni-stuttgart.de (Frank Wiedmann) writes:
>
> When looking at the routines cfftf and cfftb of fftpack from netlib,
> I found that in the documentation and in the test program, the array c
> is declared complex. In the subroutines, however, there is only
> a dimension statement, and so there the array is supposed to be real.
>
> This leads me to the following questions:
>
> 3) Is this conform with the standard?
>
>more carefully to be sure.
My reading of the standard says it is *not* conformant to the
standard. Dummy & actual arguments must match in type (as
well as other ways...), and I read a recent posting in this
group which described an implementation (DEC alpha?) on which
mismatching the data-type would cause the program to fail.
I think TOOLPACK catches this, for one.
I once asked whether I might ever run into trouble passing
a REAL array where a COMPLEX array was expected, or vice
versa -- I was warned that the standard would allow an
implementation to pass a COMPLEX argument in a different
way from a REAL one. On the ETA-10 (may it rest in peace)
it might have even been advantageous to pass the real and
imaginary parts as separate arrays. (And, yes, it would
be legal to store them separately, as long as the COMPLEX
array wasn't in COMMON or EQUIVALENCEd to a non-complex
array.)
In fact, you don't have to look far to find implementations
where some kind of type mismatch can cause trouble.
On most UNIX implementations of FORTRAN, CHARACTER arguments
are passed as *two* arguments -- an address, and a length.
If you store a character string in, say, and INTEGER, and
pass it to a subroutine which then tries to get its length
(via LEN), your program could well crash.
A meta-comment
I notice that some follow-ups say things like "FORTRAN passes
arguments by address", or the like. This is not quite accurate.
The FORTRAN standard does *not* specify how *anything* is to
be implemented. It only specifies how a standard-conforming
program should be interpreted.
The fact that most of the machines that posters here are
familiar with *do*, in fact, the same argument-passing
mechanism, that is, passing a type-indepent address
for each argument, is due mostly to the fact that most
computers in use today have pretty much the same architecture,
and *not* to the FORTRAN standard as such. Pick a different
architecture, (for example, a distributed memory parallel
processor), and the trade-offs will change, and passing
an address might be a very *poor* way to pass an argument.
In fact, I would propose adding to the FORTRAN FAQ the
following warning:
When advising people about what is "legal" FORTRAN, it it
is important *not* to assume that what works on your
machine or even on the five or 10 machines you have used,
is "standard FORTRAN".
(I suppose I could come up with a Q&A format version of this,
upon request.)
Alan McKenney E-mail: mcke...@cims.nyu.edu (INTERNET)
Courant Institute,NYU,USA ...!cmcl2!cims.nyu.edu!mckenney (UUCP)
Ob.Disclaimer: Even *I* don't take responsibility for what I post.
NYU *definitely* doesn't, and, if challenged, would probably deny I wrote it.
--
Alan McKenney E-mail: mcke...@cims.nyu.edu (INTERNET)
Courant Institute,NYU,USA ...!cmcl2!cims.nyu.edu!mckenney (UUCP)
Ob.Disclaimer: Even *I* don't take responsibility for what I post.
NYU *definitely* doesn't, and, if challenged, would probably deny I wrote it.
John D McCalpin
Frank Wiedmann <wied...@ihfrs4.ihf.uni-stuttgart.de> wrote:
>
>When looking at the routines cfftf and cfftb of fftpack from netlib,
>I found that in the documentation and in the test program, the array c
>is declared complex. In the subroutines, however, there is only
>a dimension statement, and so there the array is supposed to be real.
>
>This leads me to the following questions:
>
>1) Why does this work? [I hope it does :-)]
This works because almost all Fortran compilers store complex numbers
as if they were two contiguous real numbers. So if you pay attention
to what you are doing, you can emulate the complex arithmetic by working
on the real and imaginary parts separately.
>2) What advantages (if any) does this have?
Many compilers produce sub-optimal code for complex arithmetic.
In addition, there are a number of low-level optimizations that
can take place in the kernel of the FFT algorithm that are only
visible if the complex arithmetic is written out in terms of its
real components.
>3) Is this conform with the standard?
I am not sure if this type of storage equivalence is guaranteed
to work. It does work on almost all modern Fortran systems.
>4) Will this work with f2c?
I have run into considerable trouble with exactly this sort of
coding style and f2c in the past. I gave up rather than trying
to hunt down the details.
One more important note about FFTPACK:
The real to complex routines in FFTPACK use an unusual format for
packing that probably contributes to the confusion experienced by
f2c.
In a real to complex transform, there are N real input elements
and N/2+1 complex output elements. Thus the output appears to have
two more degrees of freedom than the input. In fact, it does not,
since the imaginary part of the first and last complex elements
are identically zero. If one wants to pack the output array into
the same amount of storage as the input array, one has basically
two choices:
(1) put the real part of the last output element into the
imaginary part of the first output element, then throw
away the last output element, or
(2) shift all of the output elements "down" toward the
beginning of the series to overwrite the imaginary part
of the first element.
In the first case, only the first element is modified, while in the
second case, all of the elements are modified. The first case is
much more common in libraries, while the second case is the one used
by FFTPACK. I think that this probably contributes to the confusion
that f2c has with the code.
--
John D. McCalpin mcca...@perelandra.cms.udel.edu
Assistant Professor mcca...@brahms.udel.edu
College of Marine Studies, U. Del. John.M...@mvs.udel.edu
Dave Love
Alan> My reading of the standard says it is *not* conformant to the
Alan> standard.
Hear, hear to that and the rest.
Alan> I think TOOLPACK catches this, for one.
Yes. NAG f90 also rejects the horrible old code that uses this
technique. (I wonder whether that means other compilers based on the
NAG frontend will too.) This is probably the best incentive to sort
out the code, not that I want to contemplate the task I'm faced with
one day.
Another thing I've come across in fft-related code which may be
widespread and is presumably worth worrying about is things like
call fred (a (1), a (2))
where A is updated by fred -- an illegal argument association or
whatever the term is.
Michel Beland
> The MIPS R3000 and R4000 line as well as the Alpha AXP architecture
> specify different register sets for passing integer and floating
> arguments. Of course it depends on the calling conventions followed
> by the particular operating system, language and/or vendor, but there
> exist many current systems in which some scalar arguments are passed
> directly in registers and not by address. An application which
> assumes that all arguments are passed by address is nonportable (and
> perhaps even broken to begin with.)
I do not quite understand what you mean here. If a scalar
argument is not passed by address, how will I make sure that the actual
argument will be modified when I change the dummy argument in the
subroutine ? I guess the answer lies in what you said about registers,
but please explain it.
This might explain what I read in the HPCC User Guide (HPCC is
the High Performance Computing Centre in Calgary) about one of the
compiling options of the Fujitsu Fortran compiler:
The "-Ab" option that selects dummy arguments to be passed by
address (as opposed to by value) produces "safer" code that
may be less prone to failure. (Failures are generally due to
coding practices that violate the ANSI standard.) There is a
performance penalty when the number of arguments is
significant.
When I read the above paragraph, I almost choked. I thought
that passing arguments by value was not standard in Fortran, so I
agree that using this option is safer... (I might add that the
Fujitsu compiler can pass arguments by value only for scalars; arrays
are always passed by address.)
Dik T. Winter
> When I read the above paragraph, I almost choked. I thought
> that passing arguments by value was not standard in Fortran, so I
> agree that using this option is safer... (I might add that the
> Fujitsu compiler can pass arguments by value only for scalars; arrays
> are always passed by address.)
>
On the other hand it requires a lot of cross-analysis to do this on the
callers end. To do it on the callees end is simple, and I know compilers
that do that.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924098
home: bovenover 215, 1025 jn amsterdam, nederland; e-mail: d...@cwi.nl
Steve Lionel
In article <beland.7...@rimbaud.CERCA.UMontreal.CA>,
|>lio...@quark.enet.dec.com (Steve Lionel) writes:
|>
|>> The MIPS R3000 and R4000 line as well as the Alpha AXP architecture
|>> specify different register sets for passing integer and floating
|>> arguments. Of course it depends on the calling conventions followed
|>> by the particular operating system, language and/or vendor, but there
|>> exist many current systems in which some scalar arguments are passed
|>> directly in registers and not by address. An application which
|>> assumes that all arguments are passed by address is nonportable (and
|>> perhaps even broken to begin with.)
|>
|> I do not quite understand what you mean here. If a scalar
|>argument is not passed by address, how will I make sure that the actual
|>argument will be modified when I change the dummy argument in the
|>subroutine ? I guess the answer lies in what you said about registers,
|>but please explain it.
Well, it turns out that I was mistaken as to what happens on AXP systems.
DEC Fortran always passes addresses of arguments (unless you specifically
say %VAL to get immediate value). I misread part of the calling standard
(for DEC OSF/1 AXP); it describes how floating arguments passed by
immediate value are passed in floating registers (for the first six
argument positions) and integer arguments are passed in the integer
registers.
However, one could conceive of a system such as I described. All you'd
need is agreement between the caller and callee as to where the argument
was passed, with a copy-in, copy-out as needed. You can't really do that
with Fortran, though, which has reference semantics.
My apologies for the misinformation, and I hope I didn't cause any needless
consternation.
(Function return values, though, are returned in different register sets,
so you can't return a floating function value and where the caller thinks
it's an integer function, and vice versa.)
--
Steve Lionel lio...@quark.enet.dec.com or
SDT Languages Group lio...@quark.zko.dec.com
Digital Equipment Corporation
110 Spit Brook Road, ZKO2-3/N30
Nashua, NH 03062-2698 "Free advice is worth every cent"
Craig Burley
>>>>> "Alan" == Alan M McKenney <mcke...@gauss.cims.nyu.edu> writes:
Alan> My reading of the standard says it is *not* conformant to the
Alan> standard.
Hear, hear to that and the rest.
Yes, it looks like my "guess" was wrong -- after reading the standard,
it does seem that one cannot do (in standard F77):
REAL A(10)
CALL FOO(A)
END
SUBROUTINE FOO(C)
COMPLEX C(5)
END
Nor can one do:
COMPLEX C(5)
CALL FOO(C)
END
SUBROUTINE FOO(A)
REAL A(10)
END
(The latter case above is the one originally described.)
However, it seems clear the standard _does_ permit:
REAL A(10)
COMPLEX C(5)
EQUIVALENCE (A,C)
CALL FOO(C)
END
SUBROUTINE FOO(C)
COMPLEX C(5)
END
And:
COMPLEX C(5)
REAL A(10)
EQUIVALENCE (A,C)
CALL FOO(A)
END
SUBROUTINE FOO(A)
REAL A(10)
END
This would appear to indicate that any implementation that passes
a COMPLEX array in a manner different than it would pass the
storage-associated REAL array must know whether the COMPLEX
array being passed by the caller is forced, by explicit storage
association (EQUIVALENCE), to associate the storage in the
standard-conforming manner.
That is, someone suggested that a particular compiler out there
might be passing COMPLEX arrays as two separate arrays, one of
the real parts, one of the imaginary parts. Presumably such a
compiler must compile the procedures accepting COMPLEX array arguments
accordingly, but must do so in a way that handles the possibility
that, in one call, the caller is using this separate-array technique,
but in another call, the (perhaps same) caller has associated COMPLEX
and REAL arrays. (There are ways of doing this, of course, such as
requiring all callers passing COMPLEX arrays that, via EQUIVALENCE or
COMMON, have been storage-associated in the F77-required way, to do
copying before and after the call to a temporary pair of arrays that
hold the real and imaginary parts, as expected by the procedure.
The F77 standard has other restrictions that makes this sort of thing
fairly straightforward.)
In summary, yes, it _is_ okay to write a library of F77 procedures
that declare arrays of complex numbers as REAL arrays, and operate
on those arrays accordingly, as long as the callers _also_ declare
the passed array arguments as REAL. And they can almost always do
this fairly easily via EQUIVALENCE (the hard part is when the caller
is itself a procedure expecting a COMPLEX array as an argument,
in which case EQUIVALENCE cannot be used in a standard F77 way to
associate the dummy REAL array argument with a (dummy) COMPLEX
array argument).
Of course, I also agree with everyone else's statements that while
_many_ F77 compilers use "pass by reference", that is purely an
implementation issue and not a requirement of the standard, and thus
should not be depended upon by code intended to be portable.
--
James Craig Burley, Software Craftsperson bur...@gnu.ai.mit.edu
Member of the League for Programming Freedom (LPF) l...@uunet.uu.net