How to test if matrix is diagonal?

Question

I need to test if one variance matrix is diagonal. If not, I'll do Cholesky LDL decomposition. But I was wondering which would be most reliable and fastest way to test is matrix diagonal? I'm using Fortran.

First thing that come to my mind is to take sum of all elements of matrix, and substract diagonal elements from that sum. If the answer is 0, matrix is diagonal. Any better ideas?

In Fortran I'll write

!A is my matrix
k=0.0d0
do i in 1:n #n is the number of rows/colums
k = k + A(i,i)
end do

if(abs(sum(A)-k) < epsilon(k)*sum(A)) then
#do cholesky LDL, which I have to write myself, haven't found any subroutines for that  in Lapack or anywhere else
end if

Answer 1

It would be much better to just traverse all the off-diagonal elements and test if they are near zero (comparing a floating-point number for inequality is prone to rounding errors and can lead to erroneous results).

First, once you find any violating element you can immediately stop traversing and this may allow for significant time decrease if violating matrices are typical.

Second, it would potentially allow for better loop unrolling by the compiler (Fortran compilers are known for good optimization strategies) and for faster on-chip execution due to less inter-instruction dependencies.

Add to this the fact that your suggested algorithm is prone to overflows and error accumulation and the "traverse-and-test" algorithm is not.

Answer 2

Search the matrix for non zero values

logical :: not_diag
integer :: i, j

not_diag = .false.

outer: do i = 2, size(A,1)
  do j = i, size(A, 2)
    if (A(i,j) > PRECISION) then
      not_diag = .true.
      exit outer
    end if
  end
end outer

if (not_diag) then
  ! DO LDL' decomposition
end if

To use the double precision LAPACK routines replace the first 's' with 'd'. So spotrf becomes dpotrf

http://www.netlib.org/lapack/double/