Fortran - FFTW and ifort

Author

FFTW and ifort

Jacek

2004-10-08, 3:57 pm

Could you help me with solving the following problem: I'm trying to
call FFTW v.3.0.1 from ifort Intel compiler (v.8.1). Unfortunately I
get the following error messsage:

Thread received signal SEGV
stopped at [void m1_32_0(const R*, const R*, R*, R*, stride,
stride):62 0x0807048e]
62 T2 = ri[WS(is, 16)];
(idb) where
>0 0x0807048e in m1_32_0(ri=0xbffffba0, ii=0xbffffba8, ro=0xbffffbd0,
io=0xbffffbd8, is=0x8184500, os=0x8184588) "m1_32.c":62
#1 0x08072c9f in m1_32(ri=0xbffffba0, ii=0xbffffba8, ro=0xbffffbd0,
io=0xbffffbd8, is=0x8184500, os=0x8184588, v=8, ivs=2, ovs=64)
"m1_32.c":673
#2 0x08111ed5 in apply(ego_=0x81844a8, ri=0xbffffba0, ii=0xbffffba8,
ro=0xbffffbd0, io=0xbffffbd8) "direct.c":47
#3 0x08110e55 in apply(ego_=0x8184610, ri=0xbffffba0, ii=0xbffffba8,
ro=0xbffffbd0, io=0xbffffbd8) "ct-dit.c":34
#4 0x08056cbc in apply(ego_=0x81848c8, ri=0xbfffdba0, ii=0xbfffdba8,
ro=0xbfffdbd0, io=0xbfffdbd8) "vrank-geq1.c":63
#5 0x0805628e in apply(ego_=0x8182be8, ri=0xbfffdba0, ii=0xbfffdba8,
ro=0xbfffdbd0, io=0xbfffdbd8) "rank-geq2.c":49
#6 0x08056c2c in fftw_dft_solve(ego_=0x8182be8, p_=0x8180978)
"solve.c":30
#7 0x0804cd32 in dfftw_execute_(p=0xbfffdc00) "f77funcs.h":27
#8 0x0804ae0f in CONVOLUTION::get_potential(distro=(...), nn=(...),
size=<no value> ) "m_pot_2dim.f90":98
#9 0xb7741866
(idb)

Subroutine which calls FFTW is named get_potential(). What am I doing
wrong?
The weird thing is that exactly the same code runs well when pgf90
compiler is used.

Thanks a lot.
Jacek

Kamaraju Kusumanchi

2004-10-09, 10:28 pm

> Subroutine which calls FFTW is named get_potential(). What am I doing
> wrong?
> The weird thing is that exactly the same code runs well when pgf90
> compiler is used.
>

It might be a wise idea to post the minimalized version of the code
which mimicks this behaviour. That way others will have a chance to look
at the code and determine if it really is a fftw bug.

thanks
raju

Jacek

2004-10-11, 8:56 am

Kamaraju Kusumanchi wrote in message
> It might be a wise idea to post the minimalized version of the code
> which mimicks this behaviour. That way others will have a chance to look
> at the code and determine if it really is a fftw bug.
>
> thanks
> raju

Here is the minimalized version of my program (main and module). It
raises the same error message as the whole one (error message comes at
the very end).

Thanks a lot,
Jacek

MAIN:
------------------------------------
program po
use con
implicit none
real(DP), allocatable :: den(:,:)
integer(8) :: nn(1:2)
integer(8) :: i,j
integer(8) :: allocate_status
nn(1:2) = 128
allocate(den(1:2*nn(1),1:2*nn(2)), STAT = allocate_status)
print*, allocate_status
do i=1, 2*nn(1)
do j=1, 2*nn(2)
den(i,j) = 0.0
if (sqrt((i-64.0)**2 + (j-64.0)**2) < 4) den(i,j) = 1.0
end do
end do
! code seems to fail at the following subroutine:
call get(den, nn)
end program po
--------------------------------

