Modify MKL_PARDISO_PIVOT

[Javi_Dopazo]

Could someone help me by telling me how I can modify the function MKL_PARDISO_PIVOT

[AryanK_Intel]

Hi Javi_Dopazo,

Thanks for posting in Intel Communities.

For more information related to MKL_PARDISO_PIVOT, please refer to the below link:

https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/mkl-pardiso-pivot.html

You could also refer to the example provided under below path based on your usage and requirement.

C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\examples\examples_core_c.zip\c\sparse_directsolvers\source\pardiso_sym_getdiag.c (or)

C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\examples\examples_core_c.zip\c\sparse_directsolvers\source\pardiso_unsym_diag_pivot.c

In addition, Could you please elaborate on your issue, it helps in understanding your issue better.

Best regards,

Sri Raj Aryan.K

[Javi_Dopazo]

First of all I would like to thank your answer. The best solution for me it is the replacement of mkl_pardiso_pivot, but in fortran I can not. And the solution, if I am not wrong, it is the redefinition in C, and then mix with fortran. 

Do you think if it is correct my solution.

Javi

[AryanK_Intel]

Hi Javi_Dopazo,

Thanks for getting back. As requested by you, For more information related to MKL_PARDISO_PIVOT in fortran, you could go through to the below link:

https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2023-2/mkl-pardiso-pivot.html

You could also refer to the Fortran examples provided under the below path based on your use case.

C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\examples\examples_core_f.zip\f\sparse_directsolvers\source

Best regards,

Sri Raj Aryan.K

[AryanK_Intel]

Hi Javi_Dopazo,

A gentle reminder:

Has the information provided helped? Could you please let us know if you need any other details?

Best regards,

Sri Raj Aryan.K

[Javi_Dopazo]

Thanks for yours answer. I had checked the examples in C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\examples\examples_core_f.zip\f\sparse_directsolvers\source

And I could not find any example about rewrite the function mkl_pardiso_pivot or how it could be used.

Could you say to me where I can fin some exemples where mkl_pardiso_pivot is used.

Thank you,

Javi

[AryanK_Intel]

Hi Javi_Dopazo,

We are looking into your issue internally and we will update you soon.

Thanks and Regards,

Sri Raj Aryan.K

[Javi_Dopazo]

Thanks, and waiting for your answer.

Javi

[Javi_Dopazo]

I have made a simple example, which works for symmetric, indefinite but it do not for real and symmetric positive definite. And I don't know why.

program Console2

implicit none


!Ax = b

real(8):: A(5), b(5), x(6)
integer:: rows(6), columns(5), iparm(64), maxfct, mnum, mtype, phase, n, nrhs, error, msglvl, perm(1), i
integer:: pt(64)

pt = 0
n = 5
mnum = 1
maxfct = 1
nrhs = 1
error = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype = -2 ! symmetric, indefinite
!mtype = 2 ! real and symmetric positive definite

A = (/10.0d0, 20.0d0, 0.0000000001d0, 40.0d0, 50.0d0/)

rows = (/1,2,3,4,5,6/)

columns = (/1,2,3,4,5/)

b = (/10.0d0, 20.0d0, 30.0d0, 40.0d0, 50.0d0/)

x = 0.0d0

iparm = 0

!iparm(1) = 1 ! no solver default
!iparm(2) = 2 ! fill-in reordering from METIS
!iparm(4) = 0 ! no iterative-direct algorithm
!iparm(5) = 0 ! no user fill-in reducing permutation
!iparm(6) = 0 ! =0 solution on the first n components of x
!iparm(8) = 2 ! numbers of iterative refinement steps
!iparm(10) = 8 ! perturb the pivot elements with 1E-13
!iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
!iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
!iparm(14) = 0 ! Output: number of perturbed pivots
!iparm(18) = -1 ! Output: number of nonzeros in the factor LU
!iparm(19) = -1 ! Output: Mflops for LU factorization
!iparm(20) = 0 ! Output: Numbers of CG Iterations
!iparm(56) = 1 !You can use the mkl_pardiso_pivot callback routine to control pivot elements which appear during numerical factorization


phase = 11 ! only reordering and symbolic factorization
CALL pardiso(pt, maxfct, mnum, mtype, phase, n, A, rows, columns, perm, nrhs, iparm, msglvl, b, x, error)

write(*,*) 'phase = 11. Error: ',error

phase = 22 ! only factorization
CALL pardiso(pt, maxfct, mnum, mtype, phase, n, A, rows, columns, perm, nrhs, iparm, msglvl, b, x, error)

