Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Petr_Parik
Beginner
636 Views

Allocatable arrays of different type and EQUIVALENCE

Jump to solution

I work on a very old code which relies on two big working arrays IW and RW, which are equivalenced in the main program and then used directly or propagated to called subroutines like this:

      program test
        implicit none
        external inner
        integer(4),parameter :: li = 100
        integer(4) iw(li)
        real(8)    rw(li/2)
        equivalence (iw,rw)
        rw = -10.20
        print *,iw
        print *,rw
!       call someroutine(iw,li,rw,lr) will also work
      end

When I replace the static arrays with dynamically allocated, then I of course cannot equivalence the arrays. I came up with this code:

      program test
        implicit none
        external inner,outer
        integer(4),parameter :: li = 100
        call outer(inner,li)
      end
      subroutine outer(proc,li)
        implicit none
        external proc
        integer(4),intent(in)  :: li
        integer(4),allocatable :: iw(:)
        allocate(iw(li))
        call proc(iw,li,iw,li/2)
        deallocate(iw)
      end
      subroutine inner(iw,li,rw,lr)
        implicit none
        integer(4),intent(in)    :: li,lr
        integer(4),intent(inout) :: iw(li)
        real(8),intent(inout)    :: rw(lr)
        rw = -10.20
        print *,iw
        print *,rw
!       call someroutine(iw,li,rw,lr) will also work
      end

So here the main program is reduced to a call to a helper subroutine outer, which allocates a single array and then passes it as two arrays to subroutine inner, which in turn contains the code of the original main program. The reason the compiler does not complain about the third parameter to inner not being a real(8) array is that in outer, inner is hidden in proc as an external subroutine with arbitrary arguments.

Is there any other, maybe more simpler or clean, way?

0 Kudos
1 Solution
FortranFan
Honored Contributor I
291 Views

@Petr_Parik ,

Also, if your interest in Fortran is limited to doing the least with this legacy codebase to get it to work with a larger problem size beyond the stack limits, you are essentially done by using the approach with the IW allocatable+target attributes, RW as pointer, and C_F_POINTER memory  mapping shown above.  No need to read any further, just "let the sleeping dogs lie" with the legacy code.

But if that ain't the case and you need to do more with Fortran, you may want to first start with a close review of a book such as this: 

https://www.amazon.com/Fortran-Scientists-Engineers-Stephen-Chapman-ebook/dp/B06XCTY8KR

And pick up on modern Fortran concepts that can open up a world of options to write better new code and improve upon existing/legacy codebases for efficient and manageable applications in scientific and technical computing.

View solution in original post

32 Replies
Petr_Parik
Beginner
167 Views

@FortranFan ,

One last question, just to be sure. To solve the issue with the array argument, it is sufficient to add just a colon after the index to have a standard-compliant code? Like so:

CALL BAR(...,RW(LBFOO:),...)

 

jimdempseyatthecove
Black Belt
160 Views

Isn't RW(...) argument a subsection of the blob?

As such, you would need to pass the lower bound (LBFOO) as well as the upper bound (UBFOO).

Using the RW(LBFOO:) would pass a slice starting at LBFOO and extending to the end of the blob.

Jim Dempsey

Petr_Parik
Beginner
153 Views

@jimdempseyatthecove ,

Remember that this is ancient code where arrays are always passed by their first element and length:

CALL FOO(RW(LBBAR),LBAR,IW(LBBAZ),LBAZ)
.
.
.
END
SUBROUTINE FOO(A,LA,IA,LIA)
DIMENSION A(LA),IA(LIA)
.
.
.

However, the moment RW becames POINTER, the compiler starts to complain about RW(LBBAR) not being an array.

Adding a colon in the argument RW(LBBAR:) makes the code compilable again. For all intents and purposes, this should have no effect at runtime.

jimdempseyatthecove
Black Belt
144 Views

While you can use (LBBAR:) today, sometime in the future, you or your replacement may be updating the code and which at that time may introduce errors (or at least confusion). While adding ":" makes sense to you today, it might not make sense to someone else later.

