PARDISO - Divide by zero exception on second call

[Michael_R_7]

The following simple test produces a FPE in PARDISO.

The exception occurs when calling PARDISO for a second time, with a smaller matrix than the first; is this allowed? 

Note that to get the exception the Fortran project settings should be set to abort when dividing by zero (/fpe:0 I think)

Is this a bug in Pardiso, or am I using it incorrectly?

Many thanks!

    include'mkl_pardiso.f90'

    program TestPardiso
    USE MKL_PARDISO
    implicit none

    ! Variables
    TYPE(MKL_PARDISO_HANDLE) :: pt(64)     ! Pointer to the address of solver internal data - DO NOT EDIT!
    INTEGER(KIND=4) :: defaultparam(64)    ! Option Settings for 
    INTEGER(KIND=4) :: param(64)           ! Option Settings for 
    
    INTEGER :: maxfct, mnum, phase, nrhs, msglvl, errflag, i, perm(11), mtype, nnz
    
    INTEGER(KIND=4) :: N
    REAL(KIND=8), ALLOCATABLE  :: A(:)    
    INTEGER(KIND=4), ALLOCATABLE :: ia(:) 
    INTEGER(KIND=4), ALLOCATABLE :: ja(:) 
    REAL(KIND=8), ALLOCATABLE :: tmp1(:), tmp2(:)
     
    ! Initialize
    param = 0
    mtype = 2
    
    ! Setup default paramaters
    CALL pardisoinit(pt, mtype, defaultparam) ! Real and symmetric positive definite

    ! Set additional options
    defaultparam( 6) = 1 ! Solver stores the solution in the right-hand side b.
    defaultparam(27) = 1 ! Check input data
        
    ! Initialize output
    errflag = 0
           
    ! Setup input options
    maxfct = 1 ! Store one matrix at a time
    mnum = 1   ! Calculate with first matrix
    phase = 12 ! Analysis, numerical factorization
    perm = 0 ! Permutation (reordering) info
    nrhs = 1
    msglvl = 1 ! Output info to screen 
    
    ! Set Matrix For 1st run
    N = 6
    nnz = 8
    ALLOCATE(a(nnz), ia(N+1), ja(nnz), tmp1(N), tmp2(N))

    a(1:4) = (/ 6745627520.0D0, 241248804.696755D0, 2.980232238769531D-008, 241248804.696755D0 /)
    a(5:8) = (/-2.980232238769531D-008, 293365306.080000D0, 1607731120.02139D0, 1607731120.02139D0 /)
	    ia = (/ 1, 2, 4, 6, 7, 8, 9 /)
	    ja = (/ 1, 2, 6, 3, 5, 4, 5, 6 /)    
        
    ! Run Pardiso
    PRINT*, "Solving matrix One!"
    param = defaultparam ! Use default options
    CALL pardiso(pt, maxfct, mnum, mtype, phase, N, a, ia, ja, perm, nrhs, param, msglvl, tmp1, tmp2, errflag)

    ! Set Matrix for 2nd run (NNZ=6)
    nnz = 6
	    a(1:3) = (/ 6745627520.00000D0, 160301274.456755D0, 160301274.456755D0 /)
    a(4:6) = (/ 293365306.080000D0, 1382876869.35472D0, 1382876869.35472D0 /)
	    ia(:N+1) = (/ 1, 2, 3, 4, 5, 6, 7 /)
    ja = (/ 1, 2, 3, 4, 5, 6 /)
    
    ! Run Pardiso
    PRINT*, "Matrix 2!"
    param = defaultparam ! Use default options
    CALL pardiso(pt, maxfct, mnum, mtype, phase, N, a(:nnz), ia(:N+1), ja(:nnz), perm(:N), nrhs, param, msglvl, tmp1(:N), tmp2(:N), errflag)
        
    end program TestPardiso

    
    
    
    

[mecej4 · Accepted solution]

After factorizing the first matrix, you are changing the structure of the matrix (i.e., the number and locations of the nonzero entries). Therefore, you need to reset the array pt() to 0 before working on the second matrix, by calling Pardiso_Init again, or by setting pt = MKL_PARDISO_HANDLE(0). (It would also be good to do a "release memory" call to Pardiso before doing that).

[Michael_R_7]

Thanks mecej4, I did wonder if that was the issue, but would have expected an error returned rather than an exception thrown! 

[mecej4]

With complex software such as Pardiso, the design of the interface involves making weighing performance versus how much checking to do on the input arguments. Typically, checks are done on the first call only, since Pardiso is designed to allow multiple calls to different 'phases' of the linear equations solution problem.