I would like to call the syevd
Lapack95 subroutine for a non-contigous array as follows:
real :: mat(15000, 15000), vec(15000)
mat=1.d0
associate(eig_vects=>mat(:10000,:10000),eig_vals=>vec(:10000))
call syevd(eig_vects,eig_evals,'V')
end associate
Is it a safe thing to do this?
Earlier I experienced problems with passing non-contigous array slices to a subroutine (the temporary array created by that subroutine could not be allocated).
May I expect similar issues when using the eig_vects
pointer as an argument?
Lapack95 is not an independant implementation of Lapack, but rather a wrapper to the the classical Lapack, with a more modern interface. In particular, the dummy array arguments of Lapack95 are assumed-shape (which can be discontiguous), whereas they are assumed-size in classical Lapack (hence contiguous).
Since the end of the call chain (classical Lapack) expects contiguous arrays, a temporary copy has to be made anyway if you want to process discontiguous arrays. The only question is "where":
associate( eig_vects=>mat(:10000,:10000) )
doesn't make any copy, eigs_vects
is just a kind of alias herecall syevd(eig_vects,eig_evals,'V')
doesn't make any copy either, as the interface of the Lapack95 syevd
can handle a discontiguous array.*syevd
; but you have no control on how it is performed. For instance, typically the Intel compiler uses the stack for temporary arrays, and it will fails with such a large matrix.As suggested by @VladimirFГероямслава , you'd better explicitly allocate a contiguous array to work with once for all:
real, allocatable :: eigs_vect(:,:)
allocate( eigs_vect, source=mat(:10000,:10000) )
Note that vec(:10000)
is contiguous, so you don't need allocating another array for it.
Note that you could also directly call the classical Lapack routines, as their interfaces can handle this kind of situation without any copy, thanks the decoupling because the size of the problem and of the "leading dimension":
call ssyevd( 'V', 'U', 10000, mat, 15000, vec &
, work, lwork, iwork, liwork, stat )