So you have three choices:

Use ...,RW(LBBAR:),LBAR.... ! ":" indicates to end of blob, fix later

or use RW(LBBAR:UBBAR),LBAR,...! LBAR is passed for legacy purposes, fix later

or use RW(LBBAR:LBBAR+LBAR-1),LBAR, ... ! add UBBAR later, then remove LBAR on next revision.

Jim Dempsey

 

jimdempseyatthecove
Black Belt
139 Views

When you revise the code (some time later), you might consider using Defining Generic Names for Procedures (intel.com)

This way you can have your subroutine FOO become a generic name with two interfaces: old/legacy argument list, and new updated argument list.

Jim Dempsey

FortranFan
Honored Contributor I
126 Views

Jim,

See my comment just prior to yours and the option therein - here's a link to it.

The use of additional variables of POINTER attribute for the work arrays of interest to OP - IA1, IA2, RA1, RA2, etc. - in this legacy code is more safe and readable, maintainable in my experience.

jimdempseyatthecove
Black Belt
118 Views

Petr_Parik,

A word of caution about being overly zealous with allocating vary large array.

At time of allocation of a huge blob, the allocation may succeed in obtaining virtual addresses, but can fail some time later at first touch when it is found that the page file has insufficient space. This failure (page file full) is not trappable.

The is a caution not to request a few terabytes of memory for your buffer.

Jim Dempsey

FortranFan
Honored Contributor I
147 Views
@Petr_Parik wrote:
.. One last question, just to be sure. To solve the issue with the array argument, it is sufficient to add just a colon after the index to have a standard-compliant code? Like so:

CALL BAR(...,RW(LBFOO:),...)

@Petr_Parik ,

While you can use an array section notation to solve the issue, note you may possibly face other issues and that's why I didn't suggest it earlier.  Note it's tough to list what problems you may run into what without knowing your entire code but it can be a performance hit in the form of an array temporary in CALLs to other fault scenarios during run-time.

If you're OPEN to the notion of making changes on the CALLER side like what you show with "CALL BAR(..,RW(..:), ..", then my suggestion will be to use the POINTER facility in modern Fortran to setup SPECIFIC objects with the POINTER attribute for EACH of your "work" ARRAYs.  See an illustration of this with your SAMPLE program where you will see 4 SPECIFIC objects for your 4 arrays - IA1, IA2, RA2, and RA2.  Given your legacy code with all the existing "book-keeping" toward the indices and array sizes, you'll see it's straightforward to do this and it WILL work with Intel Fortran, gfortran, and pretty much all the standard-conforming Fortran compilers out there.  I show this Intel Fortran and gfortran below.

      PROGRAM TEST_PTR
      USE ISO_C_BINDING
      IMPLICIT INTEGER(4) (I-N), REAL(8) (A-H,O-Z)
      PARAMETER (LI = 100)
      ALLOCATABLE IW(:)
      TARGET IW
      REAL(8), POINTER :: RW(:)
C Additional pointers for array argument to WORK* subprograms
      INTEGER(4), POINTER :: IA1(:)
      INTEGER(4), POINTER :: IA2(:)
      REAL(8), POINTER :: RA1(:)
      REAL(8), POINTER :: RA2(:)
      ALLOCATE(IW(LI))
      CALL C_F_POINTER(C_LOC(IW),RW,(/LI/2/)) ! RW => IW
      LBIA1 =  1; LIA1 = 10
      LBIA2 = 21; LIA2 = 10
      LBRA1 =  6; LRA1 = 5
      LBRA2 = 16; LRA2 = 5
