Swapping matrix names in a loop

1

Hello. I have two matrices A and B that I am taking through a loop and each iteration I would like to swap their places (like how you solve Laplace’s equation by relaxation).

Now I know you could have an even and odd loop (which is what I did before, but now I am changing that part of my code) and essentially duplicate code (or write a function and call it changing up matrix arguments), but for what I am doing I do not want to do that.

My arrays happen to be allocatable so I came up with two solutions (besides plain copying): make them pointers (or targets) and interchange what they point at, or use move_alloc in basically the same way. Here is an example of each (obviously not inside the loop I am talking about):

program main
    use iso_fortran_env, only: dp => real64
    implicit none

    integer, parameter :: m = 4
    real(dp), pointer :: A(:, :), B(:, :), temp(:, :)

    allocate(A(m, m), B(m, m))
    call random_number(A)

    ! A has random numbers I want to refer to by the name B,
    ! and I want to use the name A again for something else
    temp => A
    A    => B
    B    => temp
    ! Now B has the random numbers A had,
    ! and I can use A for something else.

    deallocate(A, B)
endprogram main

program main
    use iso_fortran_env, only: dp => real64
    implicit none

    integer, parameter :: m = 4
    real(dp), allocatable :: A(:, :), B(:, :), temp(:, :)

    allocate(A(m, m), B(m, m))
    call random_number(A)

    ! A has random numbers I want to refer to by the name B,
    ! and I want to use the name A again for something else
    call move_alloc(A, temp)
    call move_alloc(B, A)
    call move_alloc(temp, B)

    ! Now B has the random numbers A had,
    ! and I can use A for something else.

    deallocate(A, B)
endprogram main

Ideally I would like to not use pointers since they are a pain in Fortran and I have heard they can be quite bad on performance especially in Fortran. move_alloc is a thing I have only ever used once or twice before in the language for minor tasks, so I am uncertain about its performance impact (it could so happen it secretly will just keep allocating and deallocating in my code and be even worse than plain copying).

Anyways, I thought this sort of thing would be a somewhat common issue so there would be a “figured out” solution in the language but I haven’t been able to find one. If anyone has a good solution (or that one doesn’t exist, so I would have to suck it up and do the loop alternating thing) for this I would appreciate it (now that I type this I think clever use of associate might be another possible solution, but I’ll have to try it out).

1 Like
2

I’ve used even and odd loops like this:

if (odd) then
    call work(A, B)
else
    call work(B, A)
end if

So the swapping is explicit (and free) and the work function then has the correct order.

I also don’t like the pointers, like you said.

1 Like
3

The move_alloc() approach will not allocate new memory. You can check this by using c_loc(). You will see that it is the same memory just switching back and forth, which is what you want to happen.

4

It’s an old thread, you may have missed it, but there was a nice discussion here:

https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping

Can you explain how to do this please?

Every Fortran programmer :sweat_smile:

1 Like
5

You can think of the move_alloc intrinsic as simply doing some internal renaming. The data itself remains unaffected.

(move_alloc is one of the few pure intrinsic subroutines —since most of the pure intrinsics are functions).

Something like this shows that the descriptors are still at the same memory addresses:

use ISO_C_BINDING

implicit none

character(*), parameter :: FMT = '(a,i0)'
integer :: i
real, allocatable, target :: a(:), b(:), tmp(:)

a = [real :: (i, i = 1, 10)]
b = [real :: (i ** 2, i = 1, 10)]

do i = 1, 3
    print FMT, 'i=', i
    print '(a)', achar(9)//'descriptors before swap'
    print FMT, achar(9)//'a=', transfer(c_loc(a), 1_c_intptr_t)
    print FMT, achar(9)//'b=', transfer(c_loc(b), 1_c_intptr_t)
    call move_alloc(a, tmp)
    call move_alloc(b, a)
    call move_alloc(tmp, b)
    print '(a)', achar(9)//'descriptors after swap'
    print FMT, achar(9)//'a=', transfer(c_loc(a), 1_c_intptr_t)
    print FMT, achar(9)//'b=', transfer(c_loc(b), 1_c_intptr_t)
enddo
end

Which outputs something like

$ ifort -diag-disable=10448 -stand=f23 -standard-semantics swap_target.f90 && ./a.out 
i=1
	descriptors before swap
	a=157749984
	b=157750080
	descriptors after swap
	a=157750080
	b=157749984
i=2
	descriptors before swap
	a=157750080
	b=157749984
	descriptors after swap
	a=157749984
	b=157750080
i=3
	descriptors before swap
	a=157749984
	b=157750080
	descriptors after swap
	a=157750080
	b=157749984

The output is similar with gfortran (v14) and flang (v21).

1 Like
6

I can’t believe I missed that thread when searching up for what to do. It would have saved some headache for me last night. Better late than never, thanks.

7

The move_alloc() intrinsic does what is called a “shallow copy” operation. The underlying data is untouched, but the names of the arrays are switched. An important practical consequence of that is that the effort (that is, the cpu cycles) required for the shallow copy operation does not depend on the size of the array (which could be GBs or larger). For an array of size N, the operation requires O(N**0) effort instead of the O(N**1) effort that would be required by an actual copy operation. Some of the other suggestions in this thread have that same feature, such as toggling the index of an array, swapping pointer assignments, and swapping two actual arguments of a subroutine. The move_alloc() intrinsic is just one more way to achieve that same kind of efficiency.

1 Like
8

You meant “shallow move”? The allocation is “moved”, not copied —the FROM argument becomes unallocated afterwards, unlike what would happen in a copy (either shallow or deep).

(My “renaming” analogy comes from the fact that in UNIX the mv command is the one used for renaming.)

9

Yes, the from= argument in move_alloc() is deallocated during the operation, so that there is always only one of the two arguments associated with the memory, before and after. The move_alloc() subroutine is appropriately named. I’ve just never heard the term “shallow move” before.

As for the original question of whether new memory might be allocated during the operation, another important feature of move_alloc() is that any pointers that were originally associated with the from= argument before the statement become associated with the to= argument after the statement. Those pointers do not become undefined. That would be inconceivably difficult to achieve if new memory were allocated during the operation, while it is trivial to achieve if it is the same memory that is just transferred to the new name.

1 Like
10

You’re right. Borrowing from C++'s move semantics, that’ll just be a “transfer of ownership”.