using real(A[,kind=]) in passing a single precision array as a double precision to a subroutine in Fortran 90/2003

Question

I need to perform both single precision and double precision arithmetics on a variable in different parts of my code. So basically, I declare the variable as single precision first. Then I call subroutine sub_a which makes use of double precision version of the variable and performs double precision operations on that:

program main

 implicit none
 integer,parameter :: single = selected_real_kind(p=6,r=37)
 integer,parameter :: double = selected_real_kind(p=15,r=307)

 real(single),allocatable,dimension(:) :: A
 real(double),allocatable,dimension(:) :: B

 allocate(A(3),B(3))
 A=2 ! single precision
 A=A+3 ! single precision

 print '(a,1x,3(f20.15))','sqrt(A),single:',sqrt(A)
 print '(a,1x,I15)','mem_address of A before sub_a:',loc(A)

 call sub_a(real(A,kind=double),B) ! double precision
 print '(a,1x,3(f20.15))','sqrt(A),double:',B

 contains
   subroutine sub_a(a,b)
     real(double),dimension(:),intent(in) :: a
     real(double),dimension(:),intent(inout) :: b

     print '(a,1x,I15)','mem_address of A in sub_a:',loc(a)

     b=sqrt(a)

   end subroutine sub_a

end program main

As seen in the code, I also obtained the memory address of A prior to calling sub_a and the version of A inside sub_a and they are expectedly different.

My questions are:

1. Is the version of A inside the sub_a allocated in the heap memory so I should not be worried about the size limitation?
2. is there any potential issue/bug in writing this example?
3. is there a better way of capturing the purpose described in this example, specially for larger size arrays?

Many Thanks

Update: I haven't experienced any Memory issue for very large arrays, when using gfortran4.6 / ifort13.1 as compiler. I plan to use suggestion by @innoSPG as an alternative approach.

Answer 1

The temporary array at the call site can be made at the stack or at the heap, it is implementation dependent. Compilers usually have options to control the behavior.

With Intel Fortran it the option -heap-arrays n (n cannot be ommited or too low, because performance would be bad).

Gfortran would put them on the heap automatically if it doesn't know the size beforehand. It is good to use -fstack-arrays for better performance (it is included in -Ofast).

You will find similar options in other compiler's documentation.

I would not fear and avoid the temporary arrays at all cost. The possibility of much shorter (and more readable!) code is sometimes more important. Many codes have portions that are done only once during initialization or final output where the performance penalty is often irrelevant. I personally use it to make my code cleaner in these parts (not in the computational core).

Answer 2

by the nature of the call, the version of A that you have in sub_a is a temporary array created by a piece of code included by the compiler. However, if you will manipulate very large arrays, it is not a good idea.

For question 2. to my knowledge, there is no bug. The only issue is the temporary array that may be a problem if you have large arrays and limited memory on your system.

For question 3. In the case there is a memory issue, you can write sub_a to accept simple precision, and then convert each element in sub_a before using it in computation.