Array passed as argument changes dimensions?

1

I have attempted to convert a subroutine into a function (since I have now reduced the outputs to 1).

The original subroutine was this…

    subroutine girl(a, n, m)

		implicit none
		
		double precision, intent(inout) :: a(*)
		integer, intent(in)             :: m, n
		
		double precision                :: amax, temporaryStore
		integer                         :: i, i1, i2, ik, j, jk, jm, k

		save

		do i = 1, n
			amax = abs(a(i * n - n + i))
			jm   = i

			if ((i + 1) <= n) then
	  
				do j = (i + 1), n
					if (abs(a(i * n - n + j)) > amax) then
						amax = abs(a(i * n - n + j))
						jm = j
					end if
				end do
				
				if (jm /= i) then
					i1 = jm + i * n - n
					i2 = i * n - n + i
		
					do j = i, (n + m)
						temporaryStore = a(i1)
						a(i1)          = a(i2)
						a(i2)          = temporaryStore
						i1             = i1 + n
						i2             = i2 + n
					end do
				end if
			end if

			if (a(i * n - n + i) == 0d0) then
			
				write (6, '(a)') "error in subroutine girl"
				stop 'girl'
			
			end if

			do j = 1, n
			
				if (j /= i) then
					ik     = i * n - n + j
					jk     = i * n - n + i
		
					do k = (i + 1), (n + m)
						jk = jk + n
						ik = ik + n
						a(ik) = a(ik) + (-a(i * n - n + j) / a(i * n - n + i)) * a(jk)
					end do
				end if
			end do
			
			jk     = i * n - n + i

			do j = (i + 1), (n + m)
				jk    = jk + n
				a(jk) = a(jk) * (1d0 / a(i * n - n + i))
			end do
		end do
		return
	end subroutine girl

And the new function became this…

    function girl(a, n, m)

		implicit none
		
		double precision, intent(inout) :: a(:)
		integer, intent(in)             :: m, n
		double precision                :: girl(size(a))
		
		double precision                :: amax, temporaryStore
		integer                         :: i, i1, i2, ik, j, jk, jm, k

		save

		do i = 1, n
			amax = abs(a(i * n - n + i))
			jm   = i

			if ((i + 1) <= n) then
	  
				do j = (i + 1), n
					if (abs(a(i * n - n + j)) > amax) then
						amax = abs(a(i * n - n + j))
						jm = j
					end if
				end do
				
				if (jm /= i) then
					i1 = jm + i * n - n
					i2 = i * n - n + i
		
					do j = i, (n + m)
						girl(i1) = a(i2)
						girl(i2) = a(i1)
						i1       = i1 + n
						i2       = i2 + n
					end do
				end if
			end if

			if (a(i * n - n + i) == 0d0) then
			
				write (6, '(a)') "error in subroutine girl"
				stop 'girl'
			
			end if

			do j = 1, n
			
				if (j /= i) then
					ik     = i * n - n + j
					jk     = i * n - n + i
		
					do k = (i + 1), (n + m)
						jk = jk + n
						ik = ik + n
						girl(ik) = a(ik) + (-a(i * n - n + j) / a(i * n - n + i)) * a(jk)
					end do
				end if
			end do
			
			jk     = i * n - n + i

			do j = (i + 1), (n + m)
				jk    = jk + n
				girl(jk) = a(jk) * (1d0 / a(i * n - n + i))
			end do
		end do
		return
	end function girl

I get the following error…

stelcor.f90:3028:14:

 3028 |                         ha = girl (ha, mh, mh + 1)
      |                                   1
Error: Rank mismatch in argument 'a' at (1) (rank-1 and rank-2)

I’m assuming that this is due to the fact that the input is a 2 dimensional array yet appears in the original subroutine as a 1 dimensional array (which I have maintained in the new function).

Two things come to mind…

  • How has the 2-dimensional array become a 1-dimensional array; how has the system dealt with this without raising an error before?
  • And how do I resolve this error now that it is a function?
2

