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!

[sc_Intel]

Thank you for reporting this issue. I was able to reproduce the behavior on my side as well. During the investigation, I observed an inconsistency between the mkl_scalapack.h interface and the oneMKL Fortran developer reference for this function.

I’ve logged the issue internally so the team can review it in detail. I will update you as soon as there is any progress or additional information to share.

Thanks again for bringing this to our attention.

[sc_Intel]

Hello @JamesGreen29, thank you again for reporting this issue. The oneMKL ScaLAPACK header file, mkl_scalapack.h, has been updated in oneMKL version 2026.0.0, which you can download from the link Get Intel® oneAPI Math Kernel Library (oneMKL) .

[JamesGreen29]

Thank you very much for the update! I see the header file got updated, but unfortunately this doesn't seem to resolve the segfault problem for PZLATRA with GCC when using MKL 2026.0 + IMPI 2021.18 from OneAPI 2026.0.

[sc_Intel]

Hello @JamesGreen29 , The segfault is an ABI mismatch on COMPLEX*16 function returns: Intel Fortran (ifx) and gfortran use different conventions for returning a complex result. oneMKL's pzlatra is built with the Intel Fortran convention, so the `val = pzlatra(...)` form shown in the Fortran developer reference works under ifx but produces the segfault you observed when called from gfortran. For gfortran, please call PZLATRA in subroutine form, with the result passed as the first argument: call pzlatra(trace_result, N, A_local, 1, 1, descA)

Please let us know if this resolves the issue on your side.

[JamesGreen29]

Hi @sc_Intel , thanks for the reply and explanation! I can confirm that the subroutine approach works for gfortran (with both MKL 2026.0 and 2025.3)