C Set up array arguments
      IA1(1:LIA1) => IW(LBIA1:)
      IA2(1:LIA2) => IW(LBIA2:)
      RA1(1:LRA1) => RW(LBRA1:)
      RA2(1:LRA2) => RW(LBRA2:)
      CALL WORK1(IA1,LIA1,IA2,LIA2,
     &           RA1,LRA1,RA2,LRA2)
      CALL WORK2(IA1,LIA1,IA2,LIA2,
     &           RA1,LRA1,RA2,LRA2)
      END
      SUBROUTINE WORK1(IA1,LIA1,IA2,LIA2,RA1,LRA1,RA2,LRA2)
      IMPLICIT INTEGER(4) (I-N), REAL(8) (A-H,O-Z)
      DIMENSION IA1(LIA1),IA2(LIA2),RA1(LRA1),RA2(LRA2)
      IA1 = 1
      IA2 = 2
      RA1 = .123456789
      RA2 = .987654321
      END
      SUBROUTINE WORK2(IA1,LIA1,IA2,LIA2,RA1,LRA1,RA2,LRA2)
      IMPLICIT INTEGER(4) (I-N), REAL(8) (A-H,O-Z)
      DIMENSION IA1(LIA1),IA2(LIA2),RA1(LRA1),RA2(LRA2)
      PRINT *,"IA1 ="
      PRINT *,IA1
      PRINT *,"IA2 ="
      PRINT *,IA2
      PRINT *,"RA1 ="
      PRINT *,RA1
      PRINT *,"RA2 ="
      PRINT *,RA2
      END

Build and execute with Intel Fortran and gfortran:

C:\Temp>ifort test_ptr.for
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.1.2 Build 20201208_000000
Copyright (C) 1985-2020 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.26.28806.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:test_ptr.exe
-subsystem:console
test_ptr.obj

C:\Temp>test_ptr.exe
 IA1 =
           1           1           1           1           1           1
           1           1           1           1
 IA2 =
           2           2           2           2           2           2
           2           2           2           2
 RA1 =
  0.123456791043282       0.123456791043282       0.123456791043282
  0.123456791043282       0.123456791043282
 RA2 =
  0.987654328346252       0.987654328346252       0.987654328346252
  0.987654328346252       0.987654328346252

C:\Temp>gfortran test_ptr.for -o gcc-test_ptr.exe

C:\Temp>gcc-test_ptr.exe
 IA1 =
           1           1           1           1           1           1           1           1           1           1
 IA2 =
           2           2           2           2           2           2           2           2           2           2
 RA1 =
  0.12345679104328156       0.12345679104328156       0.12345679104328156       0.12345679104328156       0.12345679104328156
 RA2 =
  0.98765432834625244       0.98765432834625244       0.98765432834625244       0.98765432834625244       0.98765432834625244

C:\Temp>

 

Petr_Parik
Beginner
93 Views

@FortranFan ,

I know sub-array pointers would work. Unfortunately, that is even more error-prone work than adding ":" to calls, as I have found. So I have settled with a slightly bloated but simple "loader" approach:

 

      PROGRAM LOADER
C     allocation/equivalence stuff ...
      CALL FOO(IW,LI,RW,LI)
      END
      SUBROUTINE FOO(IW,LI,RW,LR)
      DIMENSION IW(LI),RW(LR)
C     original PROGRAM FOO code with almost no changes
C     'RW' behaves like a normal array, no need to modify calls
C     ...
      END

 

 

@jimdempseyatthecove ,

Adding ":" or upper bound would be too much error-prone work, see above.

However, interfaces for incremental refactoring are an excellent idea, thank you.

Regarding the virtual memory; allocating a large blob then failing later at runtime due OOM is no different to allocating arrays as needed and then still failing later at runtime when virtual memory is depleted. But fear not, LI is INTEGER(4) so we are talking 8 GiB blob max. Some new programs have 64-bit LI, but even in that case the maximum needed so far has been less than the physical memory (128 GiB). Generally, if I need to do computations, I need to make sure nobody else is using the computer.

jimdempseyatthecove
Black Belt
78 Views

>> allocating a large blob then failing later at runtime due OOM is no different to allocating arrays as needed and then still failing later at runtime when virtual memory is depleted

