Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- PARDISO - Divide by zero exception on second call

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted

Michael_R_7

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-14-2015
05:01 AM

5 Views

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 forINTEGER(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

Accepted Solutions

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

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-16-2015
08:52 AM

5 Views

3 Replies

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

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-16-2015
08:52 AM

6 Views

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

Michael_R_7

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-18-2015
04:41 AM

5 Views

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

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-18-2015
08:40 AM

5 Views

For more complete information about compiler optimizations, see our Optimization Notice.