MODULE:
------------------------------------
module con
implicit none
integer, parameter :: DP = kind(1.0d0)
integer(8), parameter :: FFT_ESTIMATE = 64
contains
subroutine get(distro, nn)
implicit none
interface
subroutine dfftw_plan_dft_2d(plan,nx,ny,in,out,dir,
FFT_ESTIMATE)
integer(8), intent(in) :: plan
integer(8), intent(in) :: nx, ny
complex(kind(1d0)), intent(inout) :: in(:,:), out(:,:)
integer(8), intent(in) :: dir, FFT_ESTIMATE
end subroutine dfftw_plan_dft_2d
subroutine dfftw_execute(plan)
integer(8), intent(in) :: plan
end subroutine dfftw_execute
subroutine dfftw_destroy_plan(plan)
integer(8), intent(in) :: plan
end subroutine dfftw_destroy_plan
end interface
integer(8), intent(in) :: nn(1:2)
real(DP), intent(inout) :: distro(1:2*nn(1),1:2*nn(2))
complex(DP) :: kdistro(1:2*nn(1),1:2*nn(2))
integer(8) :: plan
integer(8) :: direction
integer(8) :: nx, ny
nx = 2*nn(1)
ny = 2*nn(2)
kdistro(1:nx,1:ny) = cmplx(distro(1:nx,1:ny))
direction = -1
call dfftw_plan_dft_2d(plan, nx, ny, kdistro, kdistro,direction,
FFT_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
end subroutine get
end module con
--------------------------------------------------

ERROR MESSAGE:
--------------------------------------------------
(idb) run
0
Thread received signal SEGV
stopped at [void m1_32_0(const R*, const R*, R*, R*, stride,
stride):73 0x0806e2d4]
73 T5 = ri[WS(is, 24)];
(idb) where
>0 0x0806e2d4 in m1_32_0(ri=0xbffff4a0, ii=0xbffff4a8, ro=0xbffff4d0,
io=0xbffff4d8, is=0x8179af8, os=0x8179b80) "m1_32.c":73
#1 0x08070a17 in m1_32(ri=0xbffff4a0, ii=0xbffff4a8, ro=0xbffff4d0,
io=0xbffff4d8, is=0x8179af8, os=0x8179b80, v=8, ivs=2, ovs=64)
"m1_32.c":673
#2 0x0810fc4d in apply(ego_=0x8179aa0, ri=0xbffff4a0, ii=0xbffff4a8,
ro=0xbffff4d0, io=0xbffff4d8) "direct.c":47
#3 0x0810ebcd in apply(ego_=0x8179c08, ri=0xbffff4a0, ii=0xbffff4a8,
ro=0xbffff4d0, io=0xbffff4d8) "ct-dit.c":34
#4 0x08054a34 in apply(ego_=0x8179ec0, ri=0xbfffe4a0, ii=0xbfffe4a8,
ro=0xbfffe4d0, io=0xbfffe4d8) "vrank-geq1.c":63
#5 0x08054006 in apply(ego_=0x81781e0, ri=0xbfffe4a0, ii=0xbfffe4a8,
ro=0xbfffe4d0, io=0xbfffe4d8) "rank-geq2.c":49
#6 0x080549a4 in fftw_dft_solve(ego_=0x81781e0, p_=0x8175f70)
"solve.c":30
#7 0x0804aaaa in dfftw_execute_(p=0xbfffe500) "f77funcs.h":27
#8 0x0804a427 in CON::get(distro=(...), nn=(...)) "m_tmp.f90":36
#9 0xc07f5967
(idb)

Herman D. Knoble

2004-10-11, 4:00 pm

Dick Hendrickson

2004-10-11, 4:00 pm

Two thoughts.

First, put a "print *, 'we are here'" before and after the
calls to the three dfftw_* routines, then you'll know
which subroutine is the culprit.

Second, I'd worry about passing kdistro to both the
"in" and "out" arguments of dfftw_plan_dft_2d. The
rules of Fortran don't allow aliased arguments like
these to be changed in the subroutine. I'd guess
that the "out" argument gets changed. That's just
against the rules. I think you need a separate
array for the out argument.

HTH

Dick Hendrickson

