Matrix determinant

[Beliavsky]

It would be convenient for stdlib to have a function for computing the determinant of a square matrix. The cases of symmetric positive definite (SPD) (use Cholesky decomposition) and a general matrix should be handled separately. Here is code for the SPD case. I have not included the choldc code from Numerical Recipes, but there are public domain routines for it.

function det_symm_pos_def(amat) result(det)
! compute the determinant of a positive definite symmetric matrixx using Cholesky decomposition
real(kind=dp)    , intent(in) :: amat(:,:) ! (n,n)
real(kind=dp)                 :: det
real(kind=dp)                 :: acp(size(amat,1),size(amat,2)),p(size(amat,1))
integer                       :: ierr
character (len=*), parameter  :: msg = "in det_symm_pos_def, "
if (size(amat,1) /= size(amat,2)) then
   write (*,*) msg // "shape(amat) =",shape(amat)," need dimensions to be the same, STOPPING"
   stop
end if
acp = amat
call choldc(acp,p,ierr) ! output: acp, p, ierr
if (ierr /= 0) then
   write (*,*) msg // "returning from choldc, ierr =",ierr," STOPPING"
   stop
end if
det = product(p**2)
end function det_symm_pos_def

[awvwgk]

Sounds like something that could be done via LAPACK's ?potrf functionality. Not sure if we really want to implement such functionality ourselves and try to compete with optimized implementations like MKL.

[ivan-pi]

Actually, matrix determinant is something we likely have to implement ourselves on top of the factorizations from LAPACK (product of the factors in the diagonal). See e.g. the code used in numpy: https://github.com/numpy/numpy/blob/fb215c76967739268de71aa4bda55dd1b062bc2e/numpy/linalg/umath_linalg.c.src#L1036

Unfortunately, convenience linear algebra wrappers are kind of a daunting area for many reasons:

• involve a great amount of duplication to handle different types and deal with idiosyncrasies of the underlying LAPACK library
• little agreement on what error reporting should look like in the wrapper routines
• matrix shape and other runtime assertions