Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## Progressively increasing error in calculation

Beginner
279 Views

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?

2 Replies
Beginner
279 Views

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

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
!------------------------------------------
```

Employee
279 Views

Hi,

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

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