why.n0t
why.n0t

Reputation: 432

Share allocatable Arrays

I have some allocatable arrays which I need to share between some subroutines. I usually would just pass them as arguments or maybe write everything in a Module, but I'm afraid this isn't possible in my situation.

I only write some own subroutines and use subroutines provided and described by an FEM-Solver. So i cannot alter the arguments of this subroutines or wrap them in a Module.

As far as i know it also isn't possible to Build common blocks with array of unknown size at compile time.

Is there something else to pass my arrays?

Update:
At the moment my program environment looks like this:

I have a subroutine, provided by the FEM-program, which is called after each increment, this calls several of my subroutines where I compute some values for each node or for a subset of those.

To display these values in the post-Simulation, i have to pass them to another subroutine. This subroutine is called by the FEM-solver for each node at the end of the increment. So shifting my code to this Subroutine would produce a lot of overhead.

My idea is to compute the values once, store the Values in an array and pass this array to the second subroutine where they will be written to the database of the computation.

Update
Some Pseudo-code:
Assumed from program behaviour:

 Program FEM-solver
     *magic* 
     call ENDINC(ar1,ar2)
     *something* 
     do NodeID=1,Sum_Of_Nodes
        do valueID=1,Sum_Of_User_Computed_Values !(defined in preprocessing)
           call nodeval(NodeID,valueID,Value,ar3,...,arN)
        end do
     end do
     *voodoo* 
 end program FEM-solver    

Written and working:

Subroutine ENDINC(ar1,ar2)
  *Computation of some node values*
  *Calling of own Subroutines, which compute more values*
  *Writing an array with results values for some/each node(s)*
   nodersltArr(NodeID,rslt)=*some Value*
end Subroutine ENDINC

Needed, writng the computed Values to the Node solution database:

Subroutine nodeval(NodeID,valueID,Value,ar3,...,arN)  
  *called for each NodeID and valueID*
   value=noderslArr(NodeID,valueID)
end subroutine nodeval

Upvotes: 3

Views: 3713

Answers (3)

M. S. B.
M. S. B.

Reputation: 29401

You can pass an allocatable array to procedure that isn't declared to use allocatable arrays, as long as the array is allocated before the call. (Of course, you can't use the array as an allocatable array in the procedure in which it is declared without that property.) Perhaps that will solve your problem. Allocate the array in the code that you write, than pass it as an argument to the FEM solver.

Example code: (I'd normally put the function into a module but you say that you can't do that, so I write an example showing the case of not using a module.)

function MySum ( RegArray )

real :: MySum
real, dimension (:), intent (in) :: RegArray

MySum = sum (RegArray)

end function MySum


program TestArray

   implicit none

   interface AFunc

      function MySum ( SomeArray )

         real :: MySum
         real, dimension (:), intent (in) :: SomeArray

      end function MySum

   end interface AFunc

   real, dimension (:), allocatable :: AllocArray
   integer :: N
   real :: answer

   write (*, '("Input array size: ")', advance="no")
   read (*, *) N

   allocate ( AllocArray (1:N) )
   AllocArray = 1.0

   answer = MySum ( AllocArray )
   write (*, *) answer

end program TestArray

---------- EDIT: Second Concept ---------

Sharing an allocatable array between two subroutines, without the calling routine being "aware" of the array.

module MySubs

   real, allocatable, dimension (:,:) :: array

contains


subroutine One ( x, y, ... N, M )

   integer, intent (in) :: N, M

   if ( .NOT. allocated (array) ) allocate ( array (N, M) )


end subroutine One


subroutine Two ( .... )


end subroutine Two


end module MySubs

UPDATE: note: This approach can be used to pass information between the two routines without the main program having access the module ... for the question, without modifying the original main prpgram. Part of the example is how to allocate the arrays: the example does that by having the subroutine that would first use the array test whether the array is allocated -- if not, it allocates the array.

Upvotes: 2

user1220978
user1220978

Reputation:

The three examples below all work with gfortran. The second may fail on some compilers as it uses a F2003 feature (allocatable dummy arguments), and not all compilers are 100% F2003 compliant. However, most implement ISO TR 15581 (which includes this feature).

First version, you can use a common pointer to allocatable array.

program hip
   implicit none
   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   double precision, allocatable, dimension(:, :), target :: a
   allocate(a(100, 100))
   a(1, 1) = 3.1416d0
   p => a
   call hop
   deallocate(a)
end program

subroutine hop
   implicit none
   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   print *, size(p, 1), size(p, 2), p(1, 1)
end subroutine

Second version, allocating in a subroutine then calling another. One still needs to declare the array in main program.

program hip
   implicit none

   interface
      subroutine hip_alloc(arr)
         double precision, allocatable, dimension(:, :) :: arr
      end subroutine
   end interface

   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   double precision, allocatable, dimension(:, :) :: a
   p => null()
   print *, "a:", allocated(a)
   print *, "p:", associated(p)
   call hip_alloc(a)
   print *, "a:", allocated(a)
   print *, "p:", associated(p)
   call hop
   deallocate(a)
end program

subroutine hip_alloc(arr)
   implicit none
   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   double precision, allocatable, dimension(:, :), target :: arr
   allocate(arr(100, 100))
   arr(1, 1) = 3.1416d0
   p => arr
end subroutine

subroutine hop
   implicit none
   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   print *, size(p, 1), size(p, 2), p(1, 1)
end subroutine

Third version, here we first call a function returning a pointer, then pass this pointer to a subroutine through a common. The function does the allocation, as in second example. The pointer is deallocated in main program, but could be elsewhere.

program hip
   implicit none

   interface
      function hip_alloc(n)
         integer :: n
         double precision, dimension(:, :), pointer :: hip_alloc
      end function
   end interface

   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   p => null()
   print *, "p:", associated(p)
   p => hip_alloc(100)
   print *, "p:", associated(p)
   call hop
   deallocate(p)
end program

function hip_alloc(n)
   implicit none
   integer :: n
   double precision, dimension(:, :), pointer :: hip_alloc
   allocate(hip_alloc(n, n))
   hip_alloc(1, 1) = 3.1416d0
end function

subroutine hop
   implicit none
   double precision, dimension(:, :), pointer :: p
   common /hiphop/ p
   print *, size(p, 1), size(p, 2), p(1, 1)
end subroutine

Upvotes: 2

Kyle Kanos
Kyle Kanos

Reputation: 3264

I do not understand why writing a MODULE would not work, but have you considered CONTAINS? Everything above the CONTAINS declaration is visible to the subroutines below the CONTAINS:

PROGRAM call_both
   INTEGER,DIMENSION(2) :: a, b
   a = 1
   b = 2
   PRINT *,"main sees", a, b
   CALL subA
   CALL subB
 CONTAINS
   SUBROUTINE subA
      PRINT *,"subA sees",a,b
   END SUBROUTINE subA

   SUBROUTINE subB
      PRINT *,"subB sees",a,b
   END SUBROUTINE subB
END PROGRAM call_both

The output would be

main sees           1           1           2           2
subA sees           1           1           2           2
subB sees           1           1           2           2

This works with ALLOCATABLE arrays as well.

Upvotes: 0

Related Questions