Subroutine not returning correct numerical values in assumed shape array due to index renaming in the main program

Question

The argument of my fortran 95 subroutine is an assumed shape array with intent inout:

the_subroutine(my_argument)
real, dimension(:,:), intent(inout) :: my_argument
(...)

In the main program, I have an allocatable array. I allocate it and also rename indexes. Then I call the subroutine and pass that (correctly allocated) array to the subroutine:

allocate(the_array( 5:1005 , 5:1005 ))
call the_subroutine(my_argument=the_array)

The subroutine does certain calculations and fills the array with values. In the very last line before the end of the subroutine, I check a random value:

(...)
print*, my_argument(213,126) ! I get 2.873...
end subroutine the_subroutine

Then, in the very first line after the subroutine call, I check if the value has been correctly communicated by the subroutine to the outer world, but that is not the case:

call the_subroutine(my_argument=the_array)
print*, the_array(213,126) ! I get 3.798... A completely different value.

The problem arises from having re-indexed the array in the main program as:

allocate(the_array( 5:1005 , 5:1005 ))

where max_index - min_index = 1000-1, but the subroutine "sees" the array internally as if I had declared the normal way, i.e.:

allocate(the_array( 1:1000, 1:1000))

Or simply, allocate(the_array( 1000, 1000 ))

Therefore, the element (213,126) in the internal array is in another location as in the main program array. Is there any easy way out of this?

my_argument is assumed shape. A procedure with assumed shape arguments requires an explicit interface in the calling scope. Does the_subroutine have such an interface? — IanH
– IanH, Commented Dec 7, 2016 at 7:53
@IanH It is inside a module, so it does nor require any interface. I think. — Mephisto
– Mephisto, Commented Dec 7, 2016 at 21:12
Do you have INTENT(INOUT) on the subroutine? (You may need to show more code) — Holmz
– Holmz, Commented Dec 7, 2016 at 21:16
If it is ALLOCATABLE in the main, then should it also be ALLOCATABLE in the subroutine? — Holmz
– Holmz, Commented Dec 7, 2016 at 21:19

IanH · Accepted Answer · 2016-12-07 22:05:37Z

2

The default lower bound for an assumed shape array is one.

If you want a different lower bound, then you need to declare the array dummy argument appropriately.

subroutine the_subroutine(my_argument)
  real, dimension(5:,5:), intent(inout) :: my_argument
!                 ^  ^

(The rules for the bounds of the dummy argument are different for deferred shape arrays - array dummy arguments which also have the POINTER or ALLOCATABLE attributes.)

answered Dec 7, 2016 at 22:05

IanH

21.6k2 gold badges41 silver badges61 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Mephisto Over a year ago

Thanks. However, it is of not much value using an assumed shape array, if I have to know in advance the LBOUND value.

agentp · Accepted Answer · 2016-12-08 02:45:59Z

2

use lbound to pass the bounds to the subroutine:

  implicit none
  real,allocatable:: x(:,:)
  allocate(x(5:10,5:10))
  call sub(x,lbound(x))
  write(*,*)'o',x(5,5)
  contains

  subroutine sub(x,lb)
  implicit none
  integer lb(2)
  real, dimension(lb(1):,lb(2):)::x
  x(5,5)=42.
  end subroutine
  end

o 42.0000

edited Dec 8, 2016 at 2:45

answered Dec 7, 2016 at 22:16

agentp

7,0072 gold badges22 silver badges39 bronze badges

1 Comment

Mephisto Over a year ago

Upvoted because that is a perfectly good solution. However if you pass several arguments and several lbound values it can be very messy.

Mephisto · Accepted Answer · 2016-12-08 05:25:28Z

0

Finally, I found the solution.

First, if working in Fortran 2003 (or Fortran 95 with non-standard extensions), you may simply declare the assumed shape argument in the subroutine as ALLOCATABLE:

subroutine the_subroutine(my_argument)
real, dimension(:,:), allocatable, intent(inout) :: my_argument

Then the subroutine "sees" the renamed index correctly. However this is not allowed in the Fortran 95 standard.

In Fortran 95, the most elegant way I found for this is by making use of a pointer:

program example
implicit none
real, dimension(:,:), allocatable, target  :: the_array
real, dimension(:,:),              pointer :: the_pointer
[...]
allocate(the_array(5:1005,5:1005))
the_pointer => the_array
call the_subroutine(my_argument=the_pointer)

And in the subroutine:

subroutine the_subroutine(my_argument)
real, dimension(:,:), pointer :: my_argument

Then it works perfectly. Inside the subroutine, MY_ARGUMENT is treated exactly as if it was an assumed shape array.

edited Dec 8, 2016 at 5:25

answered Dec 8, 2016 at 5:19

Mephisto

7001 gold badge7 silver badges12 bronze badges

2 Comments

IanH Over a year ago

When you add the pointer or allocatable attribute you get what is known as a deferred shape array. By adding such an attribute you change the utility of the procedure - for example you are no longer able to pass array sections. Adding pointer is also a reasonable way of telling the compilers optimiser that it doesn't have to do any work today.

Mephisto Over a year ago

@IanH Could you elaborate a bit more on "Adding pointer is also a reasonable way of telling the compilers optimiser that it doesn't have to do any work today."? - Also, I don't agree that you are no longer able to pass array sections. You can pass array sections by doing the same trick, just pointing to the array section, e.g.: the_pointer => the_array( : , 10:20 )

agentp · Accepted Answer · 2016-12-08 15:18:08Z

0

another approach to the problem: make the array a derived type and the subroutine a method of the type:

  module example
  implicit none
  type a
  real, allocatable::y(:,:)
  end type a
  contains
  subroutine sub(this)
  type(a)::this
  write(*,*)lbound(this%y)
  end subroutine
  end module

  program p
  use example
  implicit none
  type (a)  my_array
  allocate(my_array%y(5:6,7:8))
  call sub(my_array)
  deallocate(my_array%y)
  allocate(my_array%y(2:3,1:2))
  call sub(my_array)
  end

5 7

2 1

as you see now the subroutine automatically knows the correct dimensions. Obviously now the subroutine can only work with the derived type array.

answered Dec 8, 2016 at 15:18

agentp

7,0072 gold badges22 silver badges39 bronze badges

Collectives™ on Stack Overflow

Subroutine not returning correct numerical values in assumed shape array due to index renaming in the main program

4 Answers 4

1 Comment

1 Comment

2 Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

4 Answers 4

1 Comment

1 Comment

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related