Question about matrix inverse

kong__yao · ‎02-20-2019

Recetnly, i was trying to replace the function for matrix inverse by that of MKL. However, after making it, i found that the time used for computing was not improved and be slower. I don't know the real reason.It could be appreciated that anyone can help me.

The original subroutine was syminvg.

SUBROUTINE syminvg(n,q,norm,nsing,parflg,sbrName,bigajj_limit,utest)

! -------------------------------------------------------------------------
! Purpose: Inversion symm. matrix with pivot search on the diagonal.
! Algorithm acm 150, H. Rutishauser 1963
! Change: upper triangle of matrix in vektor q(n*(n+1)/2)
!
! Author: H. Rutishauser, C. Urschl
!
! Created: 24-Jan-2003
!
! Changes: 09-Aug-2010 RD: Generic use, only one version for all purposes
! 27-Oct-2010 SL: use m_bern with ONLY
! 23-May-2011 RD: Special handling for diagonal elements with NaN
! 30-Jul-2012 RD: Evaluate the inversion (only for GRP_AIUB)
! 22-Sep-2012 RD: Change variable name "unit" to "utest"
!
! Copyright: Astronomical Institute
! University of Bern
! Switzerland
! -------------------------------------------------------------------------

I used two subroutines of intel mkl:

CALL DPPTRF('U',n,q0,INFO) CALL DPPTRI('U',n,q0,INFO)

Gennady_F_Intel · ‎02-20-2019

This is not too frequently used routines from mkl. What performance difference do you see?

kong__yao · ‎02-20-2019

The performance of Intel MKL is not better than original code. Could you tell me which routines is frequently used for matrix inversion?

kong__yao · ‎02-20-2019

The code using intel mkl was as follows:

MODULE s_SYMINVMKL
CONTAINS

! -------------------------------------------------------------------------
! Bernese GPS Software Version 5.1
! -------------------------------------------------------------------------

SUBROUTINE syminvmkl(n,q,norm,nsing,parflg,sbrName,bigajj_limit,utest)

! -------------------------------------------------------------------------
! Purpose: Inversion symm. matrix using subroutine using Intel MKL
! Change: upper triangle of matrix in vektor q(n*(n+1)/2)
!
! Author: Yao Kong, Baoqi Sun
!
! Created: 24-Jan-2019
!
! Changes: 09-Aug-2010 RD: Generic use, only one version for all purposes
! 27-Oct-2010 SL: use m_bern with ONLY
! 23-May-2011 RD: Special handling for diagonal elements with NaN
! 30-Jul-2012 RD: Evaluate the inversion (only for GRP_AIUB)
! 22-Sep-2012 RD: Change variable name "unit" to "utest"
!
! Copyright: Astronomical Institute
! University of Bern
! Switzerland
! -------------------------------------------------------------------------

USE m_bern, ONLY: i4b, r8b, lfnErr

USE s_matcop
USE s_statis
USE f_ikf

IMPLICIT NONE

! List of Parameters
! ------------------
INTEGER(i4b),INTENT(IN) :: n ! Number of parameters
REAL(r8b),DIMENSION(n*(n+1)/2) :: q ! Upper trianle, columnwise
INTEGER(i4b) :: norm ! 1: normalize the diagonale of q
INTEGER(i4b),INTENT(OUT) :: nsing ! Number of singular parameters
INTEGER(i4b),DIMENSION(n), &
OPTIONAL :: parflg ! Flag for singular parameters
! =0: parameter not singular
! =1: parameter singular
CHARACTER(LEN=*), OPTIONAL :: sbrName ! Name of calling routine to enable
! reporting of singularities
REAL(r8b), OPTIONAL :: bigajj_limit ! limit for singulatiy detection
CHARACTER(LEN=*), OPTIONAL :: utest ! Name of calling routine to check
! the numerical performance

! Local Variables
! ---------------
CHARACTER(LEN=1), DIMENSION(n) :: ls1

