Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Sourced allocation from transposed array

I'm trying to understand a bug that I encountered in my code. When trying to allocate an array using the transpose of an array, I use an allocate statement with source=transpose(original_array) in my code. However, I don't get the expected result using this method. It seems that the indexing is off by one and the first row of the source array is skipped.

Example:

program testalloc
    real*8, allocatable :: a(:, :)
    real*8, allocatable :: b(:, :)

    allocate(b(2, 3))
    b(1, :) = [1, 2, 3]
    b(2, :) = [4, 5, 6]
    call printmat(b)

    a = transpose(b)
    call printmat(a) ! Good

    deallocate(a)
    allocate(a(3, 2), source=transpose(b))
    call printmat(a) ! Bad

    deallocate(a)
    allocate(a(3, 2))
    a = transpose(b) 
    call printmat(a) ! Good

contains

    subroutine printmat(mat)
        real*8, intent(in) :: mat(:, :)
        integer :: i

        write(*,*) 'print'
        do i = 1, size(mat, 1)
            write(*,*) mat(i, :)
        end do
    end subroutine

end program

which gives

 print
   5.0000000000000000        3.0000000000000000     
   6.0000000000000000        0.0000000000000000     
   3.2114266979681025E-322   5.0000000000000000  

for the sourced allocation after compiling with gfortran (gcc version 7.3.0 (Ubuntu 7.3.0-27ubuntu1~18.04)) instead of the transposed original array. Am I doing something wrong here or is it a compiler bug?

like image 318
MSap Avatar asked Oct 27 '22 14:10

MSap


1 Answers

The sourced allocation

allocate(a(3, 2), source=transpose(b))

is a valid one. The shape specified for a is the same as that of the source expression transpose(b). As a result, a takes the value of the expression given.

It is incorrect for the compiler to give a different result. The output routine is not to blame here.

gfortran 8 appears to give the expected output.

Interestingly, with gfortran 7 the expected result appears if b is not itself allocatable.

like image 123
francescalus Avatar answered Nov 09 '22 16:11

francescalus