With the exception that one can test for failure at ALLOCATE (assuming the heap isn't larger than remaining space in page file).

>> But fear not, LI is INTEGER(4) so we are talking 8 GiB blob max

From your description (which I may have misunderstood) is your code allocates once and anything becoming unused is lost. Presumably the allocations reach a stasis with a not unreasonable amount of dead space in the blob.

Jim Dempsey

Petr_Parik
Beginner
71 Views

@jimdempseyatthecove ,

>> With the exception that one can test for failure at ALLOCATE (assuming the heap isn't larger than remaining space in page file).

You are right, but if we distinguish only between a successful computation and a failed computation, it does not matter if the program fails gracefully.

>> From your description (which I may have misunderstood) is your code allocates once and anything becoming unused is lost. Presumably the allocations reach a stasis with a not unreasonable amount of dead space in the blob.

As you pointed out earlier, modern operating systems allocate only the memory accessed by the program. But yes, there can be several gigs of unused space in the blob.

FortranFan
Honored Contributor I
354 Views
@Petr_Parik wrote:

.. If I am not mistaken this problem is nonexistent in C where you can simply recast any variable to the desired type while passing it (by reference) to a function (subroutine).

 

@Petr_Parik ,

To the extent you're using modern Fortran compilers such as Intel Fortran (IFORT) as part of Intel oneAPI or recent incarnations of GCC/gfortran, NAG compiler, etc. - all of which work with a companion C processor, you can take advantage of modern Fortran standard - Fortran 2018 - and its facilities to start introducing INCREMENTAL REFACTORING such as dynamic memory and associated management and build toward transforming such legacy codebase into a modern app.

Shown below is a simple prelude you can review - note the EQUIVALENCE statement has been replaced with a C processor based trick that maps the ALLOCATED memory to an array of floating-point type (REAL with kind of R8) with a POINTER attribute.  Then with suitable use of indices in this legacy codebase - assuming they are set up to keep the integer and real "data" of this program separate, you should be able to work this with the rest of your code.  Note the "editorial" liberty I've taken is to replace the non-standard "INTEGER*4" and "REAL*8" types and kinds with a tab more palatable options from the Fortran standard.

      PROGRAM SAMPLE
      USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : I4 => INT32,
     &                                          R8 => REAL64, 
     &                                          STDIN => INPUT_UNIT,
     &                                          STDOUT => OUTPUT_UNIT
      USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_LOC, C_F_POINTER 
      IMPLICIT INTEGER(KIND=I4)(I-N), REAL(KIND=R8)(A-H,O-Z)
      INTEGER(I4), ALLOCATABLE, TARGET :: IW(:)
      REAL(R8), POINTER :: RW(:)
      WRITE(STDOUT,*) "Enter program size parameter value: LI = "
      WRITE(STDOUT,*) "(for demo, enter an integer value >= 6)"
      READ(STDIN,*) LI
      ALLOCATE( IW(LI), STAT=ISTAT )
      IF ( ISTAT /= 0 ) THEN
         WRITE(STDOUT,*) "Error allocating memory: ISTAT=", ISTAT
         STOP
      END IF
      IW = 0
      CALL C_F_POINTER( CPTR=C_LOC(IW), FPTR=RW, SHAPE=[ LI/2 ] )
C Assuming program logic with indexes ensures integer and real data
C do not overlap
      LBEGIW = 1 ; LENDIW = 2
      LBEGRW = 2 ; LENDRW = 3
      IW(LBEGIW) = 42 ; IW(LBEGIW+1) = 43
      RW(LBEGRW) = 10.0_r8 ; RW(LBEGRW+1) = 20.0_r8 
      WRITE(STDOUT,*) "Integer data: ", IW(LBEGIW:LENDIW)
      WRITE(STDOUT,*) "Real data: ", RW(LBEGRW:LENDRW)
      RW => null()
      END

 

Upon execution with Intel Fortran,

C:\temp>ifort /standard-semantics p.for
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.1.2 Build 20201208_000000
Copyright (C) 1985-2020 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.27.29112.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
 Enter program size parameter value: LI =
 (for demo, enter an integer value >= 6)
6
 Integer data:  42 43
 Real data:  10.0000000000000 20.0000000000000

C:\temp>
Reply