moving away from EQUIVALENCE for in-place FFT of large arrays

[Alexis_R_]

Hello,

I am currently working on a project where EQUIVALENCE statements are used, mainly for the purpose of in-place, REAL to COMPLEX Fourier transform of very large arrays. I would like to move away from EQUIVALENCE, mainly because it stops me from making the arrays in question ALLOCATABLE.

The problem is also described here:

"While I know that the use of EQUIVALENCE and COMMON blocks has been heavily discouraged for some time, I have yet to find an equivalent (no pun intended) solution for their use with regards to computing fourier transforms in-place. This is a common practice in spectral CFD codes in which the data is either all REAL or all COMPLEX at any given time and is converted back and forth very efficiently using inline FFTs such as the canned FFTW libraries. For these cases, being able to EQUIVALENCE both a real and complex array not only avoids extra storage, but also eliminates unnecessary copies upon exit of the FFT algorithm."

I should point out that using TRANSFER to copy the data is not an option in my case, because of the dimensions of the arrays.

Is there an elegant, or at least compliant/portable, solution to this problem? I have read examples involving TRANSFERing derived pointer types but have yet to fully understand & get this to work for my purpose.

Any guidance would be much appreciated,
Alexis

[jimdempseyatthecove]

I can point you in the right direction (sorry I don't have example handy). Fortran now has interoperability intrinsic functions for converting call parameters from C to Fortran and Fortran to C. You can convert ("cast") the base cell of your array into a C pointer then immediately convert it back to a Fortran array descriptor of other type. You should be able to flip back and forth quite easily with little overhead (the time it takes to fill in the new array descriptor).

Jim Dempsey

[joseph-krahn · Accepted solution]

Quoting - A. Rohou

Hello,

I am currently working on a project where EQUIVALENCE statements are used, mainly for the purpose of in-place, REAL to COMPLEX Fourier transform of very large arrays. I would like to move away from EQUIVALENCE, mainly because it stops me from making the arrays in question ALLOCATABLE.

The problem is also described here:

"While I know that the use of EQUIVALENCE and COMMON blocks has been heavily discouraged for some time, I have yet to find an equivalent (no pun intended) solution for their use with regards to computing fourier transforms in-place. This is a common practice in spectral CFD codes in which the data is either all REAL or all COMPLEX at any given time and is converted back and forth very efficiently using inline FFTs such as the canned FFTW libraries. For these cases, being able to EQUIVALENCE both a real and complex array not only avoids extra storage, but also eliminates unnecessary copies upon exit of the FFT algorithm."

I should point out that using TRANSFER to copy the data is not an option in my case, because of the dimensions of the arrays.

Is there an elegant, or at least compliant/portable, solution to this problem? I have read examples involving TRANSFERing derived pointer types but have yet to fully understand & get this to work for my purpose.

Any guidance would be much appreciated,
Alexis

There are a few solutions to this problem. The best solution would be for the Fortran committee to not be so down on EQUIVALENCEs. They are very useful in some situations. I often feel that the new Fortran standards try to force people to write better code by invalidating "bad" constructs, even though they can sometimes be useful, and people can always write bad code in any case. They don't even allow equivalences for interop with C unions.

Sequence-association is probably the best way to make real and complex types equivalent, and use only valid Fortran. Sequence associaton rules, along with a defined order for the real and imaginary parts of complex types, allows you to use a dummy argument as a real/complex equivalence.

There are two disadvantages of sequence-association. 1) it requires an EXTERNAL subroutine interface, which is only a problem if you are a module 'purist'. 2) it is only valid as a dummy argument, not a global variable.

It is possible to cast pointers to allow global access of an array. The easiest way is to use ISO_C_Binding using a C pointer as an intermediate, something like this:

type(C_ptr) :: p
p = C_loc(complex_array)
call C_F_pointer(p,real_array,[2*size(complex_array,1)])

The disadvantage of pointers is that they allow memory aliasing, which can prevent some code optimizations. However, passing the pointer as a dummy argument in a routine without a pointer or target attribute removes the aliasing, and it becomes the programmers job to ensure that the routine does not access aliased memory. If you want to access the same array as both complex and real in the same routine, the target or pointer property is required to avoid problems.

[Alexis_R_]

Thanks a lot for your answer. I am not sure I understand everything in there.

Firstly, I have been using what I think is what you refer to as sequence-association, e.g. calling a subroutine with argument N being a REAL(1:N) where the subroutine's Nth dummy argument is a COMPLEX(1:N/2). However, if an interface is present (or if I use -gen-interfaces, which I like a lot), ifort will throw an error, which is fair enough. Am I missing something? What I'm missing probably has to do with me not knowing what you mean by an EXTERNAL interface.

Secondly, the intrinsic ISO_C_BINDING module way does look promising - I'm guessing this is what Jim Dempsey was pointing me towards. Could you elaborate on the circumstances under which aliasing might occur? Would ifort warn me with an "array temporary" warning in those cases? I wanted to illustrate the aliasing problem with the code below, but both sub1 and sub2 modified both the complex and real arrays, whereas I thought aliasing might happen in sub1, causing only the real array to be written to? Again, what am I missing?

