arraysfortranlapackcontiguousassociate

Is it safe to pass a pointer of a non-contiguous array to a Lapack95 routine in Fortran?


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?


Solution

  • 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":

    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 )