REAL(r8b), DIMENSION(n) :: hs1
REAL(r8b), DIMENSION(n) :: hs2
REAL(r8b), DIMENSION(n) :: scl_tot
REAL(r8b) :: xRms1,xRms0
REAL(r8b) :: xMin1,xMin0
REAL(r8b) :: xMax1,xMax0
REAL(r8b), POINTER, &
DIMENSION(:) :: q0
REAL(r8b) :: ajj
REAL(r8b) :: bigajj
REAL(r8b) :: limit = 1.0d-10 ! limit to detect a singular
! parameter

INTEGER(i4b) :: i, j, k
INTEGER(i4b) :: ij, jj, jk, kj, kk
INTEGER(i4b) :: kp1, km1
INTEGER(i4b) :: jkmax
INTEGER(i4b) :: ising
INTEGER(i4b) :: INFO
INTEGER(i4b), DIMENSION(n) :: IPIV
REAL(r8b), DIMENSION(n) :: WORK

LOGICAL :: prtFlg
LOGICAL :: isSing
LOGICAL :: doNorm = .TRUE.

! Overwrite the limitation of singular parameters
! if the optional parameter is set
IF (PRESENT(bigajj_limit)) THEN
limit = bigajj_limit
ENDIF

! Write or not write warning message
prtFlg = (PRESENT(sbrName))

! Make a copy of the input matrix
! -------------------------------
NULLIFY(q0)
#ifdef GRP_AIUB
IF (PRESENT(utest)) THEN
!! Print matrix "q" to screen formatted for MATLBA
!! write(*,'(A)') '>>q0<< q = [ ...'
!! DO i = 1,n
!! DO j = 1,n-1
!! WRITE(*,'(A,E30.20,A)') '>>q0<< ',q(IKF(i,j)),' ...'
!! ENDDO
!! IF ( i < n ) THEN
!! WRITE(*,'(A,E30.20)') '>>q0<< ',q(IKF(i,j))
!! ELSE
!! WRITE(*,'(A,E30.20,A)') '>>q0<< ',q(IKF(i,j)),' ]'
!! ENDIF
!! ENDDO
!! End of printing section
CALL matcop(1,IKF(n,n),q,q0)
ENDIF
#endif

CALL matcop(1,IKF(n,n),q,q0)
CALL DPPTRF('U',n,q0,INFO)
IF (INFO.EQ.0)THEN
CALL DPPTRI('U',n,q0,INFO)
q=q0
ELSE
!WRITE(*,*)"Matrix is not positive-definite"
CALL DSPTRF('U',n,q,IPIV,INFO)
CALL DSPTRI('U',n,q,IPIV,WORK,INFO)
END IF

Gennady_F_Intel · ‎02-21-2019

I don't understand why you attached and added the the same code to this thread at the same time. Attachment will be enough, but it is not still clear what is the performance differences do you see between mkl's and syminvg?

kong__yao · ‎02-21-2019

I am sorry for adding the same code. The performance difference is nearly zero. For a matrix A(753,753), to get the inverse of A, the time used by two function is nealy the same.

Gennady_F_Intel · ‎02-21-2019

Ok, the problem size is not too big. We expect to see the performance would be better with large problem sizes. You may try to check. From our side, we never tried this lib in the past and will be useful to know the results.

kong__yao · ‎02-22-2019

I want to know which function is frequently used!

Gennady_F_Intel · ‎02-22-2019

?getri+?getrf

Spencer_P_Intel · ‎02-22-2019

To elaborate on what Gennady has pointed out. It is very rare (but not unheard of) to have an application that truly needs the full inverse of a matrix. This is actually a rather unstable problem. However, factoring the matrix in some manner like Cholesky (A = LL^T) or the LU factorization (A=LU) is a much more stable problem. Then it is just a couple of applications of forward or backward solve to get the action of A inverse without actually computing the inverse. From, a numerical standpoint, this is much more efficient and stable. Now there are some cases where having the actual inverse is actually important and without knowing your actual problem statement, we cannot advise on whether the ?getri/?getrf functions are sufficient, but with fairly high probability they are the right tools for the job.

kong__yao · ‎03-06-2019

I used dgetri and dgetrif functions，however the performance is stiil not improved!