Thanks :)

[cpp]program cbinding

	use, intrinsic :: iso_c_binding

	implicit none

	complex	:: COMPLEX_ARRAY(2)
	type(C_ptr) :: CPTR
	real, pointer :: PTR(:)

	COMPLEX_ARRAY(1) = (1.0,0.0)
	COMPLEX_ARRAY(2) = (2.0,0.0)

	CPTR = C_loc(COMPLEX_ARRAY)
	call C_F_pointer(CPTR,PTR,[2*SIZE(COMPLEX_ARRAY,1)])

	write(*,*) 'program'
	write(*,*) COMPLEX_ARRAY
	write(*,*) PTR

	call sub1(COMPLEX_ARRAY,PTR)
	call sub2(COMPLEX_ARRAY,PTR)

	contains

		subroutine sub1(COMPLEX_ARRAY,REAL_ARRAY)
			implicit none
			complex :: COMPLEX_ARRAY(:)
			real	:: REAL_ARRAY(:)

			REAL_ARRAY(2) = 2.99

			write(*,*) 'sub1'
			write(*,*) COMPLEX_ARRAY
			write(*,*) REAL_ARRAY


		end subroutine sub1


		subroutine sub2(COMPLEX_ARRAY,REAL_ARRAY)
			implicit none
			complex, target	:: COMPLEX_ARRAY(:)
			real, pointer	:: REAL_ARRAY(:)

			REAL_ARRAY(4) = 2.99

			write(*,*) 'sub2'
			write(*,*) COMPLEX_ARRAY
			write(*,*) REAL_ARRAY

		end subroutine sub2

end program cbinding
[/cpp]

The output I get:

[plain] program
 (1.000000,0.0000000E+00) (2.000000,0.0000000E+00)
   1.000000      0.0000000E+00   2.000000      0.0000000E+00
 sub1
 (1.000000,2.990000) (2.000000,0.0000000E+00)
   1.000000       2.990000       2.000000      0.0000000E+00
 sub2
 (1.000000,2.990000) (2.000000,2.990000)
   1.000000       2.990000       2.000000       2.990000    
[/plain]

[jimdempseyatthecove]

Aliasing is not creating an array temporary.

Aliasing means you have two or more variable names (or arrays) that are the same or overlap. When the compiler (or programmer) does not know of the aliasing and assumes (or is told) no aliasing then incorrect or unexpected results occure.

Assume you have an array Ary(100)
Assume you have a subroutine to add two arrays

subroutine doAdd(A, B, n)
integer :: n
real :: A(n), B(n)
A = A + B
end subroutine doAdd

Then call

call doAdd(Ary(2:51), Ary(1:50), 50)

If the subroutine doAdd assumes (or is told) A and B are not aliased (in this case overlapped),the compilerwould likely use the SSE3 instructions. In this case incorrect results are the outcome.

Same true with multiple scalars or multiple dummy's referencing the same variable when the subroutine assumes or requires they be different.

Jim Dempsey

[Alexis_R_]

Dear Jim,

Thanks for the explanation!

[joseph-krahn]

By "EXTERNAL interface", I mean a plain Fortran77 style call with no explicit interface. In that case, sequence association is well defined by the standards, so it is not just a 'hack' that happens to work. Unfortunately, it is not allowed in explicit interfaces, which is what you get for any module procedure, and with ifort's auto-generated interfaces. With Fortran being a numeric language, I am a bit surprised that there was not a better mechanism for complex/real associations.

Intel Fortran has some extensions that allow explicit interfaces to be called with different argument types using sequence association rules. I think that will enable -gen-interfaces to work, but I usually avoid extensions, and the ISO_C_Binding is generally working in new compilers.

Joe

[nvaneck]

Quoting - krahn@niehs.nih.gov

By "EXTERNAL interface", I mean a plain Fortran77 style call with no explicit interface. In that case, sequence association is well defined by the standards, so it is not just a 'hack' that happens to work. Unfortunately, it is not allowed in explicit interfaces, which is what you get for any module procedure, and with ifort's auto-generated interfaces. With Fortran being a numeric language, I am a bit surprised that there was not a better mechanism for complex/real associations.

Intel Fortran has some extensions that allow explicit interfaces to be called with different argument types using sequence association rules. I think that will enable -gen-interfaces to work, but I usually avoid extensions, and the ISO_C_Binding is generally working in new compilers.

Joe

These answers are way to complicated for me. I read very large data sets into MicrOsiris and convert values to the required mode for each program, To do this, I equivalence 6 arrays of various types, including character, so any of the analysis and data management programs can refer to the proper mode and index to function properly. The data can be mixed mode in the converted arrays as needed, particularly when transfering data to an output file.There seems to be no efficient way to avoid EQUIVALENCE, nor do I see a need for the standard to get rid of it. Maikng it obsolesent if fine, but I hope they never drop it.

[Steven_L_Intel1]

I don't forsee that EQUIVALENCE will be deprecated in the language any time soon (if ever). Even if it is, compilers will continue to support it.