Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
5 Views

PARDISO - Divide by zero exception on second call

Jump to solution

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

    
    
    
    

 

0 Kudos

Accepted Solutions
Highlighted
Black Belt
5 Views

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).

View solution in original post

0 Kudos
3 Replies
Highlighted
Black Belt
6 Views

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).

View solution in original post

0 Kudos
Highlighted
Beginner
5 Views

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

 

0 Kudos
Highlighted
Black Belt
5 Views

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.

0 Kudos