¿Cómo probar si la matriz es diagonal?

Question

Necesito comprobar si una matriz de varianza es diagonal. Si no, haré la descomposición de Cholesky LDL. Pero me preguntaba cuál sería la forma más confiable y más rápida de probar es la matriz diagonal. Estoy usando Fortran.

Lo primero que me viene a la mente es tomar suma de todos los elementos de la matriz y restar elementos diagonales de esa suma. Si la respuesta es 0, la matriz es diagonal. Alguna mejor idea?

En Fortran Voy a escribir

!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

Buscar la matriz de valores cero no

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 

Para utilizar el doble precisión rutinas LAPACK reemplazan a la primera 's' con 'd'. Así spotrf convierte dpotrf

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

Answer 2

Sería mucho mejor simplemente recorrer todos los elementos fuera de la diagonal y prueba si están cerca de cero (comparando un número de coma flotante para la desigualdad es propenso a errores de redondeo y pueden conducir a resultados erróneos).

En primer lugar, una vez que encuentre un elemento violador, puede dejar de atravesarlo inmediatamente y esto puede permitir una disminución significativa del tiempo si las matrices violatorias son típicas.

En segundo lugar, potencialmente permitiría un mejor desenrollado del bucle por el compilador (los compiladores de Fortran son conocidos por sus buenas estrategias de optimización) y por una ejecución más rápida en el chip debido a menos dependencias entre instrucciones.

Agregue a esto el hecho de que su algoritmo sugerido es propenso a desbordamientos y acumulación de errores y el algoritmo "transversal y de prueba" no lo es.