How are temp arrays declared in subroutines managed?

[steve_o_]

I'm just learning Fortran but have not been able to find the definitive answer to the question of how arrays declared in subroutines are managed

For example see below (this is a Fortran question not specific to MKL), I declare the pivot, AF, INV, r, c and work locally in the subroutine. It looks a bit fishy to me in the context that I'm using allocate to reserve space in my main program but here I'm just using a passed value of lda. (However, it does appear to work or at least do something, I'm not totally convinced about the answer I get back is correct though!)

I don't like the idea of declaring all arrays in the main program and passing them around, I would think it adds a lot of overhead ? Also I did try dynamic allocation in the subroutine but I wasn't sure if memory was held between calls or if I was leaking memory or worse

What's the best way to create temp array storage in a subroutine ?

SUBROUTINE GetInverseMatrixTrace(Matrix, lda, trace)
INTEGER lda
REAL(8), DIMENSION(lda,lda), INTENT(IN) :: Matrix
REAL(8) trace
!! local vars
Integer info, index, workSize
INTEGER, DIMENSION(lda) ::  pivot
REAL(8), DIMENSION(lda,lda) :: AF, INV
REAL(8), DIMENSION(lda) :: r,c
REAL(8), DIMENSION(10000) :: work ! should be at least 4 x totalweights
REAL(8)  amax, rowcnd, colcnd
info = 0
workSize = 10000 ! compiler recommended minimum is 4x lda size
AF = Matrix
call dgetrf( lda, lda, AF, lda, pivot, info )
call dgeequ( lda, lda, AF, lda, r, c, rowcnd, colcnd, amax, info )
INV = AF
call dgetri( lda, INV, lda, pivot, work, workSize, info ) ! calc inverse Hessian for trace calc
trace = 0.0d0
do index=1, lda
	trace = trace + INV(index,index)
enddo
end SUBROUTINE GetInverseMatrixTrace

[mecej4 · Accepted solution]

As in other compiled languages such as C and Pascal, formal arguments that are arrays have array declarations not to allocate memory afresh, but to give the compiler the information needed to compute the addresses of array elements. In your example, Matrix(:,:) is an argument array, and you have to go back in the callers' chain to find its actual allocation/declaration. The local arrays Pivot, AF, inv, etc., are allocated on the stack and, if your OS or other reasons make large consumption of stack undesirable, you can allocate them on the heap instead. Most compilers pass arrays by reference most of the time, so there is next to no overhead in doing this.

If you want the values of local variables to be preserved, you have to give them the SAVE attribute. If you want to do dynamic allocation, see the ALLOCATE and DEALLOCATE statements.

Moving on to a different aspect, note that there may be short-cut methods to computing the trace of the inverse matrix that do not require that the full inverse be formed. See http://mathoverflow.net/questions/46553/fast-trace-of-inverse-of-a-square-matrix .

[steve_o_]

Thanks, makes me a bit more comfortable ...

I'm already doing lots of dynamic allocation elsewhere so no probs. - just scratch working arrays I was worried about, I only want them temporary

Thanks for the inverse hessian trace tip, I did look before but obviously missed it, hence I wrote my own

slightly off topic, you haven't come across a standard deviation and or chi square error routine have you ;-), rolled my own simplistic subroutine there as well but I'm wary of floating point errors with small differences so I was wondering if there was a tried and tested safe algorithm I may have missed

[jimdempseyatthecove]

Caution,

When a subroutine is not recursive, or -recursive not used, or -openmp .and. -heap-arrays not use...

Then it is undefined as to if local arrays are stack located or SAVE. Typically they are SAVE.

While this will not make a difference in single threaded programs, it will make a difference for routines compiled without recursive attribute, or -recursive , or -openmp .and. -heap-arrays, and then subsequently linked with a threaded application. IOW all threads would tromp on the same implicitly SAVE'd local arrays.

Jim Dempsey

[steve_o_]

Hmm good catch Jim

May explain why I was getting a seg fault when I was running optimized code when I run with -check all -traceback -g  -debug all  -heap-arrays it runs OK most of the time ( I get occasional NaN which I have put down to floating point error or had I been using C pointers a stray pointer error)

[mecej4]

you haven't come across a standard deviation and or chi square error routine have you 

Intel Fortran includes the MKL library, which provides some summary statistics routines. The Intel Performance Primitives (IPP) library also contains a number of statistics routines. In particular, see MeanStdDev(), which computes the mean value and the standard deviation value of a vector. For the Chi-Sq. distribution, you may look at other libraries such as IMSL, NAG and specific packages in Netlib/TOMS.

[steve_o_]

Thanks, already using MKL but knowing what to look for helps ;-), I'm interested in IPP but it needs to be OS X compat. Just bought ifort so anything else may need to wait, especially if it costs - think IPP is $$, I found a discontinued version, may look through that or download the OS X demo. No big deal just looking for optimised algorithms. 

regards

Steve

[Steven_L_Intel1]

IPP is not discontinued, but it is provided only with Intel C++, not separately anymore.

[mecej4]

Steve O., if you have not explored available software for statistical calculations, you may be pleasantly surprised after visiting the StatLib site, http://lib.stat.cmu.edu/ .

[steve_o_]

many thanks,

I previously had a quick look through the list on wikipedia numeric libs, I hadn't previously seen the one you mention

regs

Steve