This works because of storage sequence association. To make it work like before, change the array declaration from (:) back to (*) as in the original. I notice that in the function, the array has intent(inout). It looks like it should now be intent(in).

Another option would be to change the dummy argument to a 2D array. Whether this is a good idea depends on why the 2D array addressing was changed to 1D in the first place. With legacy code, that is not always an easy answer to find.

3

I prefer the simple inclusion of the *. I was unsure what this was and had assumed it might just be another form of (:). Obviously I was wrong!

And I had forgotten to remove the out from the inout; thank you.

However, I am now meeting an error when setting the return value (where I am using size).

I now have…

    function girl(a, n, m)

		implicit none
		
		double precision, intent(in) :: a(*)
		integer, intent(in)          :: m, n
		double precision             :: girl(size(a))
		
		double precision             :: amax, temporaryStore
		integer                      :: i, i1, i2, ik, j, jk, jm, k

		save

And get the following…

stelcor.f90:3315:44:

 3315 |                 double precision             :: girl(size(a))
      |                                                          1
Error: The upper bound in the last dimension must appear in the reference to the assumed size array 'a' at (1)

Obviously I am not setting the dimension of the return name correctly; how should I be tackling this one?

4

The original subroutine modified the array in-place. It was effectively a 2D array with dimension (n,m). So the size of the output array in the function should be m*n. Note that the function will require memory allocation (probably on the stack) of that output array. This means that the function will be less efficient than the original subroutine. It also means that for large values of m and n, you may exhaust the stack and you will get run time errors. To fix those run time errors, you will need to change the stack size in the OS, or change the function output declaration from automatic to allocatable. Many fortran programmers like to avoid all of those possible problems by using the subroutine interface and modifying the array in-place.

5

I will keep the subroutine as-is. It works and if, as you state, the function would either be less efficient or require a total rewrite then that is not worth the effort - far too many other aspects to do.

Thank you so much for the help! :slight_smile:

6

Actually, I (might) lie.

Looking at it, mh is the variable defining the one value, and mh is set at 4. mh becomes n in girl, with m either being mh + 1 or simply 1.

Therefore, it shouldn’t take me too long to work out what the subroutine is doing, redefine the array as a 2-dimensional array and return it as-is. This should also allow me to lose some of the additional variable, saving space.

I am going to ponder this a little longer before attempting this, but if I end up with an efficient (time-wise) function that saves on space and makes the rest of the code easier to understand, then I win in all aspects.

Again, thank you; this does help enormously.

7

If your dimensions are 4, then you will never need to worry about stack space limitations.

If you do convert to a 2D dummy argument, then, with the way the original code is written, it should be straightforward. An expression like

a(i * n - n + j)

in the 1D case will become

a(j,i)

in the 2D case (assuming the array is dimensioned as a(n,*) or a(:,:)). The problem you might run into is if the routine is ever called with a 1D array actual argument. The (*) dummy argument allows it to be called with an actual argument of any rank (1D, 2D, 3D, etc) through storage association. But changing it to a 2D array, with an explicit interface (which is a good idea for many other reasons), will limit the actual arguments to be just 2D arrays. There are work-arounds for that, involving for example pointers, but it is something that you might encounter with such a change.

With new code, you mostly never have to worry about changing ranks. But legacy code (f77 and earlier) used to do this routinely as a way to reuse memory in different ways. Before f77, the common convention was to dimension the array as (1), and then violate that upper bound freely. All of these complications regarding rank matching with different types of dummy argument declarations involve maintaining backwards compatibility with those kinds of legacy codes.

8

All arguments are 2D arrays; it is called 4 times in total, with the dimensions varying but still no bigger than 5 (twice (4, 5), once (4, 1) and once (3, 1)).

I will be re-writing this. It helps tidy the code up and allows me to work out exactly what is happening. There is some swapping going on in here, so it appears to be a sorting algorithm. With 2 loops I would assume a form of bubble sorting but that’s with a cursory glance.

Once again, thank you. The code is so much better than when I started; in large part to your help and the help of others.