Community
cancel
Showing results for
Did you mean: Beginner
107 Views

## PARDISO - Divide by zero exception on second call

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

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

```

1 Solution Black Belt
107 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).

3 Replies Black Belt
108 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). Beginner
107 Views

Thanks mecej4, I did wonder if that was the issue, but would have expected an error returned rather than an exception thrown! Black Belt
107 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. 