pzlatra segfault with GCC

[JamesGreen29]

Hello all,

I've noticed that the ScaLAPACK complex trace, pzlatra, fails with a segfault when using GCC, but it does not with intel compilers. Please see below for a minimal example, using a 2x2 BLACS grid:

program pzlatra_example
    implicit none

    ! ScaLAPACK functions
    complex*16  :: pzlatra
    external    :: pzlatra, blacs_pinfo, blacs_get, blacs_gridinit, descinit

    ! Global Matrix dimensions
    integer, parameter :: N = 4
    integer, parameter :: NB = 2  ! Block size

    ! Variables
    integer :: mypnum, nprocs, ictxt, myrow, mycol
    integer :: nprow, npcol, info
    integer :: descA(9)
    complex*16 :: A_local(NB, NB)  ! Local portion of the matrix
    complex*16 :: trace_result

    ! 1. Initialize BLACS
    call blacs_pinfo(mypnum, nprocs)
    nprow = 2 ! 2x2 process grid
    npcol = 2
    call blacs_get(-1, 0, ictxt)
    call blacs_gridinit(ictxt, 'Row-major', nprow, npcol)
    call blacs_pcoord(ictxt, mypnum, myrow, mycol)

    ! 2. Initialize the Matrix Descriptor
    ! descinit(descriptor, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
    call descinit(descA, N, N, NB, NB, 0, 0, ictxt, NB, info)

    ! 3. Fill local matrix with dummy data
    A_local = (1.0d0, 0.5d0) 

    ! 4. Call PZLATRA
    ! Syntax: PZLATRA(n, A, ia, ja, desca)
    if (mypnum == 0) print *, "Calculating Trace..."
    
    trace_result = pzlatra(N, A_local, 1, 1, descA)

    ! 5. Root process prints result
    if (myrow == 0 .and. mycol == 0) then
        print *, "The trace of the distributed matrix is: ", trace_result
    end if

    call blacs_gridexit(ictxt)
    call blacs_exit(0)
end program

Here is my compilation command for GCC (version 13.3.0)

mpif90 -o pzlatra_example.x pzlatra_example.f90 \
-L${MKLROOT}/lib/intel64 \
-lmkl_scalapack_lp64 \
-lmkl_blacs_intelmpi_lp64 \
-lmkl_gf_lp64 \
-lmkl_sequential \
-lmkl_core

 And here is my compilation command for ifx (oneapi 2025.3)

mpiifx -o pzlatra_example.x pzlatra_example.f90 \
-L${MKLROOT}/lib/intel64 \
-lmkl_intel_lp64 \
-lmkl_sequential \
-lmkl_core \
-lmkl_blacs_intelmpi_lp64 \
-lmkl_scalapack_lp64

 I notice that in the mkl_scalapack.h file, pzlatra is defined as follows:

void  pzlatra(MKL_Complex16 *, const MKL_INT *n, const MKL_Complex16 *a, const MKL_INT *ia, const MKL_INT *ja, const MKL_INT *desca);

 This looks different to the pdlatra definition:

double  pdlatra(const MKL_INT* n, const double* a, const MKL_INT* ia, const MKL_INT* ja, const MKL_INT* desca);

And this is how pzlatra is defined in the reference Fortran implementation:

COMPLEX*16         FUNCTION PZLATRA( N, A, IA, JA, DESCA )

It looks like pzlatra is defined more like a subroutine than a function in the MKL header, and it seems to me that the result of pzlatra is being sent back to the memory address given by the first argument, which causes the segfault for GCC. Is my interpretation correct? If so, what is the reason for MKL doing this, and how can I get the function to work as expected with GCC?

Thank you!