Jacek wrote:
> Kamaraju Kusumanchi wrote in message
>
>
>
>
> Here is the minimalized version of my program (main and module). It
> raises the same error message as the whole one (error message comes at
> the very end).
>
> Thanks a lot,
> Jacek
>
>
> MAIN:
> ------------------------------------
> program po
> use con
> implicit none
> real(DP), allocatable :: den(:,:)
> integer(8) :: nn(1:2)
> integer(8) :: i,j
> integer(8) :: allocate_status
> nn(1:2) = 128
> allocate(den(1:2*nn(1),1:2*nn(2)), STAT = allocate_status)
> print*, allocate_status
> do i=1, 2*nn(1)
> do j=1, 2*nn(2)
> den(i,j) = 0.0
> if (sqrt((i-64.0)**2 + (j-64.0)**2) < 4) den(i,j) = 1.0
> end do
> end do
> ! code seems to fail at the following subroutine:
> call get(den, nn)
> end program po
> --------------------------------
>
> MODULE:
> ------------------------------------
> module con
> implicit none
> integer, parameter :: DP = kind(1.0d0)
> integer(8), parameter :: FFT_ESTIMATE = 64
> contains
> subroutine get(distro, nn)
> implicit none
> interface
> subroutine dfftw_plan_dft_2d(plan,nx,ny,in,out,dir,
FFT_ESTIMATE)
> integer(8), intent(in) :: plan
> integer(8), intent(in) :: nx, ny
> complex(kind(1d0)), intent(inout) :: in(:,:), out(:,:)
> integer(8), intent(in) :: dir, FFT_ESTIMATE
> end subroutine dfftw_plan_dft_2d
> subroutine dfftw_execute(plan)
> integer(8), intent(in) :: plan
> end subroutine dfftw_execute
> subroutine dfftw_destroy_plan(plan)
> integer(8), intent(in) :: plan
> end subroutine dfftw_destroy_plan
> end interface
> integer(8), intent(in) :: nn(1:2)
> real(DP), intent(inout) :: distro(1:2*nn(1),1:2*nn(2))
> complex(DP) :: kdistro(1:2*nn(1),1:2*nn(2))
> integer(8) :: plan
> integer(8) :: direction
> integer(8) :: nx, ny
> nx = 2*nn(1)
> ny = 2*nn(2)
> kdistro(1:nx,1:ny) = cmplx(distro(1:nx,1:ny))
> direction = -1
> call dfftw_plan_dft_2d(plan, nx, ny, kdistro, kdistro,direction,
> FFT_ESTIMATE)
> call dfftw_execute(plan)
> call dfftw_destroy_plan(plan)
> end subroutine get
> end module con
> --------------------------------------------------
>
> ERROR MESSAGE:
> --------------------------------------------------
> (idb) run
> 0
> Thread received signal SEGV
> stopped at [void m1_32_0(const R*, const R*, R*, R*, stride,
> stride):73 0x0806e2d4]
> 73 T5 = ri[WS(is, 24)];
> (idb) where
>
>
> io=0xbffff4d8, is=0x8179af8, os=0x8179b80) "m1_32.c":73
> #1 0x08070a17 in m1_32(ri=0xbffff4a0, ii=0xbffff4a8, ro=0xbffff4d0,
> io=0xbffff4d8, is=0x8179af8, os=0x8179b80, v=8, ivs=2, ovs=64)
> "m1_32.c":673
> #2 0x0810fc4d in apply(ego_=0x8179aa0, ri=0xbffff4a0, ii=0xbffff4a8,
> ro=0xbffff4d0, io=0xbffff4d8) "direct.c":47
> #3 0x0810ebcd in apply(ego_=0x8179c08, ri=0xbffff4a0, ii=0xbffff4a8,
> ro=0xbffff4d0, io=0xbffff4d8) "ct-dit.c":34
> #4 0x08054a34 in apply(ego_=0x8179ec0, ri=0xbfffe4a0, ii=0xbfffe4a8,
> ro=0xbfffe4d0, io=0xbfffe4d8) "vrank-geq1.c":63
> #5 0x08054006 in apply(ego_=0x81781e0, ri=0xbfffe4a0, ii=0xbfffe4a8,
> ro=0xbfffe4d0, io=0xbfffe4d8) "rank-geq2.c":49
> #6 0x080549a4 in fftw_dft_solve(ego_=0x81781e0, p_=0x8175f70)
> "solve.c":30
> #7 0x0804aaaa in dfftw_execute_(p=0xbfffe500) "f77funcs.h":27
> #8 0x0804a427 in CON::get(distro=(...), nn=(...)) "m_tmp.f90":36
> #9 0xc07f5967
> (idb)

Kamaraju Kusumanchi

2004-10-11, 4:00 pm

Jacek wrote:
> Kamaraju Kusumanchi wrote in message
>
>
> Here is the minimalized version of my program (main and module). It
> raises the same error message as the whole one (error message comes at
> the very end).
>
> Thanks a lot,
> Jacek
>

I see couple of points which can be improved upon in the code. I have
not tested the new code though...

>
> MAIN:
> ------------------------------------
> program po
> use con
> implicit none
> real(DP), allocatable :: den(:,:)
> integer(8) :: nn(1:2)
> integer(8) :: i,j
> integer(8) :: allocate_status
> nn(1:2) = 128
> allocate(den(1:2*nn(1),1:2*nn(2)), STAT = allocate_status)
> print*, allocate_status
> do i=1, 2*nn(1)
> do j=1, 2*nn(2)
> den(i,j) = 0.0
> if (sqrt((i-64.0)**2 + (j-64.0)**2) < 4) den(i,j) = 1.0
> end do
> end do
> ! code seems to fail at the following subroutine:
> call get(den, nn)
> end program po
> --------------------------------
>
> MODULE:
> ------------------------------------
> module con
> implicit none
> integer, parameter :: DP = kind(1.0d0)
> integer(8), parameter :: FFT_ESTIMATE = 64

Why are you initializing this flag? See pg-25 of fftw manual where it
says to use any one of FFTW_ESTIMATE, FFTW_MEASURE, FFTW_PATIENT,
FFTW_EXHAUSTIVE, FFTW_EXHAUSTIVE