write(*,*) 'phase = 22. Error: ',error

if( error == -4) then

write(*,*) 'Pivot zero: ', iparm(30)

end if

phase = 33 ! only solution
CALL pardiso(pt, maxfct, mnum, mtype, phase, n, A, rows, columns, perm, nrhs, iparm, msglvl, b, x, error)

write(*,*) 'phase = 33. Error: ',error

do i = 1, 5

write(*,*) x(i)

end do

end program Console2

subroutine mkl_pardiso_pivot (ai, bi, eps)

real(8):: ai, bi, eps

if( bi < eps ) then

bi = 30.0d0

end if

end subroutine

[AryanK_Intel]

Hi Javi_Dopazo,

>> First of all I would like to thank your answer. The best solution for me it is the replacement of mkl_pardiso_pivot, but in fortran I can not. And the solution, if I am not wrong, it is the redefinition in C, and then mix with fortran.

Here is the code in C and Fortran of mkl_pardiso_pivot

In C:

int mkl_pardiso_pivot( const double*aii, double*bii, const double*eps )
{
    if ( (*bii > *eps) || ( *bii < -*eps ) )
        return 0;
    if ( *bii > 0 )
        *bii = *eps;
    else
        *bii = -*eps;
    return 1;
}

In Fortran:

FUNCTION mkl_pardiso_pivot(ai, bi, eps)
	DOUBLE, INTENT(IN) :: ai
	DOUBLE, INTENT(INOUT) :: bi
	DOUBLE, INTENT(IN) :: eps
	
	IF ( (bi > eps) .OR. (bi < (.NOT. eps)) ) THEN
		RETURN 0
	END IF
	
	IF (bi > 0) THEN
		bi = eps
	ELSE
		bi = (.NOT. eps)
	END IF
	RETURN 1
END FUNCTION mkl_pardiso_pivot

Could you please let us know if it is helpful?

Thanks and Regards,

Aryan.

[Javi_Dopazo]

I am sorry, but it isn't. The behavior is exactly like my example. It works for mytype = -2, but not for mytype = 2. And I am interested in mytype = 2. 

Can you help me?

Yours, 

Javi

[AryanK_Intel]

Hi Javi_Dopazo,

We are working on your issue internally and we will update you soon.

Thanks and Regards,

Sri Raj Aryan.K

[Javi_Dopazo]

any new Mr. Raj.

[Javi_Dopazo]

Waiting for it

Thank you.

Javi

[AryanK_Intel]

Hi Javi_Dopazo,

Apologies for the delay in response.

We would like to know regarding 

>>"it works for mtype = -2, but not for mtype = 2".

Do you mean that if MKL PARDISO detects zero, then MKL Pardiso doesn't stop execution and/or doesn't return an error = -4?

Thanks and Regards,

Aryan.

[Javi_Dopazo]

for type =2, and with a diagonal matriz, diferent zero numbers only in the main diagonal of matriz and zero the rest of matriz. If one number is negative or zero MKL_PARDISO return an error = -4. But if the number is positive and close to zero return error = 0, and if I tyr to use the diagonal and pivoting control  ( iparm(56) = 1 ), I can't. I have detected that the program doesn't acces to mkl_pardiso_pivot.

However the same problem but type = -2, I don't have this problem and I can control the pivoting using mkl_pardiso_pivot.

If you need more information, please say to me.

Thanks,

Javi

 

[AryanK_Intel]

Hi Javi_Dopazo,

For mtype=2 (i.e., real and symmetric positive definite), it is expected to have behavior as follows:

1. The MKL_PARDISO returns an error code of -4, if any number on the main diagonal is either negative or zero. 
2. The MKL_PARDISO returns an error code of 0, if the number on the main diagonal is positive and close to zero.

Notably, explicitly setting 0.0 instead of 0.000000001d, while mtype=2, then Pardiso returns -4 as expected.

Thanks & Regards,

Aryan

[Javi_Dopazo]

May I haven't explained correctily. For mytype = 2 and with iparm(56) = 1, the subroutine mkl_pardiso_pivot is not called. I attached a small program where you can see it.

Yours, 

Javi

[Javi_Dopazo]

Any new

Javi

[AryanK_Intel]

Hi,

Apologies for the wait in getting back to you.

We are working on your issue internally and we will provide you with an update shortly.

Thanks & Regards,

Aryan.

[Javi_Dopazo]

Thank you