Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7174 Discussions

Broken Fortran bindings for Nonlinear least squares solver routines using RCI.

Stijn_Schildermans
83 Views

I tried to use the dtrnlspbc_* and djacobi_* routines from Fortran. I came across several issues during my attempts:

  • The documentation for dtrnlspbc_* routines states that handle is an integer*8 while according to the include file (version 2025.0.1, but also for earlier versions) it is declared as a type(HANDLE_TR) which contains 2 additional fields m and n, which I assume are the problem dimensions.
  • When compiling with gfortran (optimization level 2 or 3) some random illegal memory access occurs within djacobi_solve and dtrnlspbc_solve, leading to undefined behavior. I was not easily able to create a minimal example reproducing this issue (undefined behavior is a lot of fun...), but I thoroughly checked to make sure it was not caused by my code. Code compiled with IFX does not suffer from this problem, nor does code compiled with gfortran using only the Djacobi routine (without RCI). I therefore suspect the issue has to do with how gfortran interprets the RCI handle.
  • I used GDB to check the contents of the handle during execution and it seems that the m and/or n fields of the handle contain a pointer to the memory region used internally by MKL while the actual handle field remains uninitialized. I guess this in itself can be OK, but mkl_rci.fi declares dtrnlspbc_solve as follows: 
     INTERFACE
        INTEGER FUNCTION DTRNLSPBC_SOLVE(handle, fvec, fjac, rci_request)
        USE MKL_RCI_TYPE
        INTEGER, INTENT(OUT)   :: rci_request
        DOUBLE PRECISION, INTENT(INOUT), DIMENSION(*):: fvec
        TYPE(HANDLE_TR), INTENT(IN):: handle
        DOUBLE PRECISION, INTENT(IN), DIMENSION(handle%m, *):: fjac
        END FUNCTION
      END INTERFACE

      Therefore, if handle%m is a random value (which may be 0), the entire fjac parameter may be optimized out, which might perhaps cause the problems I am seeing? In any case, fjac should not be declared  as dimension(handle%m,*) if handle%m does not actually contain m at runtime.

 

I wrote a C wrapper for the nonlinear solver I built for my project so that I can call the MKL procedures I need from C and as such bypass the MKL Fortran bindings. This works without problems for me. In any case, I thought it would be useful to report my adventure so it can be fixed or someone can point out to me exactly in which way I am mistaken.

0 Kudos
0 Replies
Reply