> contains
> subroutine get(distro, nn)
> implicit none
> interface
> subroutine dfftw_plan_dft_2d(plan,nx,ny,in,out,dir,
FFT_ESTIMATE)
> integer(8), intent(in) :: plan
> integer(8), intent(in) :: nx, ny
> complex(kind(1d0)), intent(inout) :: in(:,:), out(:,:)
> integer(8), intent(in) :: dir, FFT_ESTIMATE
> end subroutine dfftw_plan_dft_2d
> subroutine dfftw_execute(plan)
> integer(8), intent(in) :: plan
> end subroutine dfftw_execute
> subroutine dfftw_destroy_plan(plan)
> integer(8), intent(in) :: plan
> end subroutine dfftw_destroy_plan
> end interface
> integer(8), intent(in) :: nn(1:2)
> real(DP), intent(inout) :: distro(1:2*nn(1),1:2*nn(2))
> complex(DP) :: kdistro(1:2*nn(1),1:2*nn(2))
> integer(8) :: plan
> integer(8) :: direction
> integer(8) :: nx, ny
> nx = 2*nn(1)
> ny = 2*nn(2)
> kdistro(1:nx,1:ny) = cmplx(distro(1:nx,1:ny))
> direction = -1

Although this is correct, I prefer to use FFTW_FORWARD and FFTW_BACKWARD
instead -1 or +1 as it makes the code more readable. But that is a
personal preference anyway.

> call dfftw_plan_dft_2d(plan, nx, ny, kdistro, kdistro,direction,
> FFT_ESTIMATE)
> call dfftw_execute(plan)

Just a caveat. Creating plan destroys the original data in the
variables. So fftw might not be doing what you intend to do. I
frequently run into this problem and is very tough to realize that it is
the cause of the problem.

> call dfftw_destroy_plan(plan)
> end subroutine get
> end module con

Please let us know if you make any progress.

hth
raju

glen herrmannsfeldt

2004-10-11, 4:00 pm

Dick Hendrickson wrote:
> Two thoughts.
>
> First, put a "print *, 'we are here'" before and after the
> calls to the three dfftw_* routines, then you'll know
> which subroutine is the culprit.

From the traceback it seems to be dfftw_execute

> Second, I'd worry about passing kdistro to both the
> "in" and "out" arguments of dfftw_plan_dft_2d. The
> rules of Fortran don't allow aliased arguments like
> these to be changed in the subroutine. I'd guess
> that the "out" argument gets changed. That's just
> against the rules. I think you need a separate
> array for the out argument.

That would be true if it was a Fortran subroutine, but
it seems to be a C routine. Also interesting is that
plan is an intent(in) argument to all, and isn't set before
the first one.

It looks like it has some memory of the arguments passed
to the first one, and actually does the FFT on the
dfftw_execute call. The plan variable is used to keep track
of which one it is doing. I would expect the first call
to set plan, possibly as a C pointer, and the rest to
use it, but it shouldn't be intent(in) in that case.

Dick Hendrickson

2004-10-11, 4:00 pm

glen herrmannsfeldt wrote:
> Dick Hendrickson wrote:
>
>
>
> From the traceback it seems to be dfftw_execute
>
>
>
> That would be true if it was a Fortran subroutine, but
> it seems to be a C routine.
I don't think that matters. I think the rule is that
aliased arguments can't be changed period, and it doesn't
matter how they are changed. That allows the
compiler to do a copy IN/OUT, for example, for arrays.
That might be what's happening here, since the dimensions
and the subscripts aren't completely simple. My rule of
thumb is to usually fix the bugs I know about and hope that
fixes everything.

Dick Hendrickson

Also interesting is that
> plan is an intent(in) argument to all, and isn't set before
> the first one.
>
> It looks like it has some memory of the arguments passed
> to the first one, and actually does the FFT on the
> dfftw_execute call. The plan variable is used to keep track
> of which one it is doing. I would expect the first call
> to set plan, possibly as a C pointer, and the rest to
> use it, but it shouldn't be intent(in) in that case.
>

glen herrmannsfeldt

2004-10-11, 4:00 pm

Dick Hendrickson wrote:

> glen herrmannsfeldt wrote:
>

(snip)
[color=darkred]
[color=darkred]
[color=darkred]
> I don't think that matters. I think the rule is that
> aliased arguments can't be changed period, and it doesn't
> matter how they are changed. That allows the
> compiler to do a copy IN/OUT, for example, for arrays.

In the case of mixing Fortran and C there are many things
that can happen if you don't really know what both are doing.

C doesn't allow copy IN/OUT. Note, though, that in this
case all the arrays are passed to one routine but the
transform is actually done by the second routine. I don't
know if that is legal at all in Fortran, but I believe it
is in C. A C subroutine can keep a pointer to data passed
in a previous call, (or to a call to another routine) and
use it later. My guess is that plan is a structure pointer
returned by the first call and used by the later calls.

> That might be what's happening here, since the dimensions
> and the subscripts aren't completely simple.

That sounds suspicious to me. If the called program was
assuming a simple array very strange things could happen.

> My rule of thumb is to usually fix the bugs I know
> about and hope that fixes everything.

I like this rule, too.

-- glen

Jacek

2004-10-11, 8:56 pm

>
> Please let us know if you make any progress.
>

Unfortunately it doesn't work yet. I got rid of intent(in) and
intent(out) in specification of a FFW routines interfaces (as someone
suggested) but it didn't help.

