PARDISO works well with IFORT but fails with error -4 in phase 12 when using with IFX

[GuoYJ]

I am using MKL PARDISO for solving linear equations where the system matrix is real, symmetric, postive definite (with mtype = 2). The codes can be compiled with both IFORT and IFX. However, if IFX (tried versions 2022, 2023, 2024.0, 2024.2) was used, the code always gave me an error code -4 in the phase 12 (actually the phase 11). The fact is that if the code was compiled with IFORT, it worked correctly for the same test cases. Is there any special setting for using PARDISO with IFX? Thanks a lot! Below is my code for calling PARDISO.

BTW, I prefer to use IFX because it seems that the time for preparing the equations (including reading file from the hard disk) is more stable when the code is compiled with IFX. This is important when lots of repeating work is required.

Code for calling PARDISO (on Windows 11):

include 'mkl_pardiso.f90'
subroutine solver_mkl(a,ja,ia,nnz,b,x,n,nrhs)
! Solving linear A*x = B using Intel MKL pardiso direct sparse solver.
! For A, CSR3 format is adopted; Assuming symmetric & positive definite --> upper triangle is used.

use mkl_pardiso
use file_units

implicit none
integer, parameter :: dp = kind(1.D0)
integer, intent(in) :: n ! Number of rows of matrix A = number of equations to solve
integer, intent(in) :: nnz ! Number of non-zero elements in A
integer, intent(in) :: nrhs ! Number of right-hand sides that need to be solved for
real(dp), intent(in) :: a(nnz) ! Values of non-zero elements in A
integer, intent(in) :: ja(nnz) ! Column indices of non-zero elements in A
integer, intent(in) :: ia(n+1) ! Row index, pointing to first column index of each row in ja
real(dp), intent(inout) :: b(n,nrhs)! Single RHS, B
real(dp), intent(out) :: x(n,nrhs) ! Solution of A*x = B

type(MKL_PARDISO_HANDLE) :: pt(64) ! Internal solver memory pointer
integer :: iparm(64) ! PARDISO control parameters
integer :: maxfct, mnum, mtype, phase, msglvl, error1, error2, error3 ! All other variables
integer :: idum(1) ! Dummy variable
real(dp) :: ddum(1) ! Dummy variable
integer :: i ! Loop index

maxfct = 1
mnum = 1
mtype = 2 ! Matrix A: real, symmetric, postive definite -- upper triangle
msglvl = 0 ! Print statistical information
error1 = 0; error2 = 0; error3 = 0 ! Initialize error flag

!! PARDISO control parameters
iparm = 0
iparm(1) = 1 ! No solver default
iparm(2) = 2 ! Fill-in reordering from METIS
iparm(3) = 4 ! Number of threads
iparm(4) = 0 ! No iterative-direct algorithm
iparm(5) = 0 ! No user fill-in reducing permutation
iparm(6) = 0 ! Solution -> x
iparm(8) = 2 ! Numbers of iterative refinement steps
iparm(10) = 13 ! Perturb the pivot elements with 1E-13
iparm(11) = 1 ! Use nonsymmetric permutation and scaling MPS
iparm(13) = 0 ! Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
iparm(27) = 0 ! 输入参数检查器开关

do i = 1, 64
pt(i)%DUMMY = 0
end do

phase = 12 ! Reordering, symbolic factorization and numerical factorization
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
idum, nrhs, iparm, msglvl, ddum, ddum, error1)
!write(unit_out0,*) 'Space requred 11: ', iparm(17), 'KB'

if(error1 == 0) then
    !write(*,*) 'Number of nonzeros in factors = ', iparm(18)
    !write(*,*) 'Number of factorization MFLOPS = ', iparm(19)

    phase = 33 ! Solving equations
    CALL pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
    idum, nrhs, iparm, msglvl, b, x, error2)

    if(error2 /= 0) then
        write(unit_out0,*) ' ** ERROR was detected in phase 33 of the PARDISO solver: ', error2
        stop 2
    end if

!write(unit_out0,*) 'Space requred 33: ', iparm(17), 'KB'

else
    write(unit_out0,*) ' ** ERROR was detected in phase 12 of the PARDISO solver: ', error1
    write(unit_out0,*) ' Note: Please check if there are rigid body displacements. Please make sure      there are sufficient &
    &displacement boundary conditions.'
    stop 1
end if

phase = -1 ! Releasing internal memory
call pardiso(pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
idum, nrhs, iparm, msglvl, ddum, ddum, error3)

if(error3 /= 0) then
    write(unit_out0,*) ' ** ERROR was detected in phase -1 of the PARDISO solver: ', error3
    stop 3
end if

end subroutine solver_mkl

[Fengrui]

Using ifort or ifx shouldn't affect the call to oneMKL Pardiso. I would recommend to check the data preparation code first, to see if the input data for Pardiso is the same using ifort and ifx.

If that is not the case, could you please share the data (ia, ja, a) used in this case? We would need to investigate further on our side.

Thanks!

[Fengrui]

Hi,

Did you get a chance to check the input data for the solver?

Thanks,

Fengrui

[GuoYJ]

Dear Fengrui,

Thanks for the guide! Checked an example and found that the input data was really different -- then the issue should relate to IFX itself, NOT the PARDISO solver. I did think normal operations should be same for IFX and IFORT, thus suspected that something might be wrong with the solver.

The finding of the check is: IFX and IFORT have different treatments on the string length when using a string type dummy parameter in calling a subroutine in which the string parameter is intent(out). With IFORT, the length of the string will be truncated to the defined length in the caller when it is returned from the subroutine in case its length defined in the subroutine is longer. However, with IFX, the length will changed to the length defined in the subroutine thus if the length defined in the caller is smaller, the string will occupy the space of adjacent variables, which will lead to error!

Thanks a lot!

Yongjin

[Fengrui]

It's glad to hear that the issue has been resolved!

If the different behavior of ifort and ifx is a concern, feel free to post it in the Fortran compiler forum. Thanks!

[morskaya_svinka_1]

In the current 2025.0 documentation iparm(3) is reserved and recommended to be set to zero, but in your code it is set to 4 and commented as a number of threads. Are you using an old version of MKL where this parameter is used? The second thing I would do is initialize all the iparm values explicitly as zero because I'm bad at fortran and do not know what compilers put there by default. Then I would look through the documentation and set the default iparms for SPD matrices except for those you want to change. For example, iparm(11)=1 should not be good for SPD problems. May I ask you for what purposes do you use fortran interface and not C/C++?

[GuoYJ]

Thanks a lot for the suggestions regarding version 2025. The code is for finite element analysis. May change to C++ later. Fortran was the first programming language I learned in school back to 30+ years ago. Just began to use C++ in recent two years.

[morskaya_svinka_1]

I do not mean I consider using C/C++ better, but I'm interested whether people use Fortran nowadays just because of plenty scientific legacy codes or there are other reasons due to which Fortran is better than C at some point. As far as I heard for now, Fortran is quite similar to C (in terms of performance and compiler optimizations efficiency as well) but C (and C++ especially) gives you wider possibilities, and because of this its hard to write non-efficient code in Fortran but it can be done easily in C.