gfortran SIGABRT error: the zgesvd function from lapack

Question

I am coding a fortran 90 subroutine that calls the lapack function zgesvd. Unfortunately I get this error:

    malloc(): corrupted top size
    Program received signal SIGABRT: Process abort signal.
    Backtrace for this error:
    #0  0x7f338321cd21 in ???
    #1  0x7f338321bef5 in ???
    #2  0x7f3382f0120f in ???
    #3  0x7f3382f0118b in ???
    #4  0x7f3382ee0858 in ???
    #5  0x7f3382f4b3ed in ???
    #6  0x7f3382f5347b in ???
    #7  0x7f3382f56839 in ???
    #8  0x7f3382f58418 in ???
    #9  0x560e60a94daf in ???
    #10  0x560e60a94b75 in ???
    #11  0x560e60a94bae in ???
    #12  0x7f3382ee20b2 in ???
    #13  0x560e60a9420d in ???
    #14  0xffffffffffffffff in ???
    Aborted (core dumped)

The minimal code containing the subroutine is written below. Since Lapack calls are rather elaborated, it has to be a bit lengthy to contain the details about every argument, I apologize for that:

module my_mod
    IMPLICIT NONE

    CONTAINS
    subroutine svd_gc(A,u,s,vt,info)
        integer, parameter :: pr = SELECTED_REAL_KIND(15,3)
        complex(pr), dimension (:,:), intent(in)  :: A
        integer :: M, N
        
        integer :: lda, ldu, ldvt

        real(pr), dimension (min(SIZE(A,1),SIZE(A,2))) :: s
        complex(pr), dimension (SIZE(A,1),SIZE(A,1)) :: u
        complex(pr), dimension (SIZE(A,2),SIZE(A,2)) :: vt
        integer lwork
        complex(pr), dimension (1) :: work1
        real, dimension (5*min(SIZE(A,1),SIZE(A,2))) :: rwork
        integer :: info
        complex(pr), allocatable, dimension(:) :: work2
        
        external zgesvd

        M = SIZE(A,1)
        N = SIZE(A,2)
        lda = M
        ldu = M
        ldvt = N

        ! initialization step
        call zgesvd('All', 'All', M, N, A, lda, s, u, ldu, vt, ldvt, work1, -1, rwork, info )
        lwork = min(1000, int( work1(1) ))
        allocate( work2(lwork) )

        call zgesvd('All', 'All', M, N, A, lda, s, u, ldu, vt, ldvt, work2, lwork, rwork, info )

    end subroutine svd_gc
end module my_mod

which is called here:

program main
    use my_mod

    implicit none
    
    integer, parameter :: pr = SELECTED_REAL_KIND(15,3)
    real(pr), dimension (5,3) :: MR, MI
    complex(pr), dimension (5,3) :: M

    complex(pr), dimension(5,5) :: u
    complex(pr), dimension(3,3) :: vt
    real(pr), dimension(3) :: s
    integer i, j, info
    
    call random_vec(5,3,MR) ! generates an array of real numbers
    call random_vec(5,3,MI)
    
    M = cmplx(MR,MI)
    
    call svd_gc(M,u,s,vt,info)

    do j=1,3
        print '("---------")'
        print*,s(j)
        
        do i=1,3
            print*,vt(i,j)
        end do
        
        print '(" ")'
    end do

end program main

and I compile with gfortran... to finally run the produced file

Unfortunately this call returns a SIGABRT error, so I found a working example of zgesvd call and I took that one from Intel website: https://www.intel.com/content/www/us/en/develop/documentation/onemkl-lapack-examples/top/least-squares-and-eigenvalue-problems/singular-value-decomposition/gesvd-function/zgesvd-example/zgesvd-example-fortran.html

basically you dont need work2, you set the size of work1 to 1000, and you call:

    call zgesvd('All', 'All', M, N, A, lda, s, u, ldu, vt, ldvt, work1, -1, rwork, info )
    lwork = min(1000, int( work1(1) ))

    call zgesvd('All', 'All', M, N, A, lda, s, u, ldu, vt, ldvt, work1, lwork, rwork, info )

ok it works thank you Intel!

HOWEVER now if I dare using the 'allocate' command in any code that calls my zgesvd subroutine (svd_gc), then I get again a SIGABRT error. I am wondering, I am doing anything wrong ? or is the zgesvd function still not stable, maybe with some 'free()' calls misplaced or ill defined ? Do you see any correction or workaround ? maybe another Fortran library ?

I saw a very similar thread: gfortran error: zgesvd in lapack however, I am not using this kind of 'USV' type definition so I don't think it is related

Thanks in advance! Jee

edit: In this edit, I did change the "double" assignement 'REAL*8' to another one which I hope is more clean. Also, I added the minimal code for the main program that calls the subroutine !

Answer 1

Ok I got the error thanks to Ian Bush comment about the 'rwork' variable. 'rwork' has to be double precision as it is stated in the prototype here: http://www.netlib.org/lapack/explore-html/d3/da8/group__complex16_g_esing_gad6f0c85f3cca2968e1ef901d2b6014ee.html

Thus, replacing: real, dimension (5*min(SIZE(A,1),SIZE(A,2))) :: rwork
by: real(pr), dimension (5*min(SIZE(A,1),SIZE(A,2))) :: rwork
seems to solve the issue. Thank you for constructive comments and critics, it helped to improve !