The code seems to blow up after:

call dfftw_plan_dft_2d(plan, nx, ny, kdistro, tempo,direction, FFT_ESTIMATE

The question which comes to me is why this code works fine when Portland
Fortran compiler is used instead?

Thanks,
Jacek

glen herrmannsfeldt

2004-10-11, 8:56 pm

Jacek wrote:

> Unfortunately it doesn't work yet. I got rid of intent(in) and
> intent(out) in specification of a FFW routines interfaces (as someone
> suggested) but it didn't help.

> The code seems to blow up after:

The traceback seems to say dfftw_execute.

> call dfftw_plan_dft_2d(plan, nx, ny, kdistro, tempo,direction, FFT_ESTIMATE
>The question which comes to me is why this code works fine
> when Portland Fortran compiler is used instead?

Maybe different calling sequence for the different compilers?

http://www.tp1.rub.de/lehre/compphy...n-Examples.html

http://www.compsci.wm.edu/SciClone/.../FFTW/fftw3.pdf

It seems that it expects arrays in contiguous memory, so it might
not take allocatable arrays, or maybe only not multidimensional arrays.

If you copy all the data to a one dimensional array in column major
order (it seems that the dfftw routines reverse the subscript order
when calling C).

plan is a C pointer, and INTEGER*8 is suggested. NX and NY should be
ordinary integer variables. On a little endian machine it might not
matter, but it is best to match the type.

-- glen

Ian Bush

2004-10-12, 3:56 am

Hi Jacek,

Jacek wrote:

> Could you help me with solving the following problem: I'm trying to
> call FFTW v.3.0.1 from ifort Intel compiler (v.8.1). Unfortunately I
> get the following error messsage:
>

SNIP error message.

> Subroutine which calls FFTW is named get_potential(). What am I doing
> wrong?
> The weird thing is that exactly the same code runs well when pgf90
> compiler is used.
>

What was FFTW compiled with ? You might have a problem with inconsistent
argument passing conventions,

Ian

Jacek

2004-10-12, 3:58 pm

Hi Ian,

> What was FFTW compiled with ? You might have a problem with inconsistent
> argument passing conventions,

FFTW was built in the following way:

../configure CC="gcc-3.4" CFLAGS="-g" F77="g77-3.4" FFLAGS="-g"
--prefix=/home/guzik/bin/FFTW
make
make install

and I use it as follows:

ifort -g m.f90 j.f90 -L/home/guzik/bin/FFTW/lib -lfftw3

Do you know what I did wrong?

Thanks,
Jacek

Richard E Maine

2004-10-12, 3:58 pm

glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

> Dick Hendrickson wrote:

>
> In the case of mixing Fortran and C there are many things
> that can happen if you don't really know what both are doing.

I certainly can't disagree with Glen's statement per se, but if it is
intended to be a counter to Dick's point, then I think I'll have to
argue for Dick's position (as I understand it) here. I agree with
Glen's statement as a general caveat that mixing Fortran and C
provides ways to screw up if you don't understand what is
happening. But I claim that many such screwups do constitute
violations of the Fortran standard. That might not be much
consolation, since it doesn't mean that the compiler will necessarily
tell you about them.

The rules of the language still apply within the context of
that language. For a simple example (it is easier to do a trivial
one than to try to illustrate with the particular code in question):
If we have something like

subroutine sub(x)
real, inten(in) :: x
call some_c_thing(x)
...
end

Then x is not allowed to be altered by the C code. The Fortran
compiler is free to assume this (and to give garbage results if
the assumption is incorrect).

Yes, it is possible to write C code that violates this. And it is
possible (likely even) that such a things would be undetected by
the Fortran compiler. But it would still constitute an invalid
program.

Oh. And another widely-cited example involves C code that saves
addresses of its arguments. Something like

subroutine sub2
real :: x
x = 1.
call some_c_thing(x)
x = 42.
call some_other_c_thing
write (*,*) x
end

This had better print 42 for x. I could certainly write C code that
would make it print something else (save the address of x in the
first function; write something to that address in the second). But
that would violate the Fortran standard. Yes, even though you did
the "damage" in C, it still violates the Fortran standard when
you put it all together.

The general rule is that the requirements of the standard hold,
regardless of the language used to implement procedures. The fact that
you can do something in C doesn't mean that you may do it in a
particular context (and still have a standard-conforming program).
Note also that standard-conformance is for a program as a whole.
It is trivial to take several subroutines, each standard-conforming
individually, that do not make a standard-conforming program when
you put them together; you don't need to add C into the mix to do
this.

I make no judgement about whether this might be related to the OP's
problem. I haven't studied it well enough to make such a judgement,
and I'm not sure whether I'd have enough information even if I did
study it more. I am addressing just the general question of
standard-conformance of mixed-language programs. I'll note that this
kind of question did get quite substantial discussion when the C
interop features of F2003 were being developed. Of course, that
doesn't mean that I have correctly recalled and applied the
principles, or explained them clearly. The matter is certainly not
trivial.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

glen herrmannsfeldt

2004-10-12, 3:58 pm

Richard E Maine wrote:
> glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

(snip related to copy in/out)

(even bigger snip)
[color=darkred]
> Oh. And another widely-cited example involves C code that saves
> addresses of its arguments. Something like

> subroutine sub2
> real :: x
> x = 1.
> call some_c_thing(x)
> x = 42.
> call some_other_c_thing
> write (*,*) x
> end

> This had better print 42 for x. I could certainly write C code that
> would make it print something else (save the address of x in the
> first function; write something to that address in the second). But
> that would violate the Fortran standard. Yes, even though you did
> the "damage" in C, it still violates the Fortran standard when
> you put it all together.

That seems to be one of the things the OP's program does.

The first call returns a pointer to a C structure (with the suggestion
to use INTEGER*8), and later calls operate on the structure containing
pointers to arrays passed in the first call. The C documentation,
and also the Fortran callable intermediate program, allow the same array
to be specified for an in-place transforms.

It seems that it is failing in the dfftw_execute call that actually
does the transform. Not disagreeing that Fortran doesn't allow
modifying arrays passed as two arguments, I believe the problem
is related to using an ALLOCATABLE array. Though the document
is dated 2003, it never mentions the use of such arrays.

-- glen

Steven G. Johnson

2004-10-12, 8:56 pm

Kamaraju Kusumanchi <kk288@cornell.edu> wrote...
>
> Just a caveat. Creating plan destroys the original data in the
> variables. So fftw might not be doing what you intend to do.

Not in estimate mode, so that's not a problem here. (It is true that,
in FFTW_MEASURE or FFTW_PATIENT mode, the array gets initialized to
all zeros; see also the FFTW FAQ.)

Steven G. Johnson

2004-10-12, 8:56 pm

It seems pretty clear to me that there is some mismatch in
argument-passing conventions going on, although I'm not sure what
specifically is the culprit. (Sigh...this will be much easier once C
interoperability is standardized in Fortran.)

uoguzik@cyf-kr.edu.pl (Jacek) wrote...
> subroutine dfftw_plan_dft_2d(plan,nx,ny,in,out,dir,
FFT_ESTIMATE)
> integer(8), intent(in) :: plan
> integer(8), intent(in) :: nx, ny
> complex(kind(1d0)), intent(inout) :: in(:,:), out(:,:)
> integer(8), intent(in) :: dir, FFT_ESTIMATE

One thing that jumps out at me is that nx and ny, as well as dir and
flags (the last argument) are all of type "int" in C, which on most
systems is a 4-byte integer ...whereas you are passing integer(8),
which I think denotes an 8-byte integer? (The fact that you are on a
little-endian system probably is what has saved you here.)

A second thing is that also that the plan should be intent(out) since
the whole point of dfftw_plan_dft_2d is to initialize the plan (a
pointer). My understanding is that, with intent(in), the compiler is
allowed to not copy the value of the plan back upon return, in which
case you are passing a garbage pointer to dfftw_execute.

Another worry is: does intent(in) allow the compiler to pass arguments
by value? The C code in FFTW expects all of the Fortran-interface
parameters to be passed by reference (as for Fortran 77). To force
this, maybe you need to declare everything as intent(inout)?
Similarly for dfftw_execute and dfftw_destroy_plan.

The Fortran-callable interface in FFTW is designed to be callable from
Fortran 77, so you should be careful with non-F77 types. In
particular, all arrays should be contiguous F77-style
column-major-order blocks. I don't know enough about Fortran 90 to
say whether that is the case for the arrays as you've declared them
here; maybe someone else can chip in on that point.

Cordially,
Steven G. Johnson

Richard E Maine

2004-10-12, 8:56 pm

stevenj@alum.mit.edu (Steven G. Johnson) writes:

[what look to me like some pretty cogent observations]

> One thing that jumps out at me is that nx and ny, as well as dir and
> flags (the last argument) are all of type "int" in C, which on most
> systems is a 4-byte integer ...whereas you are passing integer(8),
> which I think denotes an 8-byte integer? (The fact that you are on a
> little-endian system probably is what has saved you here.)

This kind of thing can easily happen to work on some compilers, but
fail miserably on others. Definitely nonportable, and definitely
to be avoided.

> A second thing is that also that the plan should be intent(out) since
> the whole point of dfftw_plan_dft_2d is to initialize the plan (a
> pointer). My understanding is that, with intent(in), the compiler is
> allowed to not copy the value of the plan back upon return, in which
> case you are passing a garbage pointer to dfftw_execute.

Ooo. Nasty. I hadn't studied what was going on well enough to notice
that one, but yes, that can be a nasty bug. I hit it once in some code
of mine not too many years ago. The code had been working fine for
many years and then started getting segmentation faults with a new
compiler release. Same vendor, just a new release. Seems that the
new release added something like that as an optimization. My code
never should have specified intent(out) for the argument in question;
it should have been intent(inout), but that error had not been caught
before.

One of the things I consider nasty about that bug is that it can
hit you when the compiler gets "better" - as in adds optimizations
that improve performance and are completely "safe" for correct
code.

> Another worry is: does intent(in) allow the compiler to pass arguments
> by value?

The standard doesn't address implementation. Passing by value is
an implementation choice. So yes, it could happen (though I think
it is rare and not likely to happen in practice for severral reasons).

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

James Van Buskirk

2004-10-12, 8:56 pm

"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
news:27cfb406.0410121214.3d3ab3e7@posting.google.com...

> uoguzik@cyf-kr.edu.pl (Jacek) wrote...
dfftw_plan_dft_2d(plan,nx,ny,in,out,dir,
FFT_ESTIMATE)[color=darkred]

Hey, Steven. You should have quoted the INTERFACE statement as well
because that's what says the given code is doomed, doomed, doomed.
IN & OUT are declared in Fortran as assumed-shape dummies in (what I
would imagine is) you C function returning void.
[color=darkred]
> The Fortran-callable interface in FFTW is designed to be callable from
> Fortran 77, so you should be careful with non-F77 types. In
> particular, all arrays should be contiguous F77-style
> column-major-order blocks. I don't know enough about Fortran 90 to
> say whether that is the case for the arrays as you've declared them
> here; maybe someone else can chip in on that point.

void dfftw_plan_dft_2d() will ordinarily see addresses of descriptors
for assumed-shape arguments corresponding to the 4th & 5th arguments
(in & out). If this function is in your library I would consider it
quite a programming tour de force if you were able to make the above
invocation work with both correct (F77-style) and incorrect (as above)
interfaces for even a couple of compilers.

P.S. have you had a chance to take a look at my convolution stuff?

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

B52B

2004-10-13, 3:56 am

Jacek wrote:
>
> Hi Ian,
>
>
>
> FFTW was built in the following way:
>
> ./configure CC="gcc-3.4" CFLAGS="-g" F77="g77-3.4" FFLAGS="-g"
> --prefix=/home/guzik/bin/FFTW

As you know the Fortran - C -Fortran interoperability is non standard.
Since you use gcc and ifort it may lead to additional problems.
Why don't you use ifc (Intel C compiler free for non commercial use).

Pay attention for C-Fortran basic types correspondence.
Integer in Fortran means long integer in C (on Linux IA32), look at
the compiler documentation.

Furthermore FFTW has f77 api so you should omit intent's in your code.

Have you looked at FFTW Fortran examples?
I've seen FFTW distributions with examples packed.

And last but not least there are many implementations of FFT in
Fortran77 and Fortran90/95. They are comparable to FFTW in speed and
easier to use (sometimes faster - personal observation).

Regards,
B52B

Take

> make
> make install
>
> and I use it as follows:
>
> ifort -g m.f90 j.f90 -L/home/guzik/bin/FFTW/lib -lfftw3
>
> Do you know what I did wrong?
>
> Thanks,
> Jacek
>

Ian Bush

2004-10-13, 9:09 am

Hi Jacek,

Jacek wrote:

>
> Hi Ian,
>
>
> FFTW was built in the following way:
>
> ./configure CC="gcc-3.4" CFLAGS="-g" F77="g77-3.4" FFLAGS="-g"
> --prefix=/home/guzik/bin/FFTW
> make
> make install
>
> and I use it as follows:
>
> ifort -g m.f90 j.f90 -L/home/guzik/bin/FFTW/lib -lfftw3
>
> Do you know what I did wrong?
>

So the library was built with the gnu compiler suite while you
used the Intel compiler for your program. Do you know if these
are compatible ? It's a while since I did this sort of thing
but somewhere in the back of mind I've got a feeling that there
can be problems with this pair, but I could well be wrong,

Ian

Janne Blomqvist

2004-10-13, 9:09 am

On 2004-10-13, B52B <B52B@pl.pl.pl> wrote:
> As you know the Fortran - C -Fortran interoperability is non standard.
> Since you use gcc and ifort it may lead to additional problems.

In my experience, ifort understands gcc/g77 calling conventions. I
doubt this is the source of the problems.

> Pay attention for C-Fortran basic types correspondence.
> Integer in Fortran means long integer in C (on Linux IA32), look at
> the compiler documentation.

Since on Linux IA32 the C types long and int are the same thing, the
Fortran INTEGER is also the C integer.

--
Janne Blomqvist

Richard E Maine

2004-10-13, 3:59 pm

B52B <B52B@pl.pl.pl> writes:

> Since you use gcc and ifort it may lead to additional problems.
> Why don't you use ifc (Intel C compiler free for non commercial use).

As Janne said, I wouldn't expect this to be a problem. I've used
ifc with gcc and saw no issues. I did so only briefly and not with
the most recent release, but I'd look elsewhere.

> Furthermore FFTW has f77 api so you should omit intent's in your code.

Intent does not matter. It can't matter because dummy arguments with
intent specified do not trigger the requirement for explicit interfaces.
Thus, the compiler doesn't necessarily know the intents when compiling
the calling program.

Assumed shape on the other hand...that matters a *LOT*. See James
Buskirk's post. If the subroutine really has an f77 interface, then
it for sure doesn't use assumed shape. So it looks to me like the
interface specified probably disagrees with the actual subroutine,
and in a really important aspect. This would cause major problems.

Generally, I like providing explicit interfaces. For module procedures
they are automatic, and I do tend to write interface bodies for things
like C functions that I call (unless the function has some C-ism
that can't be expressed in a Fortran 90/95 interface - such as a void*
argument). But if you write an interface body, it does have to match
the actual procedure. Better to not write one at all than to write one
that is incorrect.

I'm now betting that the OP doesn't understand the distinctions
between assumed shape and other kinds of dummy arrays well enough
to realize that his interface body doesn't correctly describe
the subroutine in question. I could be off base, but that's what
strikes me at the moment.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Steven G. Johnson

2004-10-13, 3:59 pm

"James Van Buskirk" <not_valid@comcast.net> wrote...
> Hey, Steven. You should have quoted the INTERFACE statement as well
> because that's what says the given code is doomed, doomed, doomed.
> IN & OUT are declared in Fortran as assumed-shape dummies in (what I
> would imagine is) you C function returning void.

Thanks James! Yes, that would do it -- FFTW's Fortran interface
certainly can't handle any kind of shape "descriptors" in place of raw
arrays.

(I'm afraid my knowledge of Fortran is mostly limited to F77, which
makes responding to these kinds of F90 queries difficult...I was
hesitant to chip in, but it seemed like there was confusion in the
thread about what FFTW expected to be passed.)

> P.S. have you had a chance to take a look at my convolution stuff?

No, I haven't seen it. I got very busy this fall with a new faculty
position, moving, etcetera, and am just starting to settle in. I took
a quick look at your web page just now, and it looks like you're
making nice progress!

Cordially,
Steven G. Johnson

Jacek

2004-10-13, 8:56 pm

>
>
> Thanks James! Yes, that would do it -- FFTW's Fortran interface
> certainly can't handle any kind of shape "descriptors" in place of raw
> arrays.
>
> (I'm afraid my knowledge of Fortran is mostly limited to F77, which
> makes responding to these kinds of F90 queries difficult...I was
> hesitant to chip in, but it seemed like there was confusion in the
> thread about what FFTW expected to be passed.)
>

Thanks, people! It solves my problem. I should use arrays with
specified number of rows and columns when describing interface to
FFTW functions.

Many thanks,
Jacek

B52B

2004-10-14, 3:59 pm

Richard E Maine wrote:
> B52B <B52B@pl.pl.pl> writes:
>
>
>
>
>
> Intent does not matter. It can't matter because dummy arguments with
> intent specified do not trigger the requirement for explicit interfaces.
> Thus, the compiler doesn't necessarily know the intents when compiling
> the calling program.

I agree with you. It should doesn't matter but I remember problem where
gemm (from BLAS) had interface with intents in main program. There was
argument mismatch error during compilation.
So lack of intents was the way to overcome the problem.

>
> I'm now betting that the OP doesn't understand the distinctions
> between assumed shape and other kinds of dummy arrays well enough
> to realize that his interface body doesn't correctly describe
> the subroutine in question. I could be off base, but that's what
> strikes me at the moment.
>
Sure. He should take f77 interface as a base.

Regards,
B52B

Ian Bush

2004-10-14, 3:59 pm

Janne Blomqvist wrote:

> On 2004-10-13, B52B <B52B@pl.pl.pl> wrote:
>
> In my experience, ifort understands gcc/g77 calling conventions. I
> doubt this is the source of the problems.

Hmmmm, I knew I remembered something. ifort/ifc and gcc are fine,
but there can be some problems with ifort/ifc and g77. See

http://www.intel.com/software/produ...atibility60.pdf

That said that article is a little old, and I do not know the current situation,

Ian

Steven G. Johnson

2004-10-15, 3:56 am

Ian Bush <I.J.Bush@nospam.dl.ac.uk> wrote...
> Hmmmm, I knew I remembered something. ifort/ifc and gcc are fine,
> but there can be some problems with ifort/ifc and g77. See
>
>http://www.intel.com/software/produ...atibility60.pdf
>
> That said that article is a little old, and I do not know the current situation,

The only incompatibility that this article mentions is the
name-mangling difference (g77 adds an extra underscore to identifiers
containing an underscore, an annoyance that apparently dates back to
f2c). That's still the case (at least for g77).

It's not really relevant for FFTW, however, since FFTW includes F77
interfaces with both single and double underscores appended.

Sponsored Links