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

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

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

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

Link Copied

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

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

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

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

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

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.

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