topic Progressively increasing error in calculation in Intel® oneAPI Math Kernel Library

Progressively increasing error in calculation

S__MPay — Thu, 27 Apr 2017 15:23:46 GMT

This function calculates determinant of a matrix:

real*8 function mklDet(A)
    real*8, intent(in), dimension(:,:) :: A         !< input matrix A
    real*8, dimension(size(A,1)) :: ipiv            !! pivot indices
    integer N, info, i                              !> dimension, status
    
    N = size(A, 1)
 
    ! DGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call DGETRF(n, n, A, n, ipiv, info)
    
    !< This part of code by adopted from intel software forum topic 309460
    mklDet = 0.d0
    if (info > 0) then
        return
    else if (info < 0) then
        print *, 'Fatal error calculating determinant!'
    end if
    mklDet = 1.d0
    do i = 1, n
        if (ipiv(i).ne.i) then
            mklDet = -mklDet * a(i,i)
        else
            mklDet = mklDet * a(i,i)
        endif
    end do
    !>
end function mklDet

Now if we try to calculate determinant of a matrix with which should be zero:

!! [[6.d10,  3.d10,  5.d10];
!!  [4.d10,  4.d10,  7.d10];
!!  [12.d10, 6.d10, 10.d10]];

print *, 'zero determinant' , mklDet(reshape((/ 6.d10, 3.d10, 5.d10, 4.d10, 4.d10, 7.d10, 12.d10, 6.d10, 10.d10 /),(/3, 3/)))

as you can see the outcome is a very large number and not zero:

 zero determinant -5.329070518200751E+015

but if we use smaller numbers in the matrix

!! [[6.,  3.,  5.];
!!  [4.,  4.,  7.];
!!  [12., 6., 10.]];

print *, 'zero determinant' , mklDet(reshape((/ 6.d0, 3.d0, 5.d0, 4.d0, 4.d0, 7.d0, 12.d0, 6.d0, 10.d0 /),(/3, 3/)))

The outcome is much closer to zero:

 zero determinant  1.065814103640150E-014

How do I improve the accuracy of my calculations?

Update: When I try same

S__MPay — Mon, 01 May 2017 16:33:56 GMT

Update: When I try same function with quad precision the error is exactly zero:

    print *, 'zero determinant double precision' , dblMklDet(reshape((/ 6.d0, 3.d0, 5.d0, 4.d0, 4.d0, 7.d0, 12.d0, 6.d0, 10.d0 /),(/3, 3/)))
    print *, 'zero determinant  quad  precision' , qprMklDet(reshape((/ 6.q0, 3.q0, 5.q0, 4.q0, 4.q0, 7.q0, 12.q0, 6.q0, 10.q0 /),(/3, 3/)))

The output is:

 zero determinant double precision  1.065814103640150E-014
 zero determinant  quad  precision  0.000000000000000000000000000000000E+0000

I don't understand what has caused the problem in double precision.

Here is the source for quad precision version of my MklDet function:

      real(qp) function qprMklDet(A)
          real(qp), intent(in), dimension(:,:) :: A         !< input matrix A
          real(qp), dimension(size(A,1)) :: ipiv            !! pivot indices
          integer N, info, i                                !> dimension, status
          
          N = size(A, 1)
      
          ! DGETRF computes an LU factorization of a general M-by-N matrix A
          ! using partial pivoting with row interchanges.
          call DGETRF(n, n, A, n, ipiv, info)
          
          !< This part of code by adopted from intel software forum topic 309460
          qprMklDet = 0.d0
          if (info > 0) then
              return
          else if (info < 0) then
              print *, 'Fatal error calculating determinant!'
          end if
          qprMklDet = 1.d0
          do i = 1, n
              if (ipiv(i).ne.i) then
                  qprMklDet = -qprMklDet * a(i,i)
              else
                  qprMklDet = qprMklDet * a(i,i)
              endif
          end do
          !>
      end function qprMklDet

And qp is

      use, intrinsic :: iso_fortran_env, only: REAL64, REAL128
      implicit none
      !==========================================
      ! C  o  n  s  t  a  n  t  s
      !==========================================
      !integer, parameter :: sp = REAL32
      integer, parameter :: dp = REAL64
      integer, parameter :: qp = REAL128
      !------------------------------------------

Hi,

Konstantin_A_Intel — Mon, 08 May 2017 16:03:00 GMT

Hi,

Det calculation using LU is not stable numerically. E.g., it's clearly described at matlab page of Det function here:

https://www.mathworks.com/help/matlab/ref/det.html

So, please just avoid calling LU with singular matrices. You can first use other functions like ?gecon in order to estimate condition number of the matrix.

Regards,

Konstantin