- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
Link Copied
- « Previous
-
- 1
- 2
- Next »
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:),...)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@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:),...)
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>
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
>> 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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
>> 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@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).
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>
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page
- « Previous
-
- 1
- 2
- Next »