Community
cancel
Showing results for 
Search instead for 
Did you mean: 
kong__yao
Beginner
127 Views

Question about matrix inverse

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)

0 Kudos
10 Replies
Gennady_F_Intel
Moderator
127 Views

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

kong__yao
Beginner
127 Views

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
Beginner
127 Views

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
Moderator
127 Views

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
Beginner
127 Views

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
Moderator
127 Views

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
Beginner
127 Views

I want to know which function is frequently used!

Gennady_F_Intel
Moderator
127 Views

?getri+?getrf

Spencer_P_Intel
Employee
127 Views

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
Beginner
